對流擴(kuò)散方程

1、徐州工程學(xué)院卷課程設(shè) 計(jì)報(bào)告課程名稱偏微分方程數(shù)值解課題名稱對流擴(kuò)散方程的迎風(fēng)格式的推導(dǎo)和求解專業(yè)信息與計(jì)算科學(xué)班級10信計(jì)3姓名學(xué)號指導(dǎo)教師楊 揚(yáng)2021年 5 月 23日、實(shí)驗(yàn)?zāi)康?進(jìn)一步穩(wěn)固理論學(xué)習(xí)的結(jié)果,學(xué)習(xí)雙曲型對流擴(kuò)散方程的迎風(fēng)格式的構(gòu)造 方法,以及穩(wěn)定的條件.從而進(jìn)一步了解差分求解偏微分方程的一些根本概念, 掌握數(shù)值求解偏微分方程的根本過程.在此根底上考慮如何使用Matlab的軟件進(jìn)行上機(jī)實(shí)現(xiàn),并針對具體的題目給出相應(yīng)的數(shù)值計(jì)算結(jié)果.二、實(shí)驗(yàn)題目：u t a1u Y b1u YY f ,0 x 1, 0 t 1txxxu ( x ,0) exp( x / 2);u (0 , t)

2、 exp( t); u (1, t) exp( 1 / 2 t);其中 a1=1,b1=2, f (x, t) exp( x / 2 t) o用迎風(fēng)格式求解雙曲型對流擴(kuò)散方程,觀差分解對真解的斂散性(u exp(x/2 t)三、實(shí)驗(yàn)原理:1、用迎風(fēng)格式求解雙曲型對流擴(kuò)散方程,迎風(fēng)格式為:n 1 nn nu j u j u j u jjL a1 Lhn 1 nnu j u ju i 1 uL a1 j hnu j 1b1 j2u； u；1 h22u； u；1 h2假設(shè)令 r a1 * / h, b1 * / h2,那么迎風(fēng)格式可整理為：n 1nnn%(1 r 2 )% (r )% 1、in 1n

3、nnu j (1 r 2 )u j ( r)% 1% 12、穩(wěn)定條件：h 2 /( 2 * b 1 a1 * h ),( a1 0)o,、(*h 2 /( 2 * b 1 a1 * h ),( a1 0)fjn (a10)fjn (a10)fjn (a10)fjn (a10)四、數(shù)值實(shí)驗(yàn)的過程、相關(guān)程序及結(jié)果：本次的實(shí)驗(yàn)題目所給出的邊界條件是第一邊界條件, 直接利用所給的邊界條 件,我們可以給出界點(diǎn)處以及第 0層的函數(shù)值,根據(jù)a1的正負(fù)性,使用相應(yīng)的 1或者2式,求出其他層的函數(shù)值.誤差轉(zhuǎn)化成圖的形式,并輸出最大值. 針對三種不同的輸入對應(yīng)輸出結(jié)果：A: a1=1;b1=2;a=1;b=1;h

4、=0.1;k=0.001;結(jié)果一：1 .誤差最大值：e =7.9402e-0042 .誤差圖如下列圖所示:B: a1=-1;b1=2;a=1;b=1;h=0.1;k=0.001;結(jié)果二：1 .誤差最大值：e =0.06822 .誤差圖：C: a1=-1;b1=-0.1;a=1;b=1;h=0.1;k=0.001;結(jié)果三：1 .誤差最大值：6.2221e+0052 .誤差圖：五、實(shí)驗(yàn)結(jié)論：通過上機(jī)實(shí)現(xiàn),進(jìn)一步直觀了解流擴(kuò)散方程的穩(wěn)定具有很強(qiáng)的條件性,只要在a1,b1,h和 滿足(*)式時(shí)才是穩(wěn)定的,如結(jié)果一、二,否那么會(huì)出現(xiàn)結(jié)果三的情形,誤差相當(dāng)大. 本次實(shí)驗(yàn),熟悉并掌握了差分格式的一般構(gòu)造方法

5、,理清了具體的步驟,提升了利用計(jì)算機(jī)解決問題的水平.附：Matlab源代碼：1 . function z=ft(x)% 求下邊界z=exp(x/2);2 .function z=fx1(t)% 求左邊界z=exp(-t);3 .function z=fx2(t)% 求右邊界z=exp(1/2-t);4 .function z=f(x,t)% 求右端函數(shù)z=-exp(x/2-t);5 .function z=fu(x,t)% 求真解 z=exp(x/2-t);6 . function X,T,z=upwindL(a1,b1,a,b,h,k)% 用迎風(fēng)格式求解upwindL(1,2,1,1,0.1

6、,0.1)x=0:h:a;t=0:k:b;T,X=meshgrid(t,x);m=length(x);n=length(t);r1=a1*k/h;r2=b1*k/hA2;uu=zeros(m,n);%儲(chǔ)存數(shù)值解z=uu;%儲(chǔ)存誤差for i=1:m%求下邊界endfor j=2:n%求左右邊界uu(1,j)=fx1(t(j);uu(m,j)=fx2(t(j); end %迎風(fēng)格式求內(nèi)點(diǎn),從下往上 if(a1>0) for j=2:n for i=2:m-1%從左往右uu(i,j)=(1-r1-2*r2)*uu(i,j-1)+(r1+r2)*uu(i-1,j-1)+r2*uu(i+1,j-1)+k*f(x(i),t(j-1);%求數(shù)值解z(i,j)=abs(uu(i,j)-fu(x(i),t(j);% 求誤差 end end else for j=2:nfor i=2:m-1%從左往右uu(i,j)=(1+r1-2*r2)*uu(i,j-1)+(r2-r1)*uu(i+1,j-1)+r2*uu(i-1,j-1)+k*f(x(i),t(j-1);%求數(shù)值解z(i,j)=abs(uu(i,j)-fu(x(i),t(j);% 求誤差 end end end %主函數(shù),用于輸出 7. X,T,z=upwindL(a1,b1,a,b,h,k); me