拋物型方程隱格式的程序?qū)崿F(xiàn)

1、淮 海 工 學(xué) 院實(shí) 驗(yàn) 報(bào) 告 書課程名稱： 微分方程數(shù)值解法 實(shí)驗(yàn)名稱： 拋物型方程隱格式的程序?qū)崿F(xiàn) 班 級(jí) 信息101 姓 名： 學(xué)號(hào)： 日 期： 地點(diǎn) 數(shù)學(xué)實(shí)驗(yàn)室 指導(dǎo)教師： 成績(jī)： 數(shù) 理 科 學(xué) 系1. 實(shí)驗(yàn)?zāi)康模喊凑諏?shí)驗(yàn)內(nèi)容對(duì)-方法導(dǎo)出的幾種隱式差分格式編程實(shí)現(xiàn)求出數(shù)值解；估計(jì)數(shù)值解的誤差界；大致比較幾種不同格式的優(yōu)劣；結(jié)合格式的相容性、穩(wěn)定性和收斂性條件簡(jiǎn)單分析計(jì)算結(jié)果。2. 實(shí)驗(yàn)內(nèi)容：用-方法導(dǎo)出的隱式格式求如下初邊值問題 （2）的數(shù)值解，觀察取不同值時(shí)，數(shù)值解的誤差情況,大致比較至少三種不同格式的優(yōu)劣，分析解釋計(jì)算結(jié)果。3. 實(shí)驗(yàn)步驟：第一步: 對(duì)求解區(qū)域作網(wǎng)格剖分。我們把

