MATLAB课件10

17
n=64; [p,e,t]=poimesh(g,n); tm=cputime; u=poisolv(b,p,e,t,f); cputime-tm tm=cputime; u1=assempde(b,p,e,t,c,a,f); cputime-tm
18
1
exact=(1-p(1,:).^2-p(2,:).^2)’/4; … er=norm(u-exact,’inf’); … error=[error er]; … fprintf(‘Error:%e.Number of nodes:%d\n’,er,size(p,2)); … end pdesurf(p,t,u); pdesurf(p,t,u-exact); 11.2 解Helmholtz方程并研究反射波 11.2.1 Helmholtz方程的求解 k=60; g=‘scatterg’; b=‘scatterb’;
11
umax=max(max(u1)); umin=min(min(u1)); For j=1:nframes,… u=tri2grid(p,t,u1(:,j),tn,a2,a3); i=find(isnan(u)); u(i)=zeros(size(i);… surf(x,y,u); caxis([umin umax]); colormap(cool),… axis([-1 1 –1 1 0 1]);… mv(:,j)=getframe;… end Movie(mv,10) 11.6 求解波形传递问题
2
c=1; a=k^2; f=0; [p,e,t]=initmesh(g); [p,e,t]=refinemesh(g,p,e,t); [p,e,t]=refinemesh(g,p,e,t); pdemesh(p,e,t); u=assempde(b,p,e,t,c,a,f); h=newplot;set(get(h,’Parent’),’Renderer’,’zbuffer’)
15
x=p(1,:)’; y=p(2,:)’; R=sqrt(x.^2+y.^2); uu=-log®/2/pi; pdeplot(p,e,t,’xydata’,u-uu,’zdata’,u-uu,’mesh’,’off’); 11.8 使用矩形栅格解泊松方程
u 3x
2
2
g=‘squareg’; b=‘squareb4’; c=1
C(ic2,ic2)=C(ic2,ic2)+K(ic2,ic2)-a2*a2’; f2=F(i2);e2=K2’\f2;FC(ic2)=FC(ic2)+F(c2)-a2*e2; [i3,c3]= pdesdp(p,e,t,3);ic3=pdesubix(cp,c3); [K,F]=assempde(b,p,e,t,c,a,f,time,3); K2=K(i3,i3);d=symmmd(K3);i2=i3(d); K3=chol(K3(d,d)); B3=K(c3,i3);a3=B3/K3; C(ic3,ic3)=C(ic3,ic3)+K(ic3,ic3)-a3*a3’; f3=F(i3);e3=K3’\f3; FC(ic3)=FC(ic3)+F(c3)-a3*e3; u=zeros(np,1); u(cp)=C\FC; u(i1)=K1\(e1-a1’*u(c1));
14
11.7 点力和适应解
u ( x, y) g=‘circleg’;
2
b=‘circleb1’; c=1; a=0; f=‘circlef’; [u,p,e,t]=adaptmesh(g,b,c,a,f,’tripick’,’circlepick’,’ma xt’,2000,’par’,1e-3); pdemesh(p,e,t) pdeplot(p,e,t,’xydata’,u,’zdata’,u,’mesh’,’off’);
5
movie(hf,M,50); 11.3 最小表面问题求解 g=‘circleg’; b=‘circleb2’; c=‘1./sqrt(1+ux.^2+uy.^2)’; a=0;f=0;rtol=1e-3; [p,e,t]=initmesh(g); [p,e,t]=refinemesh(g,p,e,t); u=pdenonlin(b,p,e,t,c,a,f,’tol’,rtol); pdesurf(p,t,u);
12
u u u u ( 2 2 2 ) 0 2 t x y z
2 2 2 2
g=‘squareg’; b=‘squareb3’;
c=1;
a=0;
f=0;
d=1;
[p,e,t]=initmesh(‘squareg’);
13
x=p(1,:)’; y=p(2,:)’; u0=atan(cos(pi/2*x)); ut0=3*sin(pi*x).exp(sin(pi/2*y)); n=31; tlist=linspace(0,5,n); uu=hyperbolic(u0,ut0,tlist,b,p,e,t,c,af,d); delta=-1:0.1:1; [uxy,tn,a2,a3]=tri2grid(p,t,uu(:,1),delta,delta); gp=[tn,a2;a3];
9
u(i2)=K2\(e2-a2’*u(c2)); u(i3)=K3\(e3-a3’*u(c3)); pdesurf(p,t,u) 11.5 求解热传导方程
u u u u ( 2 2 2 ) 0 t x y z g=‘squareg’;
2 2 2
b=‘squareb1’; c=1; a=0; f=1; d=1;
fptint(‘%d’,j); … Pdeplot(p,e,t,’xydata’,uu,’colorbar’,’off’,’mesh’,’off ’), … caxis([-maxu maxu]); … axis tight,set(gca,’DataAspectRatio’,[1 1 1]);axis off, … M(:,j)=getframe(hf); … if j==m, … fprintf(‘done\n’); … end,… end
第11章 MATLAB在偏微分方程 解法中的应用
11.1 解简单泊松方程 g=‘circleg’; b=‘circleb1’; c=1;a=0;f=1; [p,e,t]=initmesh(g,’hmax’,1); error=[];er=1; while er>0.001, … [p,e,t]=refinemesh(g,p,e,t); … u=assempde(b,p,e,t,c,a,f); …
16
a=0; f=‘3*x.^2’; n=16; [p,e,t]=poimesh(g,n); pdemesh(p,e,t);axis equal tm=cputime; u=poisolv(b,p,e,t,f); cputime-tm tm=cputime; u1=assempde(b,p,e,t,c,a,f); cputime-tm
10
Hale Waihona Puke [p,e,t]=initmesh(g); u0=zeros(size(p,2),1); ix=find(sqrt(p(1,:).^2+p(2,:).^2)<0.4); u0(ix)=ones(size(ix)); nframes=20; tlist=linspace(0,0.1,nframes); u1=parabolic(u0,tlist,b,p,e,t,c,a,f,d); X=linspace(-1,1,31); y=x; [unused,tn,a2,a3]=tri2grid(p,t,u0,x,y); newplot; Mv=moviein(nframes);
pdeplot(p,e,t,’xydata’,real(u),’zdata’,real(u),’mesh’,’o 3 f’);
colormap(cool) 11.2.2 反射波的可视化研究 movie.m m=10; h=newplot;hf=get(h,’Parent’); set(hf,’Renderer’,’zbuffer’) axis tight,set(gca,’DataAspectRatio’,[1 1 1]);axis off M=moviein(m,hf); maxu=max(abs(u)); for j=1;m, … uu=real(exp(-j*2*pi/m*sqrt(-1))*u); … 4
11.4 使用子区域分解FEM问题
7
C=zeros(nc,nc); FC=zeros(nc,1); [i1,c1]=pdesdp(p,e,t,1);ic1=pdesubix(cp,c1); [K,F]=assempde(b,p,e,t,c,a,f,time,1); K1=chol(k1(d,d));B1=K(c1,i1);a1=B1/K1; C(ic1,ic1)=C(ic1,ic1)+K(ic1,ic1)-a1*a1’; f1=F(i1);e1=K1’\f1;FC(ic1)=FC(ic1)+F(c1)-a1*e1; [i2,c2]=pdesdp(p,e,t,2);ic2=pdesubix(cp,c2); [K,F]=assempde(b,p,e,t,c,a,f,time,2); K2=K(i2,i2);d=symmmd(K2);i2=i2(d); 8 K2=chol(K2(d,d));B2=K(c2,i2);a2=B2/K2;
6
• • • • • • • • • • •
2u 1 L形求解且边界上u=0; g= ‘lshapeg’; b=‘lshapeg’; c=1;a=0;f=1; time= []; [p,e,t]=initmesh(g); [p,e,t] =refinemesh(g,p,e,t); [p,e,t]=refinemesh(g,p,e,t); np=size(p,2); cp=pdesdp(p,e,t); nc=length(cp);