五点差分格式.doc

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

1、如有侵权，请联系网站删除，仅供学习与交流五点差分格式【精品文档】第 8 页微分方程数值解大作业（一） 椭圆型方程编程计算：采用五点差分格式求如下椭圆型方程其中、及边条件为：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=

2、x1+(x2-x1)*(0:M)/M; %x，y向量表示横纵坐标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(

3、i+(j-1)*m,:)=zeros(1,i-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)+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

4、 a(i+(j-1)*m,:)=zeros(1,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=

5、1;a(1,:)=t -1 zeros(1,m-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)

6、t -1 zeros(1,m-2);b(i+(j-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;

7、for i=1:m for j=1:nu1(i+(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.484746714701678 0.837645626697

8、549 1.139019509919315 0.594429507664308 0.915886065952874 1.197402289453010 0.753941698688471 1.034066839996629 1.287784599003526a = 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

9、 0 0 0 0 -1 0 -1 4b = 0.506911724444854 0.810930216216329 2.521030123526701 0.223143551314210 0 1.446918982936325 1.387270447092946 1.178654996341646 2.919669266564466运行ni(1,2,0,1,4,4)u1 =精确解 Columns 1 through 3 0.485507815781701 0.838329190404443 1.139434283188365 Columns 4 through 6 0.594707107746

10、693 0.916290731874155 1.197703191312341 Columns 7 through 9 0.753771802376380 1.034073767530539 1.287854288306638误差很小（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)fi

11、vepoints(0,1,0,2,4,4)u = 0.062500000000000 -0.000000000000000 0.062500000000000 0.562500000000000 0.250000000000000 0.062500000000000 1.562500000000000 1.000000000000000 0.562500000000000a = Columns 1 through 3 2.500000000000000 -1.000000000000000 0 -1.000000000000000 2.500000000000000 -1.0000000000

12、00000 0 -1.000000000000000 2.500000000000000 -0.250000000000000 0 0 0 -0.250000000000000 0 0 0 -0.250000000000000 0 0 0 0 0 0 0 0 0 Columns 4 through 6 -0.250000000000000 0 0 0 -0.250000000000000 0 0 0 -0.250000000000000 2.500000000000000 -1.000000000000000 0 -1.000000000000000 2.500000000000000 -1.

13、000000000000000 0 -1.000000000000000 2.500000000000000 -0.250000000000000 0 0 0 -0.250000000000000 0 0 0 -0.250000000000000 Columns 7 through 9 0 0 0 0 0 0 0 0 0 -0.250000000000000 0 0 0 -0.250000000000000 0 0 0 -0.250000000000000 2.500000000000000 -1.000000000000000 0 -1.000000000000000 2.500000000

14、000000 -1.000000000000000 0 -1.000000000000000 2.500000000000000b = 0.015625000000000 -0.187500000000000 0.140625000000000 0.750000000000000 -0.250000000000000 -0.250000000000000 2.765625000000000 0.312500000000000 0.390625000000000精确解ni(0,1,0,2,4,4)u1 =u1 = Columns 1 through 3 0.062500000000000 0 0

15、.062500000000000 Columns 4 through 6 0.562500000000000 0.250000000000000 0.062500000000000 Columns 7 through 9 1.562500000000000 1.000000000000000 0.562500000000000误差很小（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;funct

16、ion g=left(y)g=cosd(y);function g=f1(x,y)g=cosd(x)*cosd(y);数值解Pi=3.1415926fivepoints(0,pi,0,pi/2,4,4)u = 0.657818362488653 0.000000024999241 -0.657818327134387 0.504980798089256 0.000000019229497 -0.504980770894641 0.273644362624153 0.000000010432161 -0.273644347870850a = 10 -1 0 -4 0 0 0 0 0 -1 10

17、-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.558260407530267 0.000000137720159 -4.558260212764549 1.323957061211835 0.000000023374742 -1.323957028154957 0.716520423452347 0

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