1、實用文檔標準文案一維拋物方程的初邊值問題求解下歹I問題:分別用向前差分格式、向后差分格式、六點對稱格式,x 1,2 u u 7 a, t xu(x,0) sin x,u(0,t) u(1,t)0,在t 0.05,0.1和0.2時刻的數值解，并與解析解u(x,t)2e sin( x)進仃比較。1差分格式形式設空間步長h 1/ N,時間步長0, T M ，網比r /h2.(1)向前差分格式該問題是第二類初邊值問題(混合問題)，我們要求出所需次數的偏微商的函數 2u(x,t),滿足方程 a2,0 x 1,和初始條件 u(x,0) sin x, 0x1,t x及邊值條件 u(0,t) u(1,t) 0

2、, t 0。已知sin x在相應區域光滑，并且在 x 0,l與邊值相容，使問題有唯一充分光滑的 解。取空間步長h 1/N ,和時間步長T/M ,其中N,M都是正整數。用兩族平行直線 x xj jh( j 0,1,L , N)和 t tk k (k 0,1,L ,M )將矩形域G 0 x 1,t 0分割成矩形網絡，網絡格節點為(xj,tk)。以Gh表示網格內點集合，即位于矩形G的網點集合；Gh表示閉矩形G的網格集合；h Gh Gh是網格界點的集 合。向前差分格式，即k 1 kkk kujujuj 1 2uj uj 1 f a2 fi(1)hfi f(x),其中，j 1,2, ,N1,k1,2,k

3、uj0uj,Mkruj 1k kj(Xj),u0uN 01.以 r a(1 2r)uk/ h2表示網比。kruj 1 fj1)式可改寫成如下：此格式為超格式。其矩陣表達式如下:2r2ru1ju2u1ju22r1 2ruN 1uNuN uN(2)向后差分格式向后差分格式，即k 1ujkujk 1uj 1-k 12ujh2k 1uj 1fj,(2)0ujk k(Xj),u0uN0,其中 j 1,2, N 1,k 1,2,M1.式可改寫成k 1ruj 1(1k 12r)ujk 1ruj 1kujfj此種差分格式被稱為隱格式。 其矩陣表達式如下：1 2r2ru1ju22r2ruN 1uN(3)六點對稱

4、格式六點差分格式:k 1 kuj ujk 1uj 1k2ujh2kujkuj 1k2ujh2kuj_1fj(3)0k kuj j(Xj),u0uN0.將（3）式改寫成31 （1其矩陣表達式如下：k 1 r)Ujr k 12uj 1(1 r)uk 2u：1fj1 r r/2r/2 1 rU1jU2j1 r r/2r/2 1U1ju2r/2r/21 2rUn 1 j 1Unr/2r/2r/21 2rUn 1uN2利用MATLA脈解問題的過程對每種差分格式依次取N 40., =1/1600,=1/3200,=1/6400,用MATLAB求解并圖形比較數值解與精確解，用表格列出不同剖分時的L2誤差。向

5、前差分格式：t 0.05:=1/1600 :1/3200：070.605匚0.30.20.11/6400：070.605D40.30.20 100 0J 0 20.30,40506D.7n.30.91t 0.1:1/3600：050.5-1-1.5*數俏解精確解IItIIIII0.1020.30.40.50.B0.70,811911/3200:a 1_ _o505口Doo2口解解 值確 數精 Dnu1/6400:D.40.360.30.25D20.150.10 050t 0.21/1600：4 n13。X 103i*數值解超確包口.-1 0J 0 20.30.4050607n.30.911/3

6、200：0.140.120.10.080.060 040.020D 0J 0 20.30.4050607n.30.911/6400：0.140J20 10.080.060.040.020向后差分格式:t 0.05:1/1600:1/3200070.B0.50.4030.20 10070.B0.5040.30.20 101/6400t 0.1:1/1600:Q jiIIIIIIIIj0 0J 0 20.30.40506 D7 n.30.911/3200D70.60.5D.4D30.20 100.4解解 值確 數精1/64000.36 -03 -0 25 -02 D45.a 19O8O6 解解 值

7、確 數精0.35 -0.3 -0 25 0.2 t 0.21/1600:0 0J 0 20.30,40506D.7n.30.90 OJ 0 20.3D,40506D.7n.30.910.140,12 0.10.080.060 040.02 01/3200OH0.12 0 10.080.060.040.02 01/64000.140,120.10.080.060 040.020D 0J 0 20.30.4050607n.30.9六點差分格式：t 0.05:1/1600:070.60.504*數值解稀確包030.2口4EIIIIIIIIj0 0J 0 20.30.40506 D7 n.30.911

