五点差分格式.pdf

《五点差分格式.pdf》由会员分享，可在线阅读，更多相关《五点差分格式.pdf（9页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、 微分方程数值解大作业（一）椭圆型方程 编程计算：采用五点差分格式求如下椭圆型方程 2222uuxyf(x,y),(x,y);其中f(x,y)、及边条件为：1 f(x,y)0,(1,2)(0,1)，且边条件如下：222u(x,0)2lnx,u(x,1)ln(x1)1x2;u(1,y)ln(1y),u(2,y)ln(4y),0y1.，问题存在精确解为：22(,)ln()u x yxy 2f(x,y)4,(0,1)(0,2)，且边条件如下：2222u(x,0)x,u(x,2)(x2)0 x1;u(0,y)y,u(1,y)(y1),0y2.，:问题存在精确解为：2(,)()u x yxy 3f(x,

2、y)cos(xy)cos(xy),(0,)(0,)2 ，且边条件如下：u(x,0)cosx,u(x,)0 0 x;2u(0,y)cosy,u(,y)cosy,0y.2 ，问题存在精确解为：(,)coscosu x yxy.代码：主函数 1，差分解 function g=fivepoints(x1,x2,y1,y2,M,N)%变步长法 h=(x2-x1)/M;%横轴步长 k=(y2-y21/N;%纵轴步长 m=M-1;n=N-1;h1=h2;r=h1/k2;%五点中的上下两个点的系数 t=2+2*r;%五点中的中心点的系数 x=x1+(x2-x1)*(0:M)/M;%x，y 向量表示横纵坐标（y

3、=y1+(y2-y1)*(0:N)/N;a=zeros(m*n,m*n);b=zeros(m*n,1);%初始化 a，b 矩阵，a 为系数矩阵%内部的（m-2）*（n-2）个点 for i=2:m-1 for j=2:n-1 a(i+(j-1)*m,:)=zeros(1,i-1+(j-2)*m)-r zeros(1,m-2)-1 t-1 zeros(1,m-2)-r zeros(1,(n-j)*m-i);b(i+(j-1)*m)=h1*f(x(i+1),y(j+1);end end%下边缘！j=1;for i=2:m-1 a(i+(j-1)*m,:)=zeros(1,i-2)-1 t-1 ze

4、ros(1,m-2)-r zeros(1,(n-j)*m-i);b(i+(j-1)*m)=h1*f(x(i+1),y(j+1)+r*bottom(x(i+1);end;%右边缘 i=m;for j=2:n-1 a(i+(j-1)*m,:)=zeros(1,(j-1)*m-1)-r zeros(1,m-2)-1 t zeros(1,m-1)-r zeros(1,(n-j)*m-i);b(i+(j-1)*m)=h1*f(x(i+1),y(j+1)+right(y(j+1);end：%上边缘 j=n;for i=2:m-1 a(i+(j-1)*m,:)=zeros(1,i-1+(j-2)*m)-r

5、zeros(1,m-2)-1 t-1 zeros(1,m-i-1);b(i+(j-1)*m)=h1*f(x(i+1),y(j+1)+r*top(x(i+1);end%左边缘 i=1;for j=2:n-1 a(i+(j-1)*m,:)=zeros(1,i-1+(j-2)*m)-r zeros(1,m-1)t-1 zeros(1,m-2)-r zeros(1,(n-j)*m-i);b(i+(j-1)*m)=h1*f(x(i+1),y(j+1)+left(y(j+1);end;%左下角的那个点 i=1;j=1;a(1,:)=t-1 zeros(1,m-2)-r zeros(1,(n-1)*m-1)

6、;b(1)=h1*f(x(2),y(2)+r*bottom(x(2)+left(y(2);%右下角的那个点 i=m;j=1;a(i+(j-1)*m,:)=zeros(1,m-2)-1 t zeros(1,m-1)-r zeros(1,(n-2)*m);b(i+(j-1)*m)=h1*f(x(i+1),y(j+1)+r*bottom(x(i+1)+right(y(j+1);%左上角的那个点 i=1;j=n;a(i+(j-1)*m,:)=zeros(1,(n-2)*m)-r zeros(1,m-1)t-1 zeros(1,m-2);b(i+(j-1)*m)=h1*f(x(i+1),y(j+1)+r

7、*top(x(i+1)+left(y(j+1);%右上角的那个点 i=m;j=n;a(i+(j-1)*m,:)=zeros(1,(n-1)*m-1)-r zeros(1,m-2)-1 t;b(i+(j-1)*m)=h1*f(x(i+1),y(j+1)+r*top(x(i+1)+right(y(j+1);u=ab a b 2，精确解：function g=ni(x1,x2,y1,y2,M,N)m=M-1;￥n=N-1;x=x1+(x2-x1)*(0:M)/M;y=y1+(y2-y1)*(0:N)/N;for i=1:m for j=1:n u1(i+(j-1)*m)=f1(x(i+1),y(j+

8、1)end end （1）辅助函数function g=f(x,y)g=0;function g=bottom(x)g=2*log(x);function g=right(y)g=log(4+y2);function g=top(x)g=log(x2+1);function g=left(y)g=log(1+y2);function g=f1(x,y)g=log(x2+y2);运行fivepoints(1,2,0,1,4,4)【u=数值解 :a=4 -1 0 -1 0 0 0 0 0 -1 4 -1 0 -1 0 0 0 0 0 -1 4 0 0 -1 0 0 0 -1 0 0 4 -1 0

9、-1 0 0 0 -1 0 -1 4 -1 0 -1 0 0 0 -1 0 -1 4 0 0 -1 0 0 0 -1 0 0 4 -1 0 0 0 0 0 -1 0 -1 4 -1 0 0 0 0 0 -1 0 -1 4 b=0）运行ni(1,2,0,1,4,4)u1=精确解 Columns 1 through 3 、Columns 4 through 6 Columns 7 through 9 误差很小（2）辅助函数function g=f(x,y)g=-4;function g=bottom(x)！g=x2;function g=right(y)g=(y-1)2;function g=to

10、p(x)g=(x-2)2;function g=left(y)g=y2;function g=f1(x,y)g=(x-y)2;fivepoints(1,2,0,1,4,4)!fivepoints(0,1,0,2,4,4)u=￥a=Columns 1 through 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Columns 4 through 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Columns 7 through 9 0 0 0：0 0 0 0 0 0 0 0 0 0 0 0 0 0 b=精确解 ni(0,1,0,2,4,4)u1=u1=Co

11、lumns 1 through 3 0 Columns 4 through 6 Columns 7 through 9 误差很小 （3）辅助函数function g=f(x,y)g=cosd(x+y)+cosd(x-y);function g=bottom(x)g=cosd(x);function g=right(y)g=-cosd(y);function g=top(x)g=0;function g=left(y)g=cosd(y);function g=f1(x,y)g=cosd(x)*cosd(y);数值解 Pi=fivepoints(0,pi,0,pi/2,4,4)u=?a=10 -1

12、 0 -4 0 0 0 0 0 -1 10 -1 0 -4 0 0 0 0 0 -1 10 0 0 -4 0 0 0 -4 0 0 10 -1 0 -4 0 0 0 -4 0 -1 10 -1 0 -4 0 0 0 -4 0 -1 10 0 0 -4 0 0 0 -4 0 0 10 -1 0 0 0 0 0 -4 0 -1 10 -1 0 0 0 0 0 -4 0 -1 10 b=精确解 ni(0,pi,0,pi/2,4,4)u1=Columns 1 through 3 Columns 4 through 6 Columns 7 through 9 误差很小 注：（1）需要对数值解与精确解作比较，以及不同步长选取下的误差比较。（2）请在 2011-10-23 日之前提交(只需提交电子版给)。