云南財(cái)經(jīng)大學(xué)實(shí)驗(yàn)報(bào)告偏微分方程數(shù)值解法

云南財(cái)經(jīng)大學(xué)實(shí)驗(yàn)報(bào)告系（院：統(tǒng)數(shù)學(xué)專 業(yè)：班 級(jí)：學(xué) 號(hào)：姓 名：實(shí)驗(yàn)時(shí)間：2010612日指導(dǎo)教師：云南財(cái)經(jīng)大學(xué)教務(wù)處制填表說明實(shí)驗(yàn)名稱要用最簡(jiǎn)練的語言反映實(shí)驗(yàn)的內(nèi)容。理論上，驗(yàn)證定理、公式、算法，并使實(shí)驗(yàn)者獲得深刻和系統(tǒng)的理解，在實(shí)踐上，掌握使用實(shí)驗(yàn)設(shè)備的技能技巧和程序調(diào)試的方法。一般需要說明是驗(yàn)證型實(shí)驗(yàn)還是設(shè)計(jì)型實(shí)驗(yàn)，是創(chuàng)新型實(shí)驗(yàn)還是綜合型實(shí)驗(yàn)。實(shí)驗(yàn)環(huán)境實(shí)驗(yàn)用的軟硬件環(huán)境（配置。實(shí)驗(yàn)內(nèi)容（算法、程序、步驟和方法）這是實(shí)驗(yàn)報(bào)告極其重要的內(nèi)容。這部分要寫明依據(jù)何種原理、定律算法、或操作方法進(jìn)行實(shí)驗(yàn)，要寫明經(jīng)過哪幾（實(shí)驗(yàn)裝置的結(jié)構(gòu)示意圖明，這樣既可以節(jié)省許多文字說明，又能使實(shí)驗(yàn)報(bào)告簡(jiǎn)明扼要，清楚明白。結(jié)論（結(jié)果）即根據(jù)實(shí)驗(yàn)過程中所見到現(xiàn)象和測(cè)得的數(shù)據(jù)，作出結(jié)論。小結(jié)對(duì)本次實(shí)驗(yàn)的思考和建議。備注或說明可填寫實(shí)驗(yàn)成功或失敗的原因，實(shí)驗(yàn)后的心得體會(huì)等。寫客觀進(jìn)行評(píng)分、簽名，并記入成績(jī)。u

2u,0x1,t0t x2考慮擴(kuò)散方程u(x,0)sinx,0x1不同差分格式的數(shù)值結(jié)果以及邊界處理對(duì)結(jié)果u(0,t)t)0,t0影響.實(shí)驗(yàn)

u

2u,0xt0名稱 t x2u(x,0)1,0x1考察擴(kuò)散方程的初邊值問題u 差分格式的數(shù)值結(jié)果以及不同邊(xu)|

x0

0,t0界處理對(duì)結(jié)果影響.

(

u)

x1

0,t0實(shí)驗(yàn)考察一些差分格式的數(shù)值結(jié)果以及邊界處理對(duì)結(jié)果的影響.目的

u

