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

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

2、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);%計算所調(diào)用的函數(shù)x,t=meshgrid(x,t);mesh(x,t,T);%畫圖，并且把坐標(biāo)軸名稱改為x，t，Txlabel('x')ylabel('t')zlabel('T')T%輸出溫度矩陣dt=tspan(2)/ngrid(1);%t

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

4、 fprintf('n%sn','沒有穩(wěn)定性') errorend%真實值計算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%輸出溫度矩陣%誤差計算jmax=1/dx+1;%網(wǎng)格點數(shù)rms=wuchajisuan(T,Te,jmax)rms%輸出誤差function r

5、ms=wuchajisuan(T,Te,jmax)for j=1:jmax rms=(T(j)-Te(j)2/jmax)(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:n for i=2:m-1 for g=1:m-1Ue(j,i)=100-(400/(2*

6、g-1)/pi)*sin(2*g-1)*pi*xe(j)*exp(-a*(2*g-1)2*pi2*te(i) end endendfunction 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),tspan(2),n);%畫網(wǎng)格U=zeros(ngrid);U(:,1)

7、=g1(t);%邊界條件U(:,m)=g2(t);U(1,:)=f(x);%初值條件 %差分計算 for j=2:n for i=2:m-1 U(j,i)=(1-2*a*k/h2)*U(j-1,i)+a*k/h2*U(j-1,i-1)+a*k/h2*U(j-1,i+1); endend5. 將溫度隨時間變化情況用曲線表示6. 給出3000、9000、15000三個時刻的溫度分布情況，對溫度隨時間變化規(guī)律做說明。T3000=100.0000 63.4362 34.2299 15.8021 7.4641 7.4641 15.8021 34.2299 63.4362 100.0000T9000=10

8、0.0000 81.6930 65.6076 53.6839 47.3466 47.3466 53.6839 65.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ù)分析，在同一個x點上溫度隨時間的增加而增加，但增幅變小。x-T圖形仍為拋物線，但隨著時間的增加，極值變小，圖像變得平緩。7. 用計算數(shù)據(jù)說明，并結(jié)合差分方程余項，空間、時間間隔對求解精度影響。數(shù)據(jù)量較大，且原理相同，我取一個向量演示一下。保持空間間隔不變，修

9、改時間間隔，時間間隔加大，得到的誤差加大。保持時間間隔不變，修改空間間隔，空間間隔加大，得到的誤差加大。修改空間間隔的誤差在增量比修改時間間隔的大。從方差余項上來看，（沒有公式編輯器。只能從ppt里粘貼了）這個余項里的t，x都在分母上，所以與誤差成正比，且x的次數(shù)應(yīng)該是比t高，故影響較大。8. 用計算數(shù)據(jù)說明，穩(wěn)定性要求對求解精度的影響。修改穩(wěn)定性，即修改x和t分的份數(shù)（ngrid），之后看誤差。穩(wěn)定性越高，解的精度越高。即在滿足穩(wěn)定性要求（a*t/(x2)<0.5）時，a*t/(x2)越接近0，誤差越小。從概念上理解，穩(wěn)定性越好，對引入時間層誤差的抑制能力越強。所以誤差越小。二、調(diào)用M

10、ATLAB函數(shù)完成上述計算1. 編寫M文件求解上述方程，并用適當(dāng)?shù)奈淖謱Τ绦蜃龀稣f明。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')ylabel('t')zlabel('T')dt=20000/100;%t步長h3000=3000/dt;h9000=9000/

11、dt;h15000=15000/dt;%3000,9000,15000下，溫度分別在T矩陣的哪些行T3000=T(h3000,:)T9000=T(h9000,:)T15000=T(h15000,:)%輸出三個時間下的溫度分布%不再對三個時間下的溫度-長度曲線畫圖，其圖像就是三維圖的截面function c,f,s=pdefun(x,t,T,DuDx)%PDE方程函數(shù)c=100000;f=DuDx;s=0;function u0=icfun(x)%初始條件函數(shù)u0=0;function pl,ql,pr,qr=bcfun(xl,Tl,xr,Tr,t)%邊界條件函數(shù)pl = Tl-100;ql =

12、 0;pr = Tr-100;qr = 0;2. 將溫度隨時間變化情況用曲線表示。3. 給出3000、9000、15000三個時刻的溫度分布情況。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ù)分析，在同一個x點上溫度隨時間的增加而增加，但增幅變小。x-T圖形仍為拋物線，但隨著時間的增加，極值變小，圖