五点差分格式

1、精选优质文档-倾情为你奉上微分方程数值解大作业（一） 椭圆型方程编程计算：采用五点差分格式求如下椭圆型方程其中、及边条件为：1 ， 且边条件如下：问题存在精确解为： 2 ，且边条件如下：问题存在精确解为： 3 ，且边条件如下：问题存在精确解为： .代码：主函数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

2、向量表示横纵坐标y=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); endend%下边缘j=1;for i=2:m-1 a(i+(j-1)*m,:)=zeros(1,i-2

3、) -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)+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,

4、i-1+(j-2)*m) -r 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-

5、2) -r zeros(1,(n-1)*m-1);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

6、-1)*m)=h1*f(x(i+1),y(j+1)+r*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=abab2，精确解：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:nu1(i+

7、(j-1)*m)=f1(x(i+1),y(j+1) endend（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 =数值解 0.1678 0.7549 1.9315 0.4308 0.2874 1.3010 0.8471 1.6629

8、1.3526a = 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 -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 4b = 0.4854 0.6329 2.6701 0.4210 0 1.6325 1.2946 1.1646 2.4466运行ni(1,2,0,1,4,4)u1 =精确解 Columns 1 through 3

9、 0.1701 0.4443 1.8365 Columns 4 through 6 0.6693 0.4155 1.2341 Columns 7 through 9 0.6380 1.0539 1.6638误差很小（2）辅助函数function g=f(x,y)g=-4;function g=bottom(x)g=x2;function g=right(y)g=(y-1)2;function g=top(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

10、,1,0,2,4,4)u = 0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 1.0000 0.0000a = Columns 1 through 3 2.0000 -1.0000 0 -1.0000 2.0000 -1.0000 0 -1.0000 2.0000 -0.0000 0 0 0 -0.0000 0 0 0 -0.0000 0 0 0 0 0 0 0 0 0 Columns 4 through 6 -0.0000 0 0 0 -0.0000 0 0 0 -0.0000 2.0000 -1.0000 0 -1.0000 2.0000

11、 -1.0000 0 -1.0000 2.0000 -0.0000 0 0 0 -0.0000 0 0 0 -0.0000 Columns 7 through 9 0 0 0 0 0 0 0 0 0 -0.0000 0 0 0 -0.0000 0 0 0 -0.0000 2.0000 -1.0000 0 -1.0000 2.0000 -1.0000 0 -1.0000 2.0000b = 0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 2.0000 0.0000 0.0000精确解ni(0,1,0,2,4,4)u1 =u1 = Columns 1 th

12、rough 3 0.0000 0 0.0000 Columns 4 through 6 0.0000 0.0000 0.0000 Columns 7 through 9 1.0000 1.0000 0.0000误差很小（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);

13、数值解Pi=3.fivepoints(0,pi,0,pi/2,4,4)u = 0.8653 0.9241 -0.4387 0.9256 0.9497 -0.4641 0.4153 0.2161 -0.0850a = 10 -1 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 10b = 4.0267 0.0159 -4.4549 1.1835 0.4742 -1.4957 0.2347 0.0320 -0.2093精确解 ni(0,pi,0,p