微分方程數(shù)值解第二次上機(jī)報(bào)告

1、南京信息工程大學(xué) 實(shí)驗(yàn)（實(shí)習(xí)）報(bào)告實(shí)驗(yàn)課程 微分方程數(shù)值解 實(shí)驗(yàn)名稱 第二次實(shí)驗(yàn) 實(shí)驗(yàn)日期 2016 指導(dǎo)老師 專業(yè) 年級(jí) 姓名 學(xué)號(hào) 得分 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 實(shí)驗(yàn)?zāi)康模?對(duì)一維拋物方程（熱傳導(dǎo)方程）數(shù)值求解的向前差分格式、向后差分格式、CrankNicolson格式、交替顯隱式格式以及跳點(diǎn)格式進(jìn)行編程。*求解一道拋物方程定解問題實(shí)例,編程給出解三維圖

2、以及末時(shí)刻精確解與數(shù)值解比較圖。*以理論推導(dǎo)和編程實(shí)現(xiàn)分析差分方法分別對(duì)空間步長及時(shí)間步長的整體誤差精度（收斂階）（以向前差分為例）。實(shí)驗(yàn)內(nèi)容：1. 方法編程已編程向前差分格式、向后差分格式、CrankNicolson格式、交替顯隱式格式以及跳點(diǎn)格式。分別見附錄的forward.m文件、backward.m文件、Crank_Nicolson.m文件、AFB.m文件和skip.m文件。以向前差分為例探討差分方法分別對(duì)空間步長及時(shí)間步長的收斂階的程序分別見附錄中errforwardx.m文件和errforwardt.m文件。2. 實(shí)例求解下求解一個(gè)一維熱傳導(dǎo)方程定解問題實(shí)例：用偏微分方程相關(guān)知識(shí)可

3、以求得其解析解為：解析解參與編程時(shí)，其中k取60,可以達(dá)到精確解精度要求。分別運(yùn)行五個(gè)差分格式，得出解三維圖和末時(shí)刻精確解與格式數(shù)值解對(duì)比圖。由于五個(gè)格式運(yùn)行結(jié)果無法從肉眼觀察出區(qū)別，下取向前差分格式運(yùn)行結(jié)果展示。如下圖圖1、圖2和圖3。 圖1：向前差分求解結(jié)果三維圖圖2：末時(shí)刻數(shù)值解與精確解比較圖3：數(shù)值解與精確解比較放大圖結(jié)果分析：從圖1可以看出關(guān)于時(shí)間逐漸降低，在邊界上滿足零邊界值條件，符合熱傳導(dǎo)物理規(guī)律。從圖2可以看出向前差分格式在空間步長時(shí)間步長時(shí)，有很好的收斂性。末時(shí)刻差分方法解與精確解較為貼合。經(jīng)過多次放大后（如圖3）才可以看出細(xì)微差別。這是因?yàn)榇藭r(shí)網(wǎng)格比,滿足穩(wěn)定性條件。而五種

4、方法原理不同，取不同空間步長、時(shí)間步長時(shí)的收斂性也不同。下面研究差分方法的收斂階，以深入研究收斂性。3. 收斂階理論推導(dǎo) 將方程在節(jié)點(diǎn)離散化： （1）時(shí)間步長為，空間步長為.*下以向前差分差分格式為例：對(duì)充分光滑的解，由Taylor展式：（2） （3） （4）對(duì)（2）式移項(xiàng)得： （5）（3）、（4）式相加得：（6）（5）、（6）式帶入（1）式得： （7）其中： （8）（7）式舍去即得到逼近（1）式的向前格式差分方程： （9）其中，從（8）式來看，網(wǎng)格化近似后余項(xiàng)對(duì)時(shí)間步長的局部誤差階為2，對(duì)空間步長的局部誤差階為3.于是對(duì)時(shí)間、空間兩種步長的整體誤差階為1和2.4. 收斂階編程探討以向前差分格

5、式為例，先固定時(shí)間步長，變動(dòng)空間步長：固定時(shí)間步長，分別取四個(gè)空間步長，這樣取值是為了避免在右邊界處數(shù)值計(jì)算時(shí)有時(shí)為，有時(shí)取不到，以影響末時(shí)刻結(jié)果精確性。計(jì)算末時(shí)刻相對(duì)誤差，誤差與步長分別取對(duì)數(shù)后繪圖如圖4。圖4：末時(shí)刻向前差分方法相對(duì)誤差隨空間步長變化對(duì)數(shù)-對(duì)視圖圖中其實(shí)有三條直線，程序中用矩陣U存放了三條直線的斜率.分別約為：1.403、2.135、2.056.符合理論討論中的關(guān)于空間步長達(dá)到2階收斂精度。再固定空間步長，變動(dòng)時(shí)間步長：固定空間步長，取兩個(gè)空間步長，再計(jì)算末時(shí)刻相對(duì)誤差。誤差與步長分別取對(duì)數(shù)后繪圖如圖5。同時(shí)計(jì)算了這條直線的斜率約為符合理論討論中的關(guān)于時(shí)間步長達(dá)到1階收斂精

