五點(diǎn)差分格式

1、精選優(yōu)質(zhì)文檔-傾情為你奉上微分方程數(shù)值解大作業(yè)（一） 橢圓型方程編程計(jì)算：采用五點(diǎn)差分格式求如下橢圓型方程其中、及邊條件為：1 ， 且邊條件如下：問題存在精確解為： 2 ，且邊條件如下：問題存在精確解為： 3 ，且邊條件如下：問題存在精確解為： .代碼：主函數(shù)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; %五點(diǎn)中的上下兩個點(diǎn)的系數(shù)t=2+2*r; %五點(diǎn)中的中心點(diǎn)的系數(shù)x=x1+(x2-x1)*(0:M)/M; %x，y

2、向量表示橫縱坐標(biāo)y=y1+(y2-y1)*(0:N)/N;a=zeros(m*n,m*n); b=zeros(m*n,1);%初始化a，b矩陣，a為系數(shù)矩陣%內(nèi)部的（m-2）*（n-2）個點(diǎn)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;%左下角的那個點(diǎn)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);%右下角的那個點(diǎn)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);%左上角的那個點(diǎn)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);%右上角的那個點(diǎn)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）輔助函數(shù)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);運(yùn)行fivepoints(1,2,0,1,4,4)u =數(shù)值解 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運(yùn)行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誤差很?。?）輔助函數(shù)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誤差很?。?）輔助函數(shù)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、數(shù)值解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