偏微分方程的数值解法.docx
- 文档编号:16185192
- 上传时间:2023-07-11
- 格式:DOCX
- 页数:21
- 大小:16.55KB
偏微分方程的数值解法.docx
《偏微分方程的数值解法.docx》由会员分享,可在线阅读,更多相关《偏微分方程的数值解法.docx(21页珍藏版)》请在冰点文库上搜索。
偏微分方程的数值解法
偏微分方程的数值解法
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+(u0(i)+u0(i-2))/2; u1(i)=-dt*a*(tmpU1-tmpU2)/h/2+u0(i); end u0=u1; end u=u1((2*M+1): (2*M+n)); formatshort; 8.peHypbMLW用拉克斯-温德洛夫多步格式解对流方程 functionu=peMLW(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 tmpU1=-dt*a*(u0(i+1)-u0(i))/h/2+(u0(i+1)+u0(i))/2; tmpU2=-dt*a*(u0(i)-u0(i-1))/h/2+(u0(i)+u0(i-1))/2; u1(i)=-dt*a*(tmpU1-tmpU2)/h+u0(i); end u0=u1; end u=u1((M+1): (M+n)); formatshort; 9.peHypbMC用MacCormack多步格式解对流方程 functionu=peMC(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 tmpU1=-dt*a*(u0(i+1)-u0(i))/h+u0(i); tmpU2=-dt*a*(u0(i)-u0(i-1))/h+u0(i-1); u1(i)=-dt*a*(tmpU1-tmpU2)/h/2+(u0(i)+tmpU1)/2; end u0=u1; end u=u1((M+1): (M+n)); formatshort; 10.peHypb2LF用拉克斯-弗里德里希斯格式解二维对流方程的初值问题 functionu=pe2LF(a,b,dt,nx,minx,maxx,ny,miny,maxy,M) %啦-佛 formatlong; hx=(maxx-minx)/(nx-1); hy=(maxy-miny)/(ny-1); fori=1: nx+2*M forj=1: (ny+2*M) u0(i,j)=Ini2U(minx+(i-M-1)*hx,miny+(j-M-1)*hy); end end u1=u0; fork=1: M fori=k+1: nx+2*M-k forj=k+1: ny+2*M-k u1(i,j)=(u0(i+1,j)+u0(i-1,j)+u0(i,j+1)+u0(i,j-1))/4... -a*dt*(u0(i+1,j)-u0(i-1,j))/2/hx... -b*dt*(u0(i,j+1)-u0(i,j-1))/2/hy; end end u0=u1; end u=u1((M+1): (M+nx),(M+1): (M+ny)); formatshort; 11.peHypb2FL用拉克斯-弗里德里希斯格式解二维对流方程的初值问题 functionu=pe2FL(a,b,dt,nx,minx,maxx,ny,miny,maxy,M) %近似分裂 formatlong; hx=(maxx-minx)/(nx-1); hy=(maxy-miny)/(ny-1); fori=1: nx+4*M forj=1: (ny+4*M) u0(i,j)=Ini2U(minx+(i-2*M-1)*hx,miny+(j-2*M-1)*hy); end end u1=u0; fork=1: M fori=2*k+1: nx+4*M-2*k forj=2*k-1: ny+4*M-2*k+2 tmpU(i,j)=u0(i,j)-a*dt*(u0(i+1,j)-u0(i-1,j))/2/hx+... (a*dt/hx)^2*(u0(i+1,j)-2*u0(i,j)+u0(i-1,j))/2; end end fori=2*k+1: nx+4*M-2*k forj=2*k+1: nx+4*M-2*k u1(i,j)=tmpU(i,j)-b*dt*(tmpU(i,j+1)-tmpU(i,j-1))/2/hy+... (b*dt/hy)^2*(tmpU(i,j+1)-2*tmpU(i,j)+tmpU(i,j-1))/2; end end u0=u1; end u=u1((2*M+1): (2*M+nx),(2*M+1): (2*M+ny)); formatshort; 12.peParabExp用显式格式解扩散方程的初值问题 functionu=peParabExp(c,dt,n,minx,maxx,M) formatlong; h=(maxx-minx)/(n-1); forj=1: (n+2*M) u0(j)=PrIniU(minx+(j-M-1)*h); end u1=u0; fork=1: M fori=k+1: n+2*M-k u1(i)=dt*c*(u0(i+1)-2*u0(i)+u0(i-1))/h/h+u0(i); end u0=u1; end u=u1((M+1): (M+n)); formatshort; 13.peParabTD用跳点格式解扩散方程的初值问题 functionu=peParabTD(c,dt,n,minx,maxx,M) %跳点 formatlong; h=(maxx-minx)/(n-1); forj=1: (n+2*M) u0(j)=PrIniU(minx+(j-M-1)*h); end u1=u0; fork=1: M fori=k+1: n+2*M-k ifmod(n+i,2)==0 u1(i)=u0(i)+c*dt*(u0(i+1)-2*u0(i)+u0(i-1))/h/h; ifi>2 u1(i-1)=(u0(i-1)+c*dt*(u1(i)+u1(i-2))/h/h)/(1+2*c*dt/h/h); end end end u0=u1; end u=u1((M+1): (M+n)); formatshort; 14.peParabImp用隐式格式解扩散方程的初边值问题 functionu=peParabImp(c,dt,n,minx,maxx,lbu,rbu,M) formatlong; h=(maxx-minx)/(n-1); u0 (1)=lbu; u0(n)=rbu; forj=2: n-1 u0(j)=PrIniU(minx+(j-1)*h); end u1=u0; fork=1: M A=zeros(n-2,n-2); cb=-transpose(u0(2: (n-1))); cb (1)=cb (1)-dt*c*lbu/h/h; cb(n-2)=cb(n-2)-dt*c*rbu/h/h; A(1,1)=-2*dt*c/h/h-1; A(1,2)=dt*c/h/h; fori=2: n-3 A(i,i-1)=dt*c/h/h; A(i,i)=-2*dt*c/h/h-1; A(i,i+1)=dt*c/h/h; end A(n-2,n-2)=-2*dt*c/h/h-1; A(n-2,n-3)=dt*c/h/h; u1(2: (n-1))=A\cb; u0=u1; end u=u1; formatshort; 15.peParabKN用克拉克-尼科尔森格式解扩散方程的初边值问题 functionu=peParabKN(c,dt,n,minx,maxx,lbu,rbu,M) formatlong; h=(maxx-minx)/(n-1); u0 (1)=lbu; u0(n)=rbu; forj=2: n-1 u0(j)=PrIniU(minx+(j-1)*h); end u1=u0; fork=1: M A=zeros(n-2,n-2); cb=zeros(n-2,1); cb (1)=-u0 (2)-(u0(3)-2*u0 (2)+u0 (1))*dt*c*lbu/h/h/2-dt*c*lbu/h/h/2; cb(n-2)=-u0(n-1)-(u0(n)-2*u0(n-1)+u0(n-2))*dt*c*lbu/h/h/2-dt*c*rbu/h/h/2; fori=2: n-3 cb(i)=-u0(i+1)-(u0(i+2)-2*u0(i+1)+u0(i))*dt*c*lbu/h/h/2; end A(1,1)=-dt*c/h/h-1; A(1,2)=dt*c/h/h/2; fori=2: n-3 A(i,i-1)=dt*c/h/h/2; A(i,i)=-dt*c/h/h-1; A(i,i+1)=dt*c/h/h/2; end A(n-2,n-2)=-dt*c/h/h-1; A(n-2,n-3)=dt*c/h/h/2; u1(2: (n-1))=A\cb; u0=u1; end u=u1; formatshort; 16.peParabWegImp用加权隐式格式解扩散方程的初边值问题 functionu=peParabWegImp(c,sita,dt,n,minx,maxx,lbu,rbu,M) formatlong; h=(maxx-minx)/(n-1); u0 (1)=lbu; u0(n)=rbu; forj=2: n-1 u0(j)=PrIniU(minx+(j-1)*h); end u1=u0; fork=1: M A=zeros(n-2,n-2); cb=zeros(n-2,1); cb (1)=-u0 (2)-(1-sita)*(u0(3)-2*u0 (2)+u0 (1))*dt*c*lbu/h/h/2... -sita*dt*c*lbu/h/h; cb(n-2)=-u0(n-1)-(1-sita)*(u0(n)-2*u0(n-1)+u0(n-2))*dt*c*lbu/h/h/2... -sita*dt*c*rbu/h/h; fori=2: n-3 cb(i)=-u0(i+1)-(1-sita)*(u0(i+2)-2*u0(i+1)+u0(i))*dt*c*lbu/h/h; end A(1,1)=-2*sita*dt*c/h/h-1; A(1,2)=sit
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分方程 数值 解法
![提示](https://static.bingdoc.com/images/bang_tan.gif)