资源描述:
《偏微分方程的数值解法》由会员上传分享,免费在线阅读,更多相关内容在教育资源-天天文库。
1、偏微分方程的数值解法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:nyu0(j,1)=EllIni2Uxl(minx,miny+(j-1)*hy);u0(j,nx)=EllIni2Uxr(maxx,miny+(j-1)*hy);endforj=1:nxu0(1,j)=EllIni2Uyl(minx+(j-1)*hx,miny
2、);u0(ny,j)=EllIni2Uyr(minx+(j-1)*hx,maxy);endA=-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)==1ifi==1A(1,2)=1;A(1,nx-1)=1;b(1)=-u0(1,2)-u0(2,1);elseifi==(ny-3)*(nx-2)+1A(i,i+1)=1;A(i,i-nx+2)=1;b(i)=-u0(ny-1,1)-u0(ny,2);elseA(i,i+1)=1;A(i,
3、i-nx+2)=1;A(i,i+nx-2)=1;b(i)=-u0(floor(i/(nx-2))+2,1);endendelseifmod(i,nx-2)==0ifi==nx-2A(i,i-1)=1;A(i,i+nx-2)=1;b(i)=-u0(1,nx-1)-u0(2,nx);elseifi==(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);elseA(i,i-1)=1;A(i,i-nx+2)=1;A(i,i+nx-2)=1;b(i)=-u0(floor(i/(nx-2))+1,
4、nx);endendelseifi>1&&i(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);elseA(i,i-1)=1;A(i,i+1)=1;A(i,i+nx-2)=1;A(i,i-nx+2)=1;endendendendendul=Ab;fori=1:(ny-2)forj=1:(nx-2)u(i
5、,j)=ul((i-1)*(nx-2)+j);endendformatshort;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:nyu0(j,1)=EllIni2Uxl(minx,miny+(j-1)*hy);u0(j,nx)=EllIni2Uxr(maxx,miny+(j-1)*hy);endforj=1:nxu0(
6、1,j)=EllIni2Uyl(minx+(j-1)*hx,miny);u0(ny,j)=EllIni2Uyr(minx+(j-1)*hx,maxy);endA=-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)==1ifi==1A(1,nx)=1;b(1)=-u0(1,1)-u0(3,1)-u0(1,3);elseifi==(ny-3)*(nx-2)+1A(i,i-nx+1)=1;b(i)=-u0(ny,1)-u0(ny,3)-u
7、0(ny-2,1);elseA(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);endendelseifmod(i,nx-2)==0ifi==nx-2A(i,i+nx-1)=1;b(i)=-u0(1,nx-2)-u0(1,nx)-u0(3,nx);elseifi==(ny-2)*(nx-2)A(i,i-nx+1)=1;b(i)=-u0(ny,nx)-u0(ny,nx-2)-u0(ny-2,nx);elseA(i,i-nx+1)=1;A(i,