2u,0x1,t0t x2考慮擴(kuò)散方程u(x,0)sinx,0x1差分格式的數(shù)值結(jié)果以及邊界處理對(duì)結(jié)果影響.u(0,t)t)0,t0實(shí)驗(yàn)程序:內(nèi)容（算附錄1.1%向前差分的Matlab程序法、functionu=forward(xs,xt,ys,yt,M,N)程h=(xt-xs)/M;k=(yt-ys)/N;%h為空間步長(zhǎng)，k為時(shí)間步長(zhǎng)序、m=M-1;n=N;步驟sigma=k/(h*h);和方a=diag(1-2*sigma*ones(m,1))+diag(sigma*ones(m-1,1),1);法）a=a+diag(sigma*ones(m-1,1),-1);%definematrixalside=l(ys+(0:n)*k);rside=r(ys+(0:n)*k);u(:,1)=f(xs+(1:m)*h)';forj=1:nu(:,j+1)=a*u(:,j)+sigma*[lside(j);zeros(m-2,1);rside(j)];endu=[lside;u;rside];x=(0:m+1)*h;t=(0:n)*k;1subplot(2,2,1);mesh(x,t,u');view(60,30),axis([xsxtysyt-12])%用前向差分計(jì)算問題(2.11)解u(0.4,0.4)的近似值%近似值所在數(shù)據(jù)矩陣的行、列數(shù)sigmai=0.4/0.1j=0.4./k+1附錄1..2functionu=sin(pi*x);附錄1.3functionu=0*t;附錄1.4functionu=0*t;附錄2.1%向后差分的Matlab程序functionu=backward(xs,xt,ys,yt,M,N)h=(xt-xs)/M;k=(yt-ys)/N;%h為空間步長(zhǎng)，k為時(shí)間步長(zhǎng)m=M-1;n=N;sigma=k/(h*h);a=diag(1+2*sigma*ones(m,1))+diag(-sigma*ones(m-1,1),1);a=a+diag(-sigma*ones(m-1,1),-1);%definematrixalside=l(ys+(0:n)*k);rside=r(ys+(0:n)*k);u(:,1)=f(xs+(1:m)*h)';forj=1:nu(:,j+1)=a\(u(:,j)+sigma*[lside(j+1);zeros(m-2,1);rside(j+1)]);endu=[lside;u;rside];x=(0:m+1)*h;t=(0:n)*k;subplot(2,2,2);mesh(x,t,u');view(60,30),axis([xsxtysyt-12])%用后向差分計(jì)算問題(2.11)解u(0.4,0.4)的近似值%近似值所在數(shù)據(jù)矩陣的行、列數(shù)sigmai=0.4/0.1j=0.4./k+1附錄2.2functionu=f(x)u=sin(pi*x);附錄2.32functionu=l(t)u=0*t;附錄2.4functionu=r(t)u=0*t;附錄3.1%Crank-Nilson的Matlab程序functionu=crankn(xs,xt,ys,yt,M,N)h=(xt-xs)/M;k=(yt-ys)/N;%h為空間步長(zhǎng)，k為時(shí)間步長(zhǎng)m=M-1;n=N;sigma=k/(h*h);a=diag(2+2*sigma*ones(m,1))+diag(-sigma*ones(m-1,1),1);a=a+diag(-sigma*ones(m-1,1),-1);%definematrixab=diag(2-2*sigma*ones(m,1))+diag(sigma*ones(m-1,1),1);b=b+diag(sigma*ones(m-1,1),-1);%definematrixblside=l(ys+(0:n)*k);rside=r(ys+(0:n)*k);u(:,1)=f(xs+(1:m)*h)';forj=1:nsides=[lside(j)+lside(j+1);zeros(m-2,1);rside(j)+rside(j)+rside(j+1)];u(:,j+1)=a\(b*u(:,j)+sigma*sides);endu=[lside;u;rside];x=(0:m+1)*h;t=(0:n)*k;subplot(2,2,3);mesh(x,t,u');view(60,30),axis([xsxtysyt-12])%用前向差分計(jì)算問題(2.11)解u(0.4,0.4)的近似值%近似值所在數(shù)據(jù)矩陣的行、列數(shù)sigmai=0.4/0.1j=0.4./k+1附錄3.2functionu=f(x)u=sin(pi*x);附錄3.3functionu=0*t;附錄3.4functionu=0*t;4lamda=[0.25 0.01805439693317 0.01987050634483 0.01895180866708 0.0183519];uf=[0.5 0.01716771002536 0.02079884300382 0.01894069715752 0.0183518];ub=[1.0 inf 0.02269348185685 0.01889625638989 0.0183519];3fprintf('fprintf('(2.11)u(0.4,0.4)的近似值\n')fprintf(' \n')fprintf(' lamdaufubucue\n')fprintf(' \n')fprintf(' \n')步驟:前向差分法來求偏微分方程的數(shù)值解un1unun un ,axbjjj11適當(dāng)?shù)奶岢鲈谟邢迏^(qū)間上的逼近對(duì)流方程的形式u(x,0)f(x)sin(x),axbj1u(a,t)l(t)t0u,t)rt),t0(1)注:類似于常微分方程的情形,偏微分方程有無窮多個(gè)解,.為了適當(dāng)?shù)奶岢鲆粋€(gè)偏微分方程,可能會(huì)用到初始條件和邊界條件的各種結(jié)合.對(duì)于本例,我們已經(jīng)把初始條件f(t和邊界條件l(t)、r(t)結(jié)合起來,從而確定偏微分方程的唯一解.2確定前向差分法來求偏微分方程的解注:分析逼近對(duì)流方程的差分格式可知試圖用第j層的值來計(jì)算第j+1層,用網(wǎng)格形式可以表示為uc=[2.0inf0.01871857872152 0.018718578721520.0183518];ue=[4.0inf0.03493210989403 0.018009304702550.0183519];ul=[8.0a=[lamdasprintf('inf0.05276342929779 0.01520014223247uf ub uc ue];%9.2f %9.7f %9.7f0.0183519];%9.7f %9.7f\n',a')fprintf('0.250.01805440.01987050.01895180.0183519\n')fprintf('0.500.01716770.02079880.01894070.0183518\n')fprintf('1.00--0.02269350.01889630.0183519\n')fprintf('2.00--0.01871860.01871860.0183518\n')fprintf('4.00--0.03493210.01800930.0183519\n')fprintf('8.00--0.05276340.01520010.0183519\n')4為了在時(shí)間[0,1]上離散逼近對(duì)流方程,我們考慮上圖的點(diǎn)的網(wǎng)格,其中實(shí)心圓表示由初始條件和邊界條件而已知的u(x,t)的解,空心圓是網(wǎng)格點(diǎn)并將通過這種方法填成實(shí)心.其離散形式可以按時(shí)間方向向前逐步求解.使用矩陣,可以通過矩陣乘法uj1

