偏微分方程數(shù)值解上機實習(xí)

1、偏微分方程數(shù)值解上機實習(xí)數(shù)值求解二維擴散方程的初邊值問題 古典顯式格式：將原格式化為：附源程序：%-運用古典顯式差分格式求解二維擴散方程的初邊值問題;function gdxs(ti,h,t)%-ti:時間步長;%-h:空間步長;k=t/ti;m=1/h+1;r=ti/h2; %- r為網(wǎng)格比;w=ones(m,m);u=ones(m,m);for i=2:m-1 for j=2:m-1 u(i,j)=sin(pi*(i-1)*h)*sin(2*pi*h*(j-1); endendticfor l=1:k for i=2:m-1 for j=2:m-1 w(i,j)=r*u(i-1,j)+r*

2、u(i,j-1)+r*u(i+1,j)+r*u(i,j+1)+(1-4*r)*u(i,j); end end u=w;endtoct=tocumesh(u)交替方向隱式格式（P-R格式）：將原差分格式化為：代入邊界條件，轉(zhuǎn)化為三對角矩陣附追趕法源程序：%-追趕法求解三對角方程組;function x=zg(a,b,c,d)%-a:方程組系數(shù)矩陣A的下對角元素;%-b:方程組系數(shù)矩陣A的主對角元素;%-c:方程組系數(shù)矩陣A的上對角元素;%-d:追趕法所求方程的右端向量;%-l:系數(shù)矩陣A所分解成的下三角陣L中的下對角元素了l(i);%-u:系數(shù)矩陣A所分解成的下三角陣U中的主對角元素了u(i);

3、n=length(b);u(1)=b(1);y(1)=d(1);for i=1:n-1 %-追趕法求解之追過程 求解Ly=d; l(i)=a(i)/u(i); u(i+1)=b(i+1)-l(i)*c(i); y(i+1)=d(i+1)-l(i)*y(i);endx(n)=y(n)/u(n); %-追趕法求解之趕過程 求解Uz=y;for j=n-1:-1:1 if u(j)=0 break; else x(j)=(y(j)-c(j)*x(j+1)/u(j); endend%-運用P-R差分格式求解二維擴散方程的初邊值問題;function pr(ti,h,t) %-ti:時間步長h:空間步長

4、;k=t/ti+1;m=1/h+1;r=ti/h2; %- r為網(wǎng)格比;w=ones(m,m);u=ones(m,m); %-輸入初始值v=ones(m,m);for i=2:m-1 for j=2:m-1 u(i,j)=sin(pi*(i-1)*h)*sin(2*pi*h*(j-1); endend%-輸入用P-R差分格式求解的三對角矩陣b=ones(1,m-2)*(2+2*r);a=-r*ones(1,m-3);c=-r*ones(1,m-3);A=zeros(m-2,m-2);for i=1:m-2 A(i,i)=2-2*r; endfor i=1:m-3 A(i,i+1)=r; A(i

5、+1,i)=r;endp=zeros(m-2,1);p(1)=2*r;p(m-2)=2*r;ticfor l=1:k for i=2:m-1 d1=A*u(i,2:m-1)'+p; d1=d1' w(2:m-1,i)=zg(a,b,c,d1); %-調(diào)用追趕法求解 d2=A*w(2:m-1,i)+p; v(i,2:m-1)=zg(a,b,c,d2); %-調(diào)用追趕法求解 end u=v'endtoc t=tocumesh(0:0.1:1,0:0.1:1,u)局部一維格式： 將原格式化為：代入邊界條件，轉(zhuǎn)化為三對角矩陣附源程序：%-運用局部一維格式求解二維擴散方程的初邊值

6、問題;function god(ti,hi,t) %-ti為時間步長 ， hi為空間步長;m=1/hi;n=t/ti;g=ti/(hi2); %- g為網(wǎng)格比u=ones(m+1,m+1); %-輸入初始值for i=2:m for j=2:m u(i,j)=sin(pi*(i-1)*hi)*sin(2*pi*(j-1)*hi); endenda(1:m-2)=-0.5*g;b(1:m-1)=1+g;c(1:m-2)=-0.5*g; %-輸入用局部一維差分格式求解的三對角矩陣B=zeros(m-1,m+1);for i=1:m-1 B(i,i)=0.5*g; B(i,i+1)=1-g; B(i

7、,i+2)=0.5*g;endf=zeros(m-1,1);f(1,1)=0.5*g;f(m-1,1)=0.5*g;w=ones(m+1,m+1);for i=1:n for j=2:m d=B*u(:,j)+f;%-調(diào)用追趕法求解 x=zg(a,b,c,d); w(2:m,j)=x' end for j=2:m e=B*w(j,:)'+f; x=zg(a,b,c,e); %-調(diào)用追趕法求解 u(j,2:m)=x; end endumesh(u)古典顯式在t=1時運行結(jié)果：gdxs（0.0025，0.1，1）所用時間t=0 1.000000000000001.000000000

8、000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000000.999999999707990.999999999445260.999999999235500.999999999102410.999999999055020.999999999102400.999999999235500.999999999445260.99999

9、9999707991.000000000000001.000000000000000.999999999445260.999999998943480.999999998547660.999999998290520.999999998204810.999999998290520.999999998547660.999999998943480.999999999445261.000000000000001.000000000000000.999999999235500.999999998547660.999999997998510.999999997650070.999999997526020.9

10、99999997650070.999999997998510.999999998547660.999999999235501.000000000000001.000000000000000.999999999102400.999999998290520.999999997650070.999999997234010.999999997095320.999999997234010.999999997650070.999999998290520.999999999102401.000000000000001.000000000000000.999999999055020.9999999982048

11、10.999999997526020.999999997095320.999999996941990.999999997095320.999999997526020.999999998204810.999999999055021.000000000000001.000000000000000.999999999102400.999999998290520.999999997650070.999999997234010.999999997095320.999999997234010.999999997650070.999999998290520.999999999102401.000000000

12、000001.000000000000000.999999999235500.999999998547660.999999997998510.999999997650070.999999997526020.999999997650070.999999997998510.999999998547660.999999999235501.000000000000001.000000000000000.999999999445260.999999998943480.999999998547660.999999998290520.999999998204810.999999998290520.99999

13、9998547660.999999998943480.999999999445261.000000000000001.000000000000000.999999999707990.999999999445260.999999999235500.999999999102400.999999999055020.999999999102400.999999999235500.999999999445260.999999999707991.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.0

14、00000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.00000000000000P-R 格式t=1時運行結(jié)果：pr（0.0025，0.1，1）所用時間t=0.360000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.00000000

15、0000001.000000000000001.000000000000001.000000000000000.999999999544420.999999999133430.999999998807270.999999998597870.999999998525710.999999998597870.999999998807270.999999999133430.999999999544421.000000000000001.000000000000000.999999999133430.999999998351690.999999997731300.999999997332980.9999

16、99997195730.999999997332980.999999997731300.999999998351690.999999999133431.000000000000001.000000000000000.999999998807270.999999997731300.999999996877400.999999996329170.999999996140260.999999996329170.999999996877400.999999997731300.999999998807271.000000000000001.000000000000000.999999998597870.

17、999999997332980.999999996329170.999999995684680.999999995462600.999999995684680.999999996329170.999999997332980.999999998597871.000000000000001.000000000000000.999999998525710.999999997195730.999999996140260.999999995462600.999999995229100.999999995462600.999999996140260.999999997195730.999999998525

18、711.000000000000001.000000000000000.999999998597870.999999997332980.999999996329170.999999995684680.999999995462600.999999995684680.999999996329170.999999997332980.999999998597871.000000000000001.000000000000000.999999998807270.999999997731300.999999996877400.999999996329170.999999996140260.99999999

19、6329170.999999996877400.999999997731300.999999998807271.000000000000001.000000000000000.999999999133430.999999998351690.999999997731300.999999997332980.999999997195730.999999997332980.999999997731300.999999998351690.999999999133431.000000000000001.000000000000000.999999999544420.999999999133430.9999

20、99998807270.999999998597870.999999998525710.999999998597870.999999998807270.999999999133430.999999999544421.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.00000000000000局部

21、一維格式t=1時的運行結(jié)果：god（0.0025，0.1，1）所用時間t= 0.390000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000000.999999999521570.999999999089960.999999998747440.9999

22、99998527530.999999998451750.999999998527530.999999998747440.999999999089960.999999999521571.000000000000001.000000000000000.999999999089960.999999998269010.999999997617490.999999997199200.999999997055060.999999997199200.999999997617490.999999998269010.999999999089961.000000000000001.000000000000000.

23、999999998747440.999999997617490.999999996720760.999999996145030.999999995946640.999999996145030.999999996720760.999999997617500.999999998747441.000000000000001.000000000000000.999999998527530.999999997199200.999999996145030.999999995468210.999999995234990.999999995468210.999999996145030.999999997199

24、200.999999998527531.000000000000001.000000000000000.999999998451750.999999997055060.999999995946640.999999995234990.999999994989770.999999995234990.999999995946640.999999997055060.999999998451751.000000000000001.000000000000000.999999998527530.999999997199200.999999996145030.999999995468210.999999995234990.999999995468210.999999996145030.999999997199200.9999999985