一維導(dǎo)熱方程有限差分法matlab實(shí)現(xiàn)

1、第五次作業(yè)(前三題寫在作業(yè)紙上)一、用有限差分方法求解一維非定常熱傳導(dǎo)方程，初始條件和邊界條件見說明.pdf文件,熱擴(kuò)散系數(shù)a =const,:T；：2T=2:t: x1. 用Tylaor展開法推導(dǎo)出FTCS式的差分方程2. 討論該方程的相容性和穩(wěn)定性，并說明穩(wěn)定性要求對(duì)求解差分方程的影響。3. 說明該方程的類型和定解條件，如何在程序中實(shí)現(xiàn)這些定解條件。4. 編寫M文件求解上述方程，并用適當(dāng)?shù)奈淖謱?duì)程序做出說明。(部分由網(wǎng)絡(luò)搜索得至L添加，修改后得到。)function rechuandaopde%以下所用數(shù)據(jù)，除了t的范圍我根據(jù)題目要求取到了20000，其余均從pdf中得來a=0.00001

2、;%a 的取值xspan=0 1;%x的取值范圍tspan=0 20000;%t 的取值范圍ngrid=100 10;%分割的份數(shù)，前面的是t軸的，后面的是x軸的f=(x)0;% 初值g1=(t)100;% 邊界條件一g2=(t)100;% 邊界條件二T,x,t=pdesolution(a,f,g1,g2,xspan,tspan,ngrid);% 計(jì)算所調(diào)用的函數(shù)x,t=meshgrid(x,t);mesh(x,t,T);%畫圖，并且把坐標(biāo)軸名稱改為x, t, Txlabel('x')ylabel('t')zlabel('T')T%輸出溫度矩陣d

3、t=tspan(2)/ngrid(1);%t 步長h3000=3000/dt;h9000=9000/dt;h15000=15000/dt;%3000,9000,15000下，溫度分別在 T矩陣的哪些行T3000=T(h3000,:)T9000=T(h9000,:)T15000=T(h15000,:)%輸出三個(gè)時(shí)間下的溫度分布%不再對(duì)三個(gè)時(shí)間下的溫度-長度曲線畫圖，其圖像就是三維圖的截面%穩(wěn)定性討論，傅里葉級(jí)數(shù)法dx=xspan(2)/ngrid(2);%x 步長sta=4*a*dt/(dxA2)*(sin(pi/ 2)A2;if sta>0,sta<2fprintf('n%

4、sn','有穩(wěn)定性')elsefprintf('n%sn','沒有穩(wěn)定性')errorend%真實(shí)值計(jì)算xe,te,Te=truesolution(a,f,g1,g2,xspan,tspan,ngrid);xe,te=meshgrid(xe,te);mesh(xe,te,Te);%畫圖，并且把坐標(biāo)軸名稱改為xe, te, Texlabel('xe')ylabel('te')zlabel('Te')Te%諭出溫度矩陣%誤差計(jì)算jmax=1/dx+1;% 網(wǎng)格點(diǎn)數(shù)rms=wuchajisuan(

5、T,Te,jmax)rms%輸出誤差function rms=wuchajisuan(T,Te,jmax)for j=1:jmaxrms=(T(j)-Te(j)A2/jmax)A(1 /2)endfunctionUe,xe,te=truesolution(a,f,g1,g2,xspan,tspan,ngrid)n=ngrid(1);%t 份數(shù)m=ngrid(2);%x 份數(shù)Ue=zeros(ngrid);xe=linspace(xspan(1),xspan(2),m);% 畫網(wǎng)格te=linspace(tspan(1),tspan(2),n);% 畫網(wǎng)格for j=2:nfor i=2:m-1

6、for g=1:m-1Ue(j,i)=100-(400/(2*g-1)/pi)*sin(2*g-1)*pi*xe(j)*exp(-a*(2*g-1)A2*piA2*te(i)endendendfunction U,x,t=pdesolution(a,f,g1,g2,xspan,tspan,ngrid)n=ngrid(1);%t 份數(shù)m=ngrid(2);%x 份數(shù)h=range(xspan)/(m-1);%x 網(wǎng)格長度x=linspace(xspan(1),xspan(2),m);% 畫網(wǎng)格k=range(tspan)/(n-1); %t 網(wǎng)格長度t=linspace(tspan(1),tsp

7、an,n);% 畫網(wǎng)格U=zeros(ngrid);U(:,1)=g1(t);% 邊界條件U(:,m)=g2(t);U(1,:)=f(x);% 初值條件%差分計(jì)算for j=2:nfor i=2:m-1U(j,i)=(1-2*a*k/hA2)*U(j-1,i)+a*k/hA2*U(j-1,i-1)+a*k/hA2*U(j-1,i+1);endend5. 將溫度隨時(shí)間變化情況用曲線表示flnuiro604020.6 o6. 給出3000、9000、15000三個(gè)時(shí)刻的溫度分布情況，對(duì)溫度隨時(shí)間變化規(guī)律做說明。T3000=100.0000 63.4362 34.2299 15.8021 7.464

8、1 7.4641 15.802134.2299 63.4362 100.0000T9000=100.0000 81.6930 65.6076 53.6839 47.3466 47.3466 53.683965.6076 81.6930 100.0000T15000=100.0000 89.9415 81.0962 74.5310 71.0378 71.0378 74.5310 81.0962 89.9415 100.0000根據(jù)數(shù)據(jù)分析，在同一個(gè) x點(diǎn)上溫度隨時(shí)間的增加而增加，但增幅變小。x-T圖形仍為拋物線，但隨著時(shí)間的增加，極值變小，圖像變得平緩。7. 用計(jì)算數(shù)據(jù)說明，并結(jié)合差分方程余項(xiàng)，

9、空間、時(shí)間間隔對(duì)求解精度影響。數(shù)據(jù)量較大，且原理相同，我取一個(gè)向量演示一下。保持空間間隔不變，修改時(shí)間間隔，時(shí)間間隔加大，得到的誤差加大。保持時(shí)間間隔不變，修改空間間隔，空間間隔加大，得到的誤差加大。修改空間間隔的誤差在增量比修改時(shí)間間隔的大。從方差余項(xiàng)上來看,Af d2u2 dt2(沒有公式編輯器。只能從ppt里粘貼了)這個(gè)余項(xiàng)里的 t , x都在分母上，所以與誤差成正比，且 x的次數(shù) 應(yīng)該是比 t高，故影響較大。8. 用計(jì)算數(shù)據(jù)說明，穩(wěn)定性要求對(duì)求解精度的影響。修改穩(wěn)定性，即修改 x和t分的份數(shù)(ngrid ),之后看誤差。穩(wěn)定性越高，解的精度越高。即在滿足穩(wěn)定性要求 (a* t/( xA

10、2)<0.5 )時(shí)，a* t/( xA2)越接近0,誤差越小。從概念上理解，穩(wěn)定性越好，對(duì)引入時(shí)間層誤差的抑制能力越強(qiáng)。所以誤差越小。二、調(diào)用MATLAB函數(shù)完成上述計(jì)算1.編寫M文件求解上述方程，并用適當(dāng)?shù)奈淖謱?duì)程序做出說明。function pdepediaoyongm=0;x=linspace(0,1,11);%x 的網(wǎng)格t=linspace(0,20000,101);%t 的網(wǎng)格sol = pdepe(m,pdefun,icfun,bcfun,x,t);% 調(diào)用函數(shù)T=sol(:,:,1);% 解figure;% 畫圖surf(x,t,T)xlabel('x')y