6、度。圖5：末時(shí)刻向前差分方法相對(duì)誤差隨時(shí)間步長變化對(duì)數(shù)-對(duì)視圖5. 附錄所有Matlab程序如下：forward.m文件：clcclearl=1; %參數(shù)l的值a=1; %系數(shù)a的值tmax=0.1; %t最大值k=0.0002; %時(shí)間t增量h=0.02; % x增量s=a(2)*k/(h2); %網(wǎng)格比x=0:h:l;t=0:k:tmax;o=length(x);p=length(t);u=zeros(o,p);X,T=meshgrid(x,t);% u(x,0) 初始層for i=1:o if x(i)=1/2 u(i,1)=2*x(i); else u(i,1)=2-2*x(i); e

7、ndend% u(0,t)和u(l,t) 邊界條件u(1,:)=0;u(o,:)=0;% 向前差分for j=1:(p-1) for i=2:(o-1) u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j); endend% 末時(shí)刻精確解u_exact=zeros(60,o);for i=1:o for k=1:60 %取60項(xiàng)u_exact(k,i)=8*sin(k*pi/2)*exp(-0.1*(k*pi)2)*sin(k*pi*x(i)/(pi2);u_e(i)=sum(u_exact(:,i); endend% 解三維網(wǎng)格繪圖figure()mes

8、h(X,T,u)xlabel(x);ylabel(t);zlabel(u(x,t);title(向前差分求解三維圖)% 末時(shí)刻解比較figure()plot(x,u(:,p),-) %差分方法求解hold onplot(x,u_e,-.) %精確解xlabel(x);ylabel(u(x,t);title(末時(shí)刻向前差分方法解與精確解比較圖)legend(向前差分方法解,精確解（n=60）)backward.m文件：clcclearl=1; %參數(shù)l的值a=1; %系數(shù)a的值tmax=0.1; %t最大值k=0.00004; %時(shí)間t增量h=0.01; % x增量x=0:h:l;t=0:k:t

9、max;o=length(x);p=length(t);u=zeros(o,p);s=a(2)*k/(h2); %網(wǎng)格比N=o-2; %每一層內(nèi)點(diǎn)個(gè)數(shù)m,n=meshgrid(x,t);% u(x,0) 初始層for i=1:o if x(i)=1/2 u(i,1)=2*x(i); else u(i,1)=2-2*x(i); endend% u(0,t)和u(l,t) 邊界條件u(1,:)=0;u(o,:)=0;E=zeros(N,1); %邊界不為零在此改動(dòng)F=zeros(N,1); %隱式方程組右端 先聲明%隱式差分方程組系數(shù)矩陣AA=zeros(N,N); A(1,1)=1+2*s;A(

10、1,2)=-s;A(N,N)=1+2*s;A(N,N-1)=-s;for i=2:N-1 A(i,i)=1+2*s; A(i,i+1)=-s; A(i,i-1)=-s;end% 對(duì)時(shí)間層逐層求解for j=1:p-1 F=u(2:N+1,j)+E; %方程組右端更新 u(2:N+1,j+1)=AF;end% 末時(shí)刻精確解u_exact=zeros(60,o);for i=1:o for k=1:60 %取60項(xiàng)u_exact(k,i)=8*sin(k*pi/2)*exp(-0.1*(k*pi)2)*sin(k*pi*x(i)/(pi2);u_e(i)=sum(u_exact(:,i); end

11、end% 解三維網(wǎng)格繪圖figure()mesh(m,n,u)xlabel(x);ylabel(t);zlabel(u(x,t);title(隱式差分求解三維圖)% 末時(shí)刻解比較figure()plot(x,u(:,p),-) %隱式差分方法求解hold onplot(x,u_e,-.) %精確解xlabel(x);ylabel(u(x,t);title(末時(shí)刻隱式差分方法解與精確解比較圖)legend(隱式差分方法解,精確解（n=60）)Crank_Nicolson.m文件：clcclearl=1; %參數(shù)l的值a=1; %系數(shù)a的值tmax=0.1; %t最大值k=0.0002; %時(shí)間t

12、增量h=0.02; % x增量x=0:h:l;t=0:k:tmax;o=length(x);p=length(t);u=zeros(o,p);s=a(2)*k/(2*h2); %C-N網(wǎng)格比N=o-2; %每一層內(nèi)點(diǎn)個(gè)數(shù)m,n=meshgrid(x,t);% u(x,0) 初始層for i=1:o if x(i)=1/2 u(i,1)=2*x(i); else u(i,1)=2-2*x(i); endend% u(0,t)和u(l,t) 邊界條件u(1,:)=0;u(o,:)=0;E=zeros(N,1); %非齊次邊界問題 在此改動(dòng)F=zeros(N,1); %隱式方程組右端 先聲明%方程組

13、 隱式一端系數(shù)矩陣A A=zeros(N,N); A(1,1)=1+2*s;A(1,2)=-s;A(N,N)=1+2*s;A(N,N-1)=-s;for i=2:N-1 A(i,i)=1+2*s; A(i,i+1)=-s; A(i,i-1)=-s;end%方程組 顯式一端系數(shù)矩陣T T=zeros(N,N); T(1,1)=1-2*s;T(1,2)=s;T(N,N)=1-2*s;T(N,N-1)=s;for i=2:N-1 T(i,i)=1-2*s; T(i,i+1)=s; T(i,i-1)=s;end% 對(duì)時(shí)間層逐層求解for j=1:p-1 F=u(2:N+1,j)+E; %方程組右端更新

14、 u(2:N+1,j+1)=inv(A)*T*F;end% 末時(shí)刻精確解u_exact=zeros(60,o);for i=1:o for k=1:60 %取60項(xiàng)u_exact(k,i)=8*sin(k*pi/2)*exp(-0.1*(k*pi)2)*sin(k*pi*x(i)/(pi2);u_e(i)=sum(u_exact(:,i); endend% 解三維網(wǎng)格繪圖figure()mesh(m,n,u)xlabel(x);ylabel(t);zlabel(u(x,t);title(C-N差分求解三維圖)% 末時(shí)刻解比較figure()plot(x,u(:,p),-) %隱式差分方法求解h

15、old onplot(x,u_e,-.) %精確解xlabel(x);ylabel(u(x,t);title(末時(shí)刻C-N差分方法解與精確解比較圖)legend(C-N差分方法解,精確解（n=60）)AFB.m文件：clcclear%一維交替顯隱式格式l=1; %參數(shù)l的值a=1; %系數(shù)a的值tmax=0.1; %t最大值k=0.0002; %時(shí)間t增量h=0.02; % x增量x=0:h:l;t=0:k:tmax;o=length(x);p=length(t);u=zeros(o,2*p-1); %共2p-1層 前p層為初始層+隱式矯正層 后p-1層為顯式預(yù)測(cè)層s=a(2)*k/(2*h2

16、); %交替顯隱式格式網(wǎng)格比N=o-2; %每一層內(nèi)點(diǎn)個(gè)數(shù)m,n=meshgrid(x,t);% u(x,0) 初始層for i=1:o if x(i)=1/2 u(i,1)=2*x(i); else u(i,1)=2-2*x(i); endend% u(0,t)和u(l,t) 邊界條件u(1,:)=0;u(o,:)=0;E=zeros(N,1); %非齊次邊界問題 在此改動(dòng)F=zeros(N,1); %隱式方程組右端 先聲明%方程組 隱式一端系數(shù)矩陣A A=zeros(N,N); A(1,1)=1+2*s;A(1,2)=-s;A(N,N)=1+2*s;A(N,N-1)=-s;for i=2:

17、N-1 A(i,i)=1+2*s; A(i,i+1)=-s; A(i,i-1)=-s;end% 時(shí)間層循環(huán)for j=1:p-1 for i=2:(o-1) u(i,j+p)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j); %顯式預(yù)估 end F=u(2:N+1,j+p)+E; %方程組右端更新 u(2:N+1,j+1)=AF; %隱式矯正 end% 末時(shí)刻精確解u_exact=zeros(60,o);for i=1:o for k=1:60 %取60項(xiàng)u_exact(k,i)=8*sin(k*pi/2)*exp(-0.1*(k*pi)2)*sin(k*pi*x(i

18、)/(pi2);u_e(i)=sum(u_exact(:,i); endend% 解三維網(wǎng)格繪圖figure()mesh(m,n,u(:,1:p)xlabel(x);ylabel(t);zlabel(u(x,t);title(交替顯隱式格式求解三維圖)% 末時(shí)刻解比較figure()plot(x,u(:,p),-) %隱式差分方法求解hold onplot(x,u_e,-.) %精確解xlabel(x);ylabel(u(x,t);title(末時(shí)刻交替顯隱式格式解與精確解比較圖)legend(交替顯隱式格式解,精確解（n=60）)skip.m文件：clcclear%一維跳點(diǎn)格式l=1; %參

19、數(shù)l的值a=1; %系數(shù)a的值tmax=0.1; %t最大值k=0.0002; %時(shí)間t增量h=0.02; % x增量x=0:h:l;t=0:k:tmax;o=length(x);p=length(t);u=zeros(o,p);s=a(2)*k/(h2); %跳點(diǎn)格式網(wǎng)格比N=o-2; %每一層內(nèi)點(diǎn)個(gè)數(shù)m,n=meshgrid(x,t);% u(x,0) 初始層for i=1:o if x(i)=1/2 u(i,1)=2*x(i); else u(i,1)=2-2*x(i); endend% u(0,t)和u(l,t) 邊界條件u(1,:)=0;u(o,:)=0;%跳點(diǎn)格式循環(huán)for j=1

20、:p-1 for i=2:(o-1) if rem(i+j+1,2)=0 u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j); %(i+j+1)先在偶數(shù)格點(diǎn)顯式 end if rem(i+j+1,2)=0 u(i,j+1)=u(i+1,j)/(1+2*s)+s*u(i+1,j+1)+s*u(i-1,j+1); %(i+j+1)再在奇數(shù)格點(diǎn)隱式（實(shí)際上顯式） end end end% 末時(shí)刻精確解u_exact=zeros(60,o);for i=1:o for k=1:60 %取60項(xiàng)u_exact(k,i)=8*sin(k*pi/2)*exp(-0.1

21、*(k*pi)2)*sin(k*pi*x(i)/(pi2);u_e(i)=sum(u_exact(:,i); endend% 解三維網(wǎng)格繪圖figure()mesh(m,n,u(:,1:p)xlabel(x);ylabel(t);zlabel(u(x,t);title(跳點(diǎn)格式求解三維圖)% 末時(shí)刻解比較figure()plot(x,u(:,p),-) %跳點(diǎn)格式求解hold onplot(x,u_e,-.) %精確解xlabel(x);ylabel(u(x,t);title(末時(shí)刻跳點(diǎn)格式解與精確解比較圖)legend(跳點(diǎn)格式解,精確解（n=60）)errforwardx.m文件:clcc

22、learformat long%向前差分格式收斂階l=1; %參數(shù)l的值a=1; %系數(shù)a的值tmax=0.1; %t最大值k=0.00004; %時(shí)間t增量h0=0.0125; %初始空間步長err=zeros(1,4); %存放相對(duì)誤差ke=1;m=log( 1 2 4 8*h0); for hh= 1 2 4 8 %対空間步長循環(huán) h=hh*h0; % x增量s=a(2)*k/(h2); %網(wǎng)格比x=0:h:l;t=0:k:tmax;o=length(x);p=length(t);u=zeros(o,p);% u(x,0) 初始層for i=1:o if x(i)=1/2 u(i,1)=

23、2*x(i); else u(i,1)=2-2*x(i); endend% u(0,t)和u(l,t) 邊界條件u(1,:)=0;u(o,:)=0;% 向前差分for j=1:(p-1) for i=2:(o-1) u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j); endend% 末時(shí)刻精確解u_exact=zeros(60,o);u_e=zeros(1,o);for ii=1:o for kk=1:60 %取60項(xiàng)u_exact(kk,ii)=8*sin(kk*pi/2)*exp(-0.1*(kk*pi)2)*sin(kk*pi*x(ii)/(pi

24、2);u_e(ii)=sum(u_exact(:,ii); endenderr(ke)=norm(u(:,p)-u_e,2)/norm(u_e,2); %末時(shí)刻相對(duì)誤差ke=ke+1; end; %各段斜率計(jì)算len=length(err);U=zeros(1,len-1);for flag=1:len-1 U(flag)=(log(err(flag+1)-log(err(flag)/(m(flag+1)-m(flag);end % 末時(shí)刻相對(duì)誤差plot(m,log(err),-) xlabel(log（h）);ylabel(log（error）);title(末時(shí)刻向前差分方法相對(duì)誤差隨空間步長變化對(duì)數(shù)-對(duì)數(shù)圖)legend(向前差分方法)errforwardt.m文件:clcclearformat longl=1; %參數(shù)l的值a=1; %系數(shù)a的值tmax=0.1; %t最大值h=0.01; %x增量k0=0.00