Au s來得到在時(shí)間t 處的值j j j1注: s 表示由問題強(qiáng)加的邊界條.j后向差分法來求偏微分方程的數(shù)值解作為供比較的選擇方案,通過使用隱式方法,可以把有限差分方法重新做成具有較好誤差的放大的性質(zhì).用:

1(u(x,t)u(x,tk))ku(x,x), (x,在點(diǎn)

)處將差分公式帶入到方程t k 2 0 j j(1),有:wi1,j

)wiji1,

wi,j

.(2)其中kh2

,并可表示m*m維矩陣方程:5后向差分方程(2)可以看成矩陣迭代w A1w b,其中A=j j定理:設(shè)h是時(shí)間步長(zhǎng),c>0(1),h,k,后向差分穩(wěn)定Crank-Nicolson差分法來求偏微分方程的數(shù)值解1Crank-Nicolson是顯示方法和隱式方法的結(jié)合,它是無條件穩(wěn)定的,并只有誤差階O(h2k2.在方程(1)1 1 1 2u中利用混合差分:

(

2w

)

2w

))代替 ,并用后向差分

1

h2 w

i1,)

ij i1,

2 i1,

i,j1

i1,

x2k再令k

ij,

i,j1

[w

2w w

2w

],或者h(yuǎn)2

i,j

i1,j

ij i1,

i1,j

i,j

i1,j1(2)w (2)w ,這就導(dǎo)出了如下所建i1,j ij i1,j i1,ji,ji1,j的模型:令w (w ,wj 1j 2A=

w )',Crank-NicolsonAwmj

j