11、label('t')zlabel('T')dt=20000/100;%t 步長h3000=3000/dt;h9000=9000/dt;h15000=15000/dt;%3000,9000,15000下，溫度分別在 T矩陣的哪些行T3000=T(h3000,:)T9000=T(h9000,:)T15000=T(h15000,:)%輸出三個(gè)時(shí)間下的溫度分布%不再對(duì)三個(gè)時(shí)間下的溫度-長度曲線畫圖，其圖像就是三維圖的截面function c,f,s=pdefun(x,t,T,DuDx)%PDE 方程函數(shù)c=100000;f=DuDx;s=0;function u0=ic

12、fun(x)% 初始條件函數(shù)u0=0;function pl,ql,pr,qr=bcfun(xl,Tl,xr,Tr,t)% 邊界條件函數(shù)pl = Tl-100;ql = 0;pr = Tr-100;qr = 0;2.將溫度隨時(shí)間變化情況用曲線表示。12010080604020024X 103.給出3000、9000、15000三個(gè)時(shí)刻的溫度分布情況。T3000 =100.0000 67.1058 39.8611 21.1973 10.9885 7.8279 10.9885 21.1973 39.8611 67.1058 100.0000T9000=100.0000 83.4839 68.6032 56.8191 49.2705 46.6732 49.2705 56.8191 68.6032 83.4839 100.0000T15000=100.0000 90.8310 82.5601 75.9972 71.7845 70.3330 71.7845 75.9972 82.5601 90.8310 100.0000根據(jù)數(shù)據(jù)分析，在同一個(gè) x點(diǎn)上溫度隨時(shí)間的增加而增加，但增幅變小。 x-T圖形仍為拋物 線，但隨著時(shí)間的增加，極