2、求解區(qū)域分割為一些均勻的矩形，其邊與坐標(biāo)軸平行，取空間步長為和時(shí)間步長為，其中，為正整數(shù)，用兩族平行直線：平行于 軸作豎直線和平行于 軸作水平線,將矩形區(qū)域分割成矩形網(wǎng)格，得網(wǎng)格剖分圖。第二步: 差分法的目的是：求初邊值問題（2）的解在節(jié)點(diǎn)的近似值。為此，需要構(gòu)造逼近微分方程定解問題的差分格式。采用-方法，構(gòu)造如下差分格式以表示網(wǎng)比。由于本是實(shí)驗(yàn)中, 所以上述差分格式簡(jiǎn)化為利用Taylor展式，可得差分格式的截?cái)嗾`差為：利用Fourier方法判斷差分格式的穩(wěn)定性條件是：當(dāng)，當(dāng)且僅當(dāng)時(shí)穩(wěn)定，當(dāng)，對(duì)所有的格式均穩(wěn)定。利用最大值原理可以證明差分格式的收斂條件是：。第三步: 有限差分方程的解法。根據(jù)問

3、題的規(guī)模和計(jì)算機(jī)的容量和速度，選取適當(dāng)?shù)慕夥?。這是數(shù)值計(jì)算的關(guān)鍵一步，有限差分法解方程的主要計(jì)算都集中在這里。由于本實(shí)驗(yàn)問題的有限差分方程為-方法導(dǎo)出的隱式差分格式，可以用Thomas方法求解。對(duì)每個(gè)時(shí)間層需要求解的方程為：初邊值條件為。第四步: 根據(jù)前面的分析，編制程序，上機(jī)計(jì)算，求出數(shù)值解，必要時(shí)畫出解的幾何圖形，分析觀察解的性質(zhì)。為了方便比較可以用傅立葉級(jí)數(shù)方法求出定解問題的解析解為：。另外，上述差分程組為三對(duì)角方程組，可用Thomas法（追趕法）求解。若有三對(duì)角方程組：，則Thomas法（追趕法）公式為， 其中為所求三對(duì)角方程組的解。四 實(shí)驗(yàn)數(shù)據(jù)記錄及分析（或程序及運(yùn)行結(jié)果）：程序如下

4、：clearcita=input(cita=);M=input(M=);J=input(J=);u=zeros(J+1,M+1);d=zeros(1,J);e=zeros(1,J-1);f=zeros(1,J);x=zeros(1,J);uu=zeros(1,J);for n=1:M u(n,1)=0; u(n,M)=0;endfor j=1:J u(1,j)=sin(j*pi/J);endtt=1/M;xx=1/J;r=tt/(xx*xx); for n=1:M for j=2:J-1; d(1,j)=(1-cita)*r*u(n,j+1)+(1-2*(1-cita)*r)*u(n,j)+(

5、1-cita)*r*u(n,j-1); a=cita*r; b=1+2*cita*r; c=cita*r; end e(1,2)=c/b; f(1,2)=d(1,2)/b; for i=3:J-2 e(1,i)=c/(b-e(1,i-1)*a); end for i=3:J-1 f(1,i)=(d(1,i)+f(1,i-1)*a)/(b-e(1,i-1)*a); end x(1,J-1)=f(1,J-1); for i=J-2:-1:2 x(1,i)=f(1,i)+e(1,i)*x(1,i+1); end for i=2:J-1 u(n+1,i)=x(1,i); endend for j=1:

6、J-1 uu(1,j)=exp(-pi*pi)*sin(pi*j/J); h(j)=uu(1,j)-u(M,j);endxuhmesh(u) ;運(yùn)行結(jié)果如下：cita=0M=10J=10x = 1.0e+010 * Columns 1 through 9 0 -2.4999 3.6685 -3.2575 2.0465 -0.9383 0.3102 -0.0706 0.0099 Column 10 0u = 1.0e+010 * Columns 1 through 9 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0

7、 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0 0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0 -0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0 0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0 -0

8、.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 0 0.0001 -0.0001 0.0001 -0.0000 0.0000 -0.0000 -0.0000 0.0000 0 -0.0023 0.0031 -0.0024 0.0012 -0.0004 0.0001 -0.0000 -0.0000 0 0.0754 -0.1067 0.0887 -0.0505 0.0200 -0.0053 0.0009 -0.0001 0 -2.4999 3.6685 -3.2575 2.0465 -0.9383 0.3102 -0.0706

9、 0.0099 Columns 10 through 11 0.0000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0h = 1.0e+009 * 0.0000 -0.7543 1.0667 -0.8874 0.5047 -0.2001 0.0534 -0.0087 0.0007cita=0.5M=10J=10x = 1.0e-003 * Columns 1 through 9 0 -0.9215 0.2009 0.2907 0.0655 -0.0568 -0.0596 -0.0252 -0.0037 Column 10 0u = Columns 1 th

10、rough 9 0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090 0.5878 0.3090 0 0.1335 0.2338 0.2981 0.3250 0.3148 0.2704 0.1975 0.1041 0 0.0163 0.0425 0.0664 0.0813 0.0846 0.0760 0.0571 0.0305 0 0.0139 0.0174 0.0183 0.0185 0.0178 0.0156 0.0117 0.0063 0 -0.0029 0.0012 0.0046 0.0061 0.0059 0.0049 0.0034 0.0

11、017 0 0.0037 0.0017 0.0005 0.0004 0.0008 0.0010 0.0009 0.0006 0 -0.0023 -0.0002 0.0009 0.0009 0.0005 0.0002 0.0000 -0.0000 0 0.0019 0.0001 -0.0005 -0.0003 0.0000 0.0002 0.0002 0.0001 0 -0.0014 0.0001 0.0005 0.0002 0.0000 -0.0001 -0.0001 -0.0000 0 0.0011 -0.0002 -0.0004 -0.0001 0.0001 0.0001 0.0001 0

12、.0000 0 -0.0009 0.0002 0.0003 0.0001 -0.0001 -0.0001 -0.0000 -0.0000 Columns 10 through 11 0.0000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0h = 0.0000 -0.0011 0.0002 0.0004 0.0002 -0.0000 -0.0000 -0.0000 -0.0000cita=1M=10J=10x = 1.0e-003 * Columns 1 through 9 0 0.1294 0.2432 0.3277 0.3726 0.3726 0.32

13、77 0.2432 0.1294 Column 10 0u = Columns 1 through 9 0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090 0.5878 0.3090 0 0.1834 0.3264 0.4211 0.4628 0.4508 0.3887 0.2847 0.1503 0 0.0757 0.1406 0.1869 0.2098 0.2075 0.1808 0.1333 0.0706 0 0.0334 0.0625 0.0839 0.0950 0.0945 0.0828 0.0613 0.0325 0 0.0150 0.

14、0282 0.0379 0.0430 0.0429 0.0377 0.0279 0.0149 0 0.0068 0.0127 0.0171 0.0195 0.0195 0.0171 0.0127 0.0068 0 0.0031 0.0058 0.0078 0.0088 0.0088 0.0078 0.0058 0.0031 0 0.0014 0.0026 0.0035 0.0040 0.0040 0.0035 0.0026 0.0014 0 0.0006 0.0012 0.0016 0.0018 0.0018 0.0016 0.0012 0.0006 0 0.0003 0.0005 0.0007 0.0008 0.0008 0.0007 0.0005 0.0003 0 0.0001 0.0002 0.0003 0.0004 0.0004 0.0003 0.0002 0.0001 Columns 10 through 11 0.0000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0h = 1.0e-003 * 0.0160 -0.2551 -0.4947 -0.6737 -0.7703 -0.7729 -0.6810 -0.5061 -0.26956、思考題1.什么是隱式格式? 寫出隱格式的截?cái)?/p>