8、/3200070.60.50.4D30.2數倍解精確解0.10 20.30.40506 D7 C1.80.91/6400070.60.50.4數值解精確解0.30.20.10.102030 40.60.60.70,80 9t 0.1:1/1600:0.36 -03 -0 25 -02 0.4解解 值確 數精1/3200口40.35 -0,3 -0.25 0.2 0.15 -0J -0.05 -D，0 0J 0 20.30.4解解值確數精5.61/6400II 數值解精確解0.40.36030 250 20.150 10.050D 0J 0 20.30.40506t 0.21/1600:0.14

9、0J20.10.080.060.040.0201/32000.140,120.10.080.060 040.02D 0J 0 20.30.4050607n.30.9101/64000.140J20 10.080.060.040.020L2誤差:t=1/1600 t=1/3200t=1/6400向前差分格式t=0.058.888396579076e+210.00140.00034640856587426t=0.18.83785296707723e+590.00170.000422936658240724t=0.21.16572934157246e+1360.001261490899424590.

10、000315223617970853向后差分格式t=0.050.004830905543619830.002765170691253280.00172971295195573t=0.10.005903735046069590.003377972228420410.00211264169361075t=0.20.004408530271260620.002520544967200810.00157579425393478六點差分格式t=0.050.05260005558733240.05251690331377400.0524758656498487t=0.10.0443229427820120

11、0.04425557966346000.0442226561119014t=0.20.02653759444290490.02649801781713880.02647889456215013方法總結及分析本文向前差分格式，向后差分格式及六點差分格式都是使用三對角系數矩陣，計算簡單。根據matlab作，特別明顯的是，t 0.05:1/1600:時，圖像解析解波動特別大，與真1解差距很大。這是因為 r 1 一差分格式不穩定。根據誤差比較，解本題時向后差分格式2誤差最小。衡量一個差分格式是否經濟實用，有點多方面因素決定，應考慮不同的情況決定。附件程序向前差分格式：function u , err

12、= xq( t , t1 )% t是時刻值；% t1是時間步長；N = 40;h = 1 / N;x=h : h : (1-h)'r = t1 / (hA2);u(:,1)=sin(pi * x);%t=0 時刻R = zeros(N-1 , N-1);for i = 2 : N-2R(i , i) = 1 - 2 * r;R(i , i+1) = r;R(i, i-1) = r；endR(1 , 1) = 1 -2 * r;R(1 , 2) = r;R(N-1 , N-1) = 1 - 2 * r;R(N-1, N-2) = r;k = t / t1;for i = 1 : ku(:

13、 , i+1) =R * u(: , i);endu=0 ; u(： , i+1) ; 0;for i=1 : kfor j=1 : N+1up(j , i)=exp(-pi*pi*t) * sin(pi*(j-1)*h);endendx=0 : h:1;plot(x , u ,'b.' , x , up(: , k) ,'r-');legend('數值解','精確解');err=norm(u - up(: , k);end向后差分格式:function u , err = xh( t , t1 )N = 40;h = 1 / N

14、;x=h:h:(1-h)'r = t1 / (hA2);u(:,1)=sin(pi * x);%t=0 時刻R = zeros(N-1 , N-1);for i = 2 : N-2R(i ,i) = 1 + 2 * r;R(i , i+1) = -r;R(i, i-1) = -r；endR(1 , 1) = 1 + 2 * r;R(1 , 2) = -r;R(N-1 , N-1) = 1 + 2 * r;R(N-1, N-2) = -r;k = t / t1;for i = 1 : ku(:,i+1) =inv(R) * u(:,i);endu=0;u(:,i+1);0;for i=1

15、:kfor j=1:N+1upQ,i)=exp(-pi*pi*t)*sin(pi*(j-1)*h);endendx=0:h:1;plot(x,u, 'b.' ,x,up(:,k),'r-');legend('數值解','精確解');err=norm(u-up(:,k);end六點差分格式：function u , err = ld( t , t1 )N = 40;h = 1 / N;x=h:h:(1-h)'r = t1 / (hA2);u(:,1)=sin(pi * x);%t=0 時刻R = zeros(N-1 , N-1);for i = 2 : N-2R(i ,i) = 1 + r;R(i , i+1) = -r/2;R(i, i-1) = -r/2;endR(1 , 1) = 1 + r;R(1 , 2) = -r/2;R(N-1 , N-1) = 1 + 2 * r;R(N-1, N-2) = -r;for i = 2 : N-2R1(i ,i) = 1 - r;R1(i , i+1) = r/2;R1(i, i-1) = r/2;endR1(1 ,1) = 1 - r;R1(1 , 2) = r/2;R1(N-1 , N-1) = 1 -2 * r;R1(N-1, N-2) =