偏微分方程的数值解法.docx
- 文档编号:3247052
- 上传时间:2022-11-21
- 格式:DOCX
- 页数:18
- 大小:16.61KB
偏微分方程的数值解法.docx
《偏微分方程的数值解法.docx》由会员分享,可在线阅读,更多相关《偏微分方程的数值解法.docx(18页珍藏版)》请在冰豆网上搜索。
偏微分方程的数值解法
偏微分方程的数值解法
1.peEllip5用五点差分格式解拉普拉斯方程
functionu=peEllip5(nx,minx,maxx,ny,miny,maxy)
formatlong;
hx=(maxx-minx)/(nx-1);
hy=(maxy-miny)/(ny-1);
u0=zeros(nx,ny);
forj=1:
ny
u0(j,1)=EllIni2Uxl(minx,miny+(j-1)*hy);
u0(j,nx)=EllIni2Uxr(maxx,miny+(j-1)*hy);
end
forj=1:
nx
u0(1,j)=EllIni2Uyl(minx+(j-1)*hx,miny);
u0(ny,j)=EllIni2Uyr(minx+(j-1)*hx,maxy);
end
A=-4*eye((nx-2)*(ny-2),(nx-2)*(ny-2));
b=zeros((nx-2)*(ny-2),1);
fori=1:
(nx-2)*(ny-2)
ifmod(i,nx-2)==1
ifi==1
A(1,2)=1;
A(1,nx-1)=1;
b
(1)=-u0(1,2)-u0(2,1);
else
ifi==(ny-3)*(nx-2)+1
A(i,i+1)=1;
A(i,i-nx+2)=1;
b(i)=-u0(ny-1,1)-u0(ny,2);
else
A(i,i+1)=1;
A(i,i-nx+2)=1;
A(i,i+nx-2)=1;
b(i)=-u0(floor(i/(nx-2))+2,1);
end
end
else
ifmod(i,nx-2)==0
ifi==nx-2
A(i,i-1)=1;
A(i,i+nx-2)=1;
b(i)=-u0(1,nx-1)-u0(2,nx);
else
ifi==(ny-2)*(nx-2)
A(i,i-1)=1;
A(i,i-nx+2)=1;
b(i)=-u0(ny-1,nx)-u0(ny,nx-1);
else
A(i,i-1)=1;
A(i,i-nx+2)=1;
A(i,i+nx-2)=1;
b(i)=-u0(floor(i/(nx-2))+1,nx);
end
end
else
ifi>1&&i A(i,i-1)=1; A(i,i+nx-2)=1; A(i,i+1)=1; b(i)=-u0(1,i+1); else ifi>(ny-3)*(nx-2)&&i<(ny-2)*(nx-2) A(i,i-1)=1; A(i,i-nx+2)=1; A(i,i+1)=1; b(i)=-u0(ny,mod(i,(nx-2))+1); else A(i,i-1)=1; A(i,i+1)=1; A(i,i+nx-2)=1; A(i,i-nx+2)=1; end end end end end ul=A\b; fori=1: (ny-2) forj=1: (nx-2) u(i,j)=ul((i-1)*(nx-2)+j); end end formatshort; 2.peEllip5m用工字型差分格式解拉普拉斯方程 functionu=peEllip5m(nx,minx,maxx,ny,miny,maxy) formatlong; hx=(maxx-minx)/(nx-1); hy=(maxy-miny)/(ny-1); u0=zeros(nx,ny); forj=1: ny u0(j,1)=EllIni2Uxl(minx,miny+(j-1)*hy); u0(j,nx)=EllIni2Uxr(maxx,miny+(j-1)*hy); end forj=1: nx u0(1,j)=EllIni2Uyl(minx+(j-1)*hx,miny); u0(ny,j)=EllIni2Uyr(minx+(j-1)*hx,maxy); end A=-4*eye((nx-2)*(ny-2),(nx-2)*(ny-2)); b=zeros((nx-2)*(ny-2),1); fori=1: (nx-2)*(ny-2) ifmod(i,nx-2)==1 ifi==1 A(1,nx)=1; b (1)=-u0(1,1)-u0(3,1)-u0(1,3); else ifi==(ny-3)*(nx-2)+1 A(i,i-nx+1)=1; b(i)=-u0(ny,1)-u0(ny,3)-u0(ny-2,1); else A(i,i-nx+3)=1; A(i,i+nx-1)=1; b(i)=-u0(floor(i/(nx-2))+1,1)-u0(floor(i/(nx-2))+3,1); end end else ifmod(i,nx-2)==0 ifi==nx-2 A(i,i+nx-1)=1; b(i)=-u0(1,nx-2)-u0(1,nx)-u0(3,nx); else ifi==(ny-2)*(nx-2) A(i,i-nx+1)=1; b(i)=-u0(ny,nx)-u0(ny,nx-2)-u0(ny-2,nx); else A(i,i-nx+1)=1; A(i,i+nx-3)=1; b(i)=-u0(floor(i/(nx-2)),nx)-u0(floor(i/(nx-2))+2,nx); end end else ifi>1&&i A(i,i+nx-1)=1; A(i,i+nx-3)=1; b(i)=-u0(1,i)-u0(1,i+2); else ifi>(ny-3)*(nx-2)&&i<(ny-2)*(nx-2) A(i,i-nx+3)=1; A(i,i-nx+1)=1; b(i)=-u0(ny,mod(i,(nx-2)))-u0(ny,mod(i,(nx-2))+2); else A(i,i-nx+3)=1; A(i,i+nx-1)=1; A(i,i+nx-3)=1; A(i,i-nx+1)=1; end end end end end ul=A\b; fori=1: (ny-2) forj=1: (nx-2) u(i,j)=ul((i-1)*(nx-2)+j); end end formatshort; 3.peHypbYF用迎风格式解对流方程 functionu=peYF(a,dt,n,minx,maxx,M) formatlong; h=(maxx-minx)/(n-1); ifa>0 forj=1: (n+M) u0(j)=IniU(minx+(j-M-1)*h); end else forj=1: (n+M) u0(j)=IniU(minx+(j-1)*h); end end u1=u0; fork=1: M ifa>0 fori=(k+1): n+M u1(i)=-dt*a*(u0(i)-u0(i-1))/h+u0(i); end else fori=1: n+M-k u1(i)=-dt*a*(u0(i+1)-u0(i))/h+u0(i); end end u0=u1; end ifa>0 u=u1((M+1): M+n); else u=u1(1: n); end formatlong; 4.peHypbLax用拉克斯-弗里德里希斯格式解对流方程 functionu=peHypbLax(a,dt,n,minx,maxx,M) formatlong; h=(maxx-minx)/(n-1); forj=1: (n+2*M) u0(j)=IniU(minx+(j-M-1)*h); end u1=u0; fork=1: M fori=k+1: n+2*M-k u1(i)=-dt*a*(u0(i+1)-u0(i-1))/h/2+(u0(i+1)+u0(i-1))/2; end u0=u1; end u=u1((M+1): (M+n)); formatshort; 4.peHypbLaxW用拉克斯-温德洛夫格式解对流方程 functionu=peLaxW(a,dt,n,minx,maxx,M) formatlong; h=(maxx-minx)/(n-1); forj=1: (n+2*M) u0(j)=IniU(minx+(j-M-1)*h); end u1=u0; fork=1: M fori=k+1: n+2*M-k u1(i)=dt*dt*a*a*(u0(i+1)-2*u0(i)+u0(i-1))/2/h/h-... dt*a*(u0(i+1)-u0(i-1))/h/2+u0(i); end u0=u1; end u=u1((M+1): (M+n)); formatshort; 6.peHypbBW用比姆-沃明格式解对流方程 functionu=peBW(a,dt,n,minx,maxx,M) formatlong; h=(maxx-minx)/(n-1); forj=1: (n+2*M) u0(j)=IniU(minx+(j-2*M-1)*h); end u1=u0; fork=1: M fori=2*k+1: n+2*M u1(i)=u0(i)-dt*a*(u0(i)-u0(i-1))/h-a*dt*(1-a*dt/h)*... (u0(i)-2*u0(i-1)+u0(i-2))/2/h; end u0=u1; end u=u1((2*M+1): (2*M+n)); formatshort; 7.peHypbRich用Richtmyer多步格式解对流方程 functionu=peRich(a,dt,n,minx,maxx,M) formatlong; h=(maxx-minx)/(n-1); forj=1: (n+4*M) u0(j)=IniU(minx+(j-2*M-1)*h); end u1=u0; fork=1: M fori=2*k+1: n+4*M-2*k tmpU1=-dt*a*(u0(i+2)-u0(i))/h/4+(u0(i+2)+u0(i))/2; tmpU2=-dt*a*(u0(i)-u0(i-2))/h/4
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分方程 数值 解法