(s

j

s其中jB=6s (wj 1

,0 0,w )'mjforward程序步驟matlabforward.mM文件matlabf.mM3matlabl.mM文件4matlabr.mM文件5最后在matlab命令窗口中輸入forward(xs,xt,ys,yt,M,N)backward程序步驟matlabbackward.mM文件matlabf.mM3matlabl.mM文件4matlabr.mM文件5最后在matlab命令窗口中輸入backward(xs,xt,ys,yt,M,N)crankn程序步驟matlabcrankn.mM文件matlabf.mM3matlabl.mM文件4matlabr.mM文件5最后在matlab命令窗口中輸入crankn(xs,xt,ys,yt,M,N)方法前向差分法注:用有限差分方法來求偏微分方程的近似解,想法是對(duì)自變量建立網(wǎng)格并把偏微分方程離散化,把連續(xù)問題變成有限多個(gè)方程的離散化問題.若這個(gè)偏微分方程是線性的,那么它的離散方程也是線性的.后向差分法利用后向差分公式:u

1(u(xtu(xtkku(xx在點(diǎn)xt

)處將差分公式帶入到方t k 2 0 j j程(1),有:wij就是后向差分方法

)wiji1,

wi,j

.(2)其中kh2

,并可表示m*m維矩陣方程:,這Crank-Nilson差分方法Crank-Nicolson是顯示方法和隱式方法的結(jié)合,它是無條件穩(wěn)定的,并只有誤差階O(h2實(shí)驗(yàn)27

k2).u

2u,0xt0t x2u(x,0)1,0x1考察擴(kuò)散方程的初邊值問題u 差分格式的數(shù)值結(jié)果以及不同邊界(xu)|

x0

0,t0程序:1

(

u)

x1

0,t0functionu=forwardmend1(xs,xt,ys,yt,M,N)h=(xt-xs)/M;k=(yt-ys)/N;m=M-1;n=N;sigma=k/(h*h);u(1:5,1)=f1(xs+(1:M/2)*h)';forj=1:nfori=1:(M/2-1)ifi==(M/2-1)u(M/2,j+1)=2*sigma*(u(M/2-1,j)+u(M/2,j));end

end

elseend

u(i+1,j+1)=sigma*u(i,j)+(1-2*sigma)*u(i+1,j)+sigma*u(i+2,j);u(1,j+1)=(1-2*sigma*(1+h))*u(1,j)+2*sigma*u(2,j);x=(0:M/2-1)*h;t=(0:n)*k;subplot(3,3,1);mesh(x,t,u');view(60,30),axis([xsxtysyt-12])%用前向差分計(jì)算問題(2.12)解u(0.2,0.005)的近似值%近似值所在數(shù)據(jù)矩陣的列數(shù)sigmaj=0.005./k+1%forwardmend1(0,1,0,0.005,10,2)附錄2functionu=forwardmend2(xs,xt,ys,yt,M,N)h=(xt-xs)/M;k=(yt-ys)/N;m=M-1;n=N;sigma=k/(h*h);u(1:5,1)=f1(xs+(1:M/2)*h)';forj=1:nfori=1:(M/2-1)ifi==(M/2-1)8u(M/2,j+1)=2*sigma*(u(M/2-1,j)+u(M/2,j));elseu(i+1,j+1)=sigma*u(i,j)+(1-2*sigma)*u(i+1,j)+sigma*u(i+2,j);u(1,j+1)=(1-2*sigma*(1+h))*u(1,j)+2*sigma*u(2,j);end

end

endx=(0:M/2-1)*h;t=(0:n)*k;subplot(3,3,2);mesh(x,t,u');view(60,30),axis([xsxtysyt-12])%近似值所在數(shù)據(jù)矩陣的列數(shù)sigmaj=0.05./k+1

%forwardmend2(0,1,0,0.05,10,20)附錄3functionu=forwardmend3(xs,xt,ys,yt,M,N)h=(xt-xs)/M;k=(yt-ys)/N;m=M-1;n=N;sigma=k/(h*h);u(1:5,1)=f1(xs+(1:M/2)*h)';forj=1:nfori=1:(M/2-1)ifi==(M/2-1)u(M/2,j+1)=2*sigma*(u(M/2-1,j)+u(M/2,j));elseu(i+1,j+1)=sigma*u(i,j)+(1-2*sigma)*u(i+1,j)+sigma*u(i+2,j);u(1,j+1)=(1-2*sigma*(1+h))*u(1,j)+2*sigma*u(2,j);end

end

endx=(0:M/2-1)*h;t=(0:n)*k;subplot(3,3,3);mesh(x,t,u');view(60,30),axis([xsxtysyt-12])%近似值所在數(shù)據(jù)矩陣的列數(shù)sigmaj=0.1./k+1

%forwardmend3(0,1,0,0.1,10,40)附錄4functionu=forwardmend4(xs,xt,ys,yt,M,N)h=(xt-xs)/M;k=(yt-ys)/N;9m=M-1;n=N;sigma=k/(h*h);u(1:5,1)=f1(xs+(1:M/2)*h)';forj=1:nfori=1:(M/2-1)ifi==(M/2-1)u(M/2,j+1)=2*sigma*(u(M/2-1,j)+u(M/2,j));end

end

elseend

u(i+1,j+1)=sigma*u(i,j)+(1-2*sigma)*u(i+1,j)+sigma*u(i+2,j);u(1,j+1)=(1-2*sigma*(1+h))*u(1,j)+2*sigma*u(2,j);x=(0:M/2-1)*h;t=(0:n)*k;subplot(3,3,4);mesh(x,t,u');view(60,30),axis([xsxtysyt-12])%用前向差分計(jì)算問題(2.12)解u(0.2,0.005)的近似值%近似值所在數(shù)據(jù)矩陣的列數(shù)sigmaj=0.25./k+1%forwardmend3(0,1,0,0.25,10,100)附錄5functionu=forwardmend5(xs,xt,ys,yt,M,N)h=(xt-xs)/M;k=(yt-ys)/N;m=M-1;n=N;sigma=k/(h*h);u(1:5,1)=f1(xs+(1:M/2)*h)';forj=1:nfori=1:(M/2-1)ifi==(M/2-1)u(M/2,j+1)=2*sigma*(u(M/2-1,j)+u(M/2,j));end

end

elseend

u(i+1,j+1)=sigma*u(i,j)+(1-2*sigma)*u(i+1,j)+sigma*u(i+2,j);u(1,j+1)=(1-2*sigma*(1+h))*u(1,j)+2*sigma*u(2,j);x=(0:M/2-1)*h;t=(0:n)*k;subplot(3,3,5);mesh(x,t,u');view(60,30),axis([xsxtysyt-12])%用前向差分計(jì)算問題(2.12)解u(0.2,0.005)的近似值%近似值所在數(shù)據(jù)矩陣的列數(shù)10sigmaj=0.5./k+1%forwardmend5(0,1,0,0.5,10,200)附錄6functionu=forwardmend6(xs,xt,ys,yt,M,N)h=(xt-xs)/M;k=(yt-ys)/N;m=M-1;n=N;sigma=k/(h*h);u(1:5,1)=f1(xs+(1:M/2)*h)';forj=1:nfori=1:(M/2-1)ifi==(M/2-1)u(M/2,j+1)=2*sigma*(u(M/2-1,j)+u(M/2,j));end

end

elseend

u(i+1,j+1)=sigma*u(i,j)+(1-2*sigma)*u(i+1,j)+sigma*u(i+2,j);u(1,j+1)=(1-2*sigma*(1+h))*u(1,j)+2*sigma*u(2,j);x=(0:M/2-1)*h;t=(0:n)*k;subplot(3,3,6);mesh(x,t,u');view(60,30),axis([xsxtysyt-12])%用前向差分計(jì)算問題(2.12)解u(0.2,0.005)的近似值%近似值所在數(shù)據(jù)矩陣的列數(shù)sigmaj=1./k+1%forwardmend6(0,1,0,1,10,400)附錄7t=[0.005 1.00000000000000 0.9984 0.0016];sf1=[0.0500.912545015841260.91205.450158412599881e-004];ue=[0.1000.835216321768380.83420.00101632176838];error=[0.250a=[t sprintf('0.64451734712260.6454ue error];%9.3f %9.7f-8.826528773999698e-004];%9.7f %9.7f\n',a')fprintf(' (1)(2.12)的近似值u(0.2,t)\n')fprintf(' \n')fprintf('tsf1ueerror\n')fprintf(' fprintf('0.0051.00000000.99840000.0016000\n')fprintf('0.0500.91254500.91200000.0005450\n')fprintf('0.1000.83521630.83420000.0010163\n')fprintf('0.2500.64451730.6454000-0.0008827\n')fprintf('0.5000.4211401110.4212000-0.0000599\n')fprintf('1.0000.17894480.1794000-0.0004552\n')fprintf(' 步驟:邊界條件用用一階方法處理在內(nèi)點(diǎn)用向前差分格式u u (u 2u u ),利用公式(2.9)i,jij i1,j ij iju 2a2au 2ah(t和公式(2.10)o，j0j 1j nu 2ah)]u 2au 2ah(t)，得到二階精度的邊界處理M.jMj mj n 1u ( u 0,j2 1j 0j 1

由于初邊值問題(2.12)是關(guān)于x

對(duì)稱的因此我們僅考慮1 (u 0.9u )1M，j1 2 mj Mj

i1 5 1 i0

) 2x [0, 就可以了,基于上述考慮,形成計(jì)算格式2

i,j

ui1,j

u uij i1,ju0j1

h))u0j

2u1j利用對(duì)稱性，u4j

u ,所以6j

M,j2

(umjM2

uM2

)

M,j2

2(umj

u )Mj2上述算法稱為算法

u(x,t)4 secll

e2

cos(x1),0x1 初邊值問題(2.12)的解析解

l1

32 l l 2l

,其中

l為方程tan程序步驟

12的正根.1matlabforwardmend1.mM2matlabforwardmend2.mM3matlabforwardmend3.mM4matlabforwardmend4.mM5matlabforwardmend5.mM6matlabforwardmend6.mM7matlab1-table.mM文件方法前向差分法注:用有限差分方法來求偏微分方程的近似解,想法是對(duì)自變量建立網(wǎng)格并把偏微分方程離散化,把連續(xù)問題變成有限多個(gè)方程的離散化問題.若這個(gè)偏微分方程是線性的,那么它的離散方程也是線性的.12實(shí)驗(yàn)1對(duì)于不同的，采用不同差分格式計(jì)算問題(2.11)解u(0.4,0.4)的近似值h0.0025時(shí),此時(shí)其h20.25其數(shù)值解為前向差分法：forward(0,1,0,2.5,10,1000) sigma=0.250 0000.006477383466170.006318870610100.006164236839730.006013387227070.0050.012320715508080.012019206138730.011725075228830.011438142214630.0110.016958010072590.016543018027760.016138181561130.015743252148000.0150.019935336457790.019447484050250.018971570240900.018507302871430.0180.020961253212850.020448294835320.019947889442790.019459729841840.0180.019935336457790.019447484050250.018971570240900.018507302871430.0180.016958010072590.016543018027760.016138181561130.015743252148000.0150.012320715508080.012019206138730.011725075228830.011438142214630.0110.006477383466170.006318870610100.006164236839730.006013387227070.0050000后向差分法：backward(0,1,0,2.5,10,1000)0sigma=0.25000結(jié)論0.007111887979470.006942005024570.006776180094550.006614316254640.0060.013527614812080.013204478229550.012889060469030.012581177149630.012（結(jié)0.018619164454440.018174405104410.017740269801440.017316504766990.016果）0.021888140552660.021365294579120.020854937921950.020356772246580.0190.023014552949940.022464800159670.021928179413760.021404377024700.0200.021888140552660.021365294579120.020854937921950.020356772246580.0190.018619164454440.018174405104410.017740269801440.017316504766990.0160.013527614812080.013204478229550.012889060469030.012581177149630.0120.007111887979470.006942005024570.006776180094550.006614316254640.0060000Crank-Nilson差分法：crankn(0,1,0,2.5,10,1000)sigma=0.2500000.006791107314030.006626925990430.006466713902750.006310375091010.0060.012917453727740.012605162292420.012300420792450.012003046701140.0110.017779349769380.017349517483890.016930076792910.016520776470020.0160.020900879179580.020395581022840.019902498918110.019421337530990.0180.021976484910700.021445182986900.020926725780330.020420802758020.0190.020900879179580.020395581022840.019902498918110.019421337530990.0180.017779349769380.017349517483890.016930076792910.016520776470020.0160.012917453727740.012605162292420.012300420792450.012003046701140.0110.006791107314030.006626925990430.006466713902750.006310375091010.006000013其圖形為:h0.005時(shí),此時(shí)h2前向差分法：

0.5其數(shù)值解為forward(0,1,0,0.5,10,100) sigmaforward(0,1,0,0.5,10,100) sigma=0.50 0000.006167025676150.005865189955470.005578127126450.005305114152340.0050.011730379910930.011156254252910.010610228304680.010090926768540.0090.016145482829670.015355266653880.014603726410630.013888969165020.0130.018980153396840.018051198568340.017167710025360.016327462489490.0150.019956914307020.018980153396840.018051198568340.017167710025360.0160.018980153396840.018051198568340.017167710025360.016327462489490.0150.016145482829670.015355266653880.014603726410630.013888969165020.0130.011730379910930.011156254252910.010610228304680.010090926768540.0090.006167025676150.005865189955470.005578127126450.005305114152340.0050000向差分法：backward(0,1,0,0.5,10,100) sigma=0.50 0 0 00.007435657788320.007088711550080.006757953750800.006442629041300.0060.014143461585040.013483530623670.012854391903030.012254608663600.0110.019466804818520.018558487774540.017692552613980.016867021807020.0160.022884601563170.021816810837450.020798843003820.019828373336530.018140.024062294060420.022939552448930.021869197726380.020848785531450.024062294060420.022939552448930.021869197726380.020848785531450.01987598550.022884601563170.021816810837450.020798843003820.019828373336530.01890318550.019466804818520.018558487774540.017692552613980.016867021807020.01608001000.014143461585040.013483530623670.012854391903030.012254608663600.01168281110.007435657788320.007088711550080.006757953750800.006442629041300.00614201730000Crank-Nilson差分法：crankn(0,1,0,0.5,10,100)sigma=0.500000.006787225177600.006462969825980.006154205566860.005860192323180.0050.012910069465440.012293299135230.011705994613960.011146748191410.0100.017769186204260.016920274672680.016111919347790.015342182682700.0140.020888931192210.019890975834680.018940697157520.018035817437740.0170.021963922053320.020914609693400.019915427561860.018963980719040.0180.020888931192210.019890975834680.018940697157520.018035817437740.0170.017769186204260.016920274672680.016111919347790.015342182682700.0140.012910069465440.012293299135230.011705994613960.011146748191410.0100.006787225177600.006462969825980.006154205566860.005860192323180.0050000當(dāng)h0.1,0.01時(shí),此時(shí)h2前向差分法：

1.0其數(shù)值解為forward(0,1,0,1,10,100) sigma=1.0 1.0e+02900000000000000-0.000000000000000.00000000000000-0.000000000000000.00000000000-0.000000000000000.00000000000000-0.000000000000000.00000000000000-0.000000000000.00000000000000-0.000000000000000.00000000000000-0.000000000000000.00000000000-0.000000000000000.00000000000000-0.000000000000000.00000000000000-0.000000000000.00000000000000-0.000000000000000.00000000000000-0.000000000000000.00000000000-0.000000000000000.00000000000000-0.000000000000000.00000000000000-0.000000000000.00000000000000-0.000000000000000.00000000000000-0.000000000000000.00000000000-0.000000000000000.00000000000000-0.000000000000000.00000000000000-0.000000000000.00000000000000-0.000000000000000.00000000000000-0.000000000000000.000000000000000后向差分法：backward(0,1,0,1,10,100) sigma=1.00 0000.010712915820500.009757758438270.008887762336120.008095334583570.0070.020377176799220.018560359494290.016905528570110.015398241414590.0140.028046777736690.025546143245390.023268463879900.021193861090090.0190.032970964655910.030031292505180.027353719824210.024914877975780.0220.034667723832380.031576769614250.028761403087560.026197053013040.0230.032970964655910.030031292505180.027353719824210.024914877975780.0220.028046777736690.025546143245390.023268463879900.021193861090090.0190.020377176799220.018560359494290.016905528570110.015398241414590.0140.010712915820500.009757758438270.008887762336120.008095334583570.0070000Crank-Nilson差分法：crankn(0,1,0,1,10,100)sigma=1.000000.009085207834590.008237380037360.007468671175740.006771697903750.0060.017281092225970.015668427923460.014206256779510.012880534835490.0110.023785382905830.021565740916050.019553234988890.017728535273560.0160.027961394584350.025352048930430.022986206322160.020841143157090.0180.029400350142460.026656721757390.024169127626300.021913674739630.0190.027961394584350.025352048930430.022986206322160.020841143157090.0180.023785382905830.021565740916050.019553234988890.017728535273560.0160.017281092225970.015668427923460.014206256779510.012880534835490.0110.009085207834590.008237380037360.007468671175740.006771697903750.006000016當(dāng)h0.1,0.02時(shí),此時(shí)h2前向差分法：

2.0其數(shù)值解為forward(0,1,0,2,10,100)0forward(0,1,0,2,10,100)0sigma=2.001.0e+065*000.00000000000000-0.000000000000000.00000000000000-0.000000000000000.0000.000000000000000.00000000000000-0.000000000000000.00000000000000-0.0000.00000000000000-0.000000000000000.00000000000000-0.000000000000000.0000.000000000000000.00000000000000-0.000000000000000.00000000000000-0.0000.00000000000000-0.000000000000000.00000000000000-0.000000000000000.0000.000000000000000.00000000000000-0.000000000000000.00000000000000-0.0000.000000000000000.000000000000000.00000000000000-0.000000000000000.0000.000000000000000.00000000000000-0.000000000000000.00000000000000-0.0000.000000000000000.000000000000000.00000000000000-0.000000000000000.0000000向差分法：backward(0,1,0,2,10,100) sigma=2.00 0 0 00.012367982961450.010343077902360.008649693392020.007233552379890.0060.023525301577830.019673703275180.016452694528880.013759034253710.0110.032379799765350.027078529496670.022645191292580.018937685989950.0150.038064737548520.031832720583820.026621018954240.022262585074880.0180.040023633607810.033470903188610.027990995801120.023408267220130.0190.038064737548520.031832720583820.026621018954240.022262585074880.018170.032379799765350.027078529496670.022645191292580.018937685989950.032379799765350.027078529496670.022645191292580.018937685989950.01583717910.023525301577830.019673703275180.016452694528880.013759034253710.01150638410.012367982961450.010343077902360.008649693392020.007233552379890.00604926410000Crank-Nilson差分法：crankn(0,1,0,2,10,100)sigma=2.000000.009008287832400.007401940360130.006082034912110.004997493477700.0040.017134781687340.014079327225470.011568717870990.009505797474220.0070.023584003725670.019378531445530.015922974120660.013083607783170.0100.027724659159920.022780829989540.018718578721520.015380703403460.0120.029151431786540.023953182170790.019681878417110.016172228610950.0130.027724659159920.022780829989540.018718578721520.015380703403460.0120.023584003725670.019378531445530.015922974120660.013083607783170.0100.017134781687340.014079327225470.011568717870990.009505797474220.0070.009008287832400.007401940360130.006082034912110.004997493477700.0040000當(dāng)h0.1,0.04時(shí),此時(shí)h2前向差分法：

4.0其數(shù)值解為forward(0,1,0,4,10,100)sigma=4.01.0e+099*00000.000000000000000.000000000000000.000000000000000.000000000000000.000180.000000000000000.000000000000000.000000000000000.000000000000000.000000000000000.000000000000000.000000000000000.000000000000000.00000000000.000000000000000.000000000000000.000000000000000.000000000000000.00000000000.000000000000000.000000000000000.000000000000000.000000000000000.00000000000.000000000000000.000000000000000.000000000000000.000000000000000.00000000000.000000000000000.000000000000000.000000000000000.000000000000000.00000000000.000000000000000.000000000000000.000000000000000.000000000000000.00000000000.000000000000000.000000000000000.000000000000000.000000000000000.00000000000.000000000000000.000000000000000.000000000000000.000000000000000.00000000000000后向差分法：backward(0,1,0,4,10,100) sigma=4.00 0000.042559200355570.030584072085590.021978454893920.015794249966860.0110.080952409652960.058174362103670.041805505489920.030042448701940.0210.111421433064910.080070140234440.057540341932490.041349883240040.0290.130983750289700.094128095157590.067642728799560.048609703105020.0340.137724465418670.098972136297710.071123774077140.051111266546370.0360.130983750289700.094128095157590.067642728799560.048609703105020.0340.111421433064910.080070140234440.057540341932490.041349883240040.0290.080952409652960.058174362103670.041805505489920.030042448701940.0210.042559200355570.030584072085590.021978454893920.015794249966860.0110000Crank-Nilson差分法：crankn(0,1,0,4,10,100)sigma=4.000000.028599387779310.019234716889690.012936442440020.008700494213850.0050.054399268219120.036586605674070.024606575760510.016549323434140.0110.074874169263660.050357142581200.033868046001470.022778189570780.0150.088019864941660.059198371513640.039814275927250.026777367807260.0180.092549562968710.062244851383010.041863207122900.028155390713860.0180.088019864941660.059198371513640.039814275927250.026777367807260.0180.074874169263660.050357142581200.033868046001470.022778189570780.0150.054399268219120.036586605674070.024606575760510.016549323434140.0110.028599387779310.019234716889690.012936442440020.008700494213850.005000019當(dāng)h0.1,0.08時(shí),此時(shí)h2前向差分法：

8.0其數(shù)值解為forward(0,1,0,8,10,100)0forward(0,1,0,8,10,100)0sigma=8.001.0e+131*000.00000000000000