數(shù)學(xué)建模教材20第二十章偏微分方程的數(shù)值解

第二十章偏微分方程的數(shù)值解自然科學(xué)與工程技術(shù)中種種運(yùn)動(dòng)發(fā)展過(guò)程與平衡現(xiàn)象各自遵守一定的規(guī)律。這些規(guī)律的定量表述一般地呈現(xiàn)為關(guān)于含有未知函數(shù)及其導(dǎo)數(shù)的方程。我們將只含有未知多元函數(shù)及其偏導(dǎo)數(shù)的方程，稱之為偏微分方程。方程中出現(xiàn)的未知函數(shù)偏導(dǎo)數(shù)的最高階數(shù)稱為偏微分方程的階。如果方程中對(duì)于未知函數(shù)和它的所有偏導(dǎo)數(shù)都是線性的，這樣的方程稱為線性偏微分方程，否則稱它為非線性偏微分方程。初始條件和邊界條件稱為定解條件，未附加定解條件的偏微分方程稱為泛定方程。對(duì)于一個(gè)具體的問(wèn)題，定解條件與泛定方程總是同時(shí)提出。定解條件與泛定方程作為一個(gè)整體，稱為定解問(wèn)題?！?偏微分方程的定解問(wèn)題各種物理性質(zhì)的定常(即不隨時(shí)間變化)過(guò)程，都可用橢圓型方程來(lái)描述。其最典型、最簡(jiǎn)單的形式是泊松(Poisson)方程(1)特別地，當(dāng)f(x,y)三0時(shí)，即為拉普拉斯(Laplace)方程，又稱為調(diào)和方程S2u S2u八(2)——+——=0(2)Sx2 Sy2帶有穩(wěn)定熱源或內(nèi)部無(wú)熱源的穩(wěn)定溫度場(chǎng)的溫度分布，不可壓縮流體的穩(wěn)定無(wú)旋流動(dòng)及靜電場(chǎng)的電勢(shì)等均滿足這類方程。Poisson方程的第一邊值問(wèn)題為If^u-+S2u=f(x,y) (x,y)￡△.SX2Sy2 (3)%(x,y)I =3(x,y) r=Sa▼ (x,y)e「其中a為以r為邊界的有界區(qū)域，r為分段光滑曲線，aUr稱為定解區(qū)域，f(x,y),我x,y)分別為Ar上的已知連續(xù)函數(shù)。第二類和第三類邊界條件可統(tǒng)一表示成TOC\o"1-5"\h\z\o"CurrentDocument"包口? 、 一\o"CurrentDocument"口+au口 =我x,y) (4)OSn 口(x,y)er其中n為邊界r的外法線方向。當(dāng)a=0時(shí)為第二類邊界條件，aw0時(shí)為第三類邊界條件。在研究熱傳導(dǎo)過(guò)程，氣體擴(kuò)散現(xiàn)象及電磁場(chǎng)的傳播等隨時(shí)間變化的非定常物理問(wèn)題時(shí)，常常會(huì)遇到拋物型方程。其最簡(jiǎn)單的形式為一維熱傳導(dǎo)方程Su「S2u_———a =0 (a>0) (5)St Sx2方程(5)可以有兩種不同類型的定解問(wèn)題：初值問(wèn)題(也稱為Cauchy問(wèn)題)-360-

蜀u S2u八余—-a =0?蜀u S2u八余—-a =0?31 Sx2加(X,0)=①(X)初邊值問(wèn)題蜀u S2ua10\o"CurrentDocument"St SX2t>0,一8<X<+8(6)u(x,0)=①(x)*u(0,t)=g(t),u(1,t)=g(t),0<t<T(7)其中我X),g?),g2(t)為已知函數(shù)，且滿足連接條件①(0)=g1(0),問(wèn)題(7)中的邊界條件u(0,t)第三類邊界條件為叭1)=g2(0)g1(t'u(1,t)=g2(t)稱為第一類邊界條件。第二類和x=0其中入(t)>0入(t)>0。當(dāng)入(t)=X(t)三0時(shí)，為第二類邊界條件，否則稱為第三類邊界1條件。2 1 2Su Su=0市+aSX物理中常見(jiàn)的一維振動(dòng)與波動(dòng)問(wèn)題可用二階波動(dòng)方程S2Su Su=0市+aSX物理中常見(jiàn)的一維振動(dòng)與波動(dòng)問(wèn)題可用二階波動(dòng)方程S2u S2u=a2(9)St2(10)描述，它是雙曲型方程的典型形式。方程(10)的初值問(wèn)題為*——=a2——St2 SX2U(x,0)=我x)t>0,-8<X<+8一8<X<+8(11)=。(X)t=0一8<X<+8邊界條件一般也有三類，最簡(jiǎn)單的初邊值問(wèn)題為-361-d2d2ut>0,0<x<l—=。(x)0<x<latt=0u(l,t)=g(t)0<t<T如果偏微分方程定解問(wèn)題的解存在，唯一且連續(xù)依賴于定解數(shù)據(jù)（即出現(xiàn)在方程和定解條件中的已知函數(shù)），則此定解問(wèn)題是適定的。可以證明，上面所舉各種定解問(wèn)題都是適定的?！?偏微分方程的差分解法差分方法又稱為有限差分方法或網(wǎng)格法，是求偏微分方程定解問(wèn)題的數(shù)值解中應(yīng)用最廣泛的方法之一。它的基本思想是：先對(duì)求解區(qū)域作網(wǎng)格剖分，將自變量的連續(xù)變化區(qū)域用有限離散點(diǎn)（網(wǎng)格點(diǎn)）集代替；將問(wèn)題中出現(xiàn)的連續(xù)變量的函數(shù)用定義在網(wǎng)格點(diǎn)上離散變量的函數(shù)代替；通過(guò)用網(wǎng)格點(diǎn)上函數(shù)的差商代替導(dǎo)數(shù)，將含連續(xù)變量的偏微分方程定解問(wèn)題化成只含有限個(gè)未知數(shù)的代數(shù)方程組（稱為差分格式）。如果差分格式有解，且當(dāng)網(wǎng)格無(wú)限變小時(shí)其解收斂于原微分方程定解問(wèn)題的解，則差分格式的解就作為原問(wèn)題的近似解（數(shù)值解）。因此，用差分方法求偏微分方程定解問(wèn)題一般需要解決以下問(wèn)題：（i）選取網(wǎng)格；（ii）對(duì)微分方程及定解條件選擇差分近似，列出差分格式；（iii）求解差分格式；（iv）討論差分格式解對(duì)于微分方程解的收斂性及誤差估計(jì)。下面我們只對(duì)偏微分方程的差分解法作一簡(jiǎn)要的介紹。2.1橢圓型方程第一邊值問(wèn)題的差分解法以Poisson方程（1）為基本模型討論第一邊值問(wèn)題的差分方法?？紤]Poisson方程的第一邊值問(wèn)題（3）蜀2u"2u^r-+—=f（x,y） （x,y）￡△.ax2ay2%（x,y）I =我x,y） 「=Sa▼ （x,y）e「取h,T分別為x方向和y方向的步長(zhǎng)，以兩族平行線x=xk=kh,y=y,=j（k,j=0,±1比2,L）將定解區(qū)域剖分成矩形網(wǎng)格。節(jié)點(diǎn)的全球記為R=｛（xk,y/Ixk=kh,y=j,i,j?為整數(shù)｝。定解區(qū)域內(nèi)部的節(jié)點(diǎn)稱為內(nèi)點(diǎn)，記內(nèi)點(diǎn)集rI△為a。邊界r與網(wǎng)格線的交點(diǎn)稱為邊界點(diǎn)，邊界點(diǎn)全體記為r。與節(jié)點(diǎn)hT hT（xk,yj沿x方向或y方向只差一個(gè)步長(zhǎng)的點(diǎn)（xk±1,yj和（X卜,y±±1）稱為節(jié)點(diǎn)（xk,1）的相鄰節(jié)點(diǎn)。如果一個(gè)內(nèi)點(diǎn)的四個(gè)相鄰節(jié)點(diǎn)均屬于AUr，稱為正則內(nèi)點(diǎn)，正則內(nèi)點(diǎn)的全體記為a（1），至少有一個(gè)相鄰節(jié)點(diǎn)不屬于aUr的內(nèi)點(diǎn)稱為非正則內(nèi)點(diǎn)，非正則內(nèi)點(diǎn)的全體記為A（2）。我們的問(wèn)題是要求出問(wèn)題（3）在全體內(nèi)點(diǎn)上的數(shù)值解。為簡(jiǎn)便記，記（k,j）=（xk,yj,u（k,j）=u（xk,y）,fkj=f（x^,yj。對(duì)正則內(nèi)點(diǎn)-362-

(k,j)￡A(D，由二階中心差商公式S2ud.x2=u(k+1,j)-2u(k,j)+u(k-1,j)+O(h2S2ud.x2+O(T2)u(k,j+1)2u(k,j)+u+O(T2)T2Poisson方程(1)在點(diǎn)(k,j)處可表示為u(k+1,j)一2u(k,j)+u(k-1,j) u(k,j+1)一2u(k,j)+u(k,j一1)TOC\o"1-5"\h\z廿一十 T2 (12)=f+O(h2+T2)在式(12)中略去O(h2+T2),即得與方程(1)相近似的差分方程=fk,j(13)\o"CurrentDocument"u—2u +uu—2u +u=fk,j(13)-k+1,j kk,j k—1,j+-k,j+1 kTj k,j—1h2 T2式(13)中方程的個(gè)數(shù)等于正則內(nèi)點(diǎn)的個(gè)數(shù)，而未知數(shù)ukj則除了包含正則內(nèi)點(diǎn)處解u的近似值，還包含一些非正則內(nèi)點(diǎn)處u的近似值，因而方程個(gè)數(shù)少于未知數(shù)個(gè)數(shù)。在非正則內(nèi)點(diǎn)處Poisson方程的差分近似不能按式(13)給出，需要利用邊界條件得到。邊界條件的處理可以有各種方案，下面介紹較簡(jiǎn)單的兩種。(1)直接轉(zhuǎn)移。用最接近非正則內(nèi)點(diǎn)的邊界點(diǎn)上的u值作為該點(diǎn)上u值的近似，這就是邊界條件的直接轉(zhuǎn)移。例如，點(diǎn)尸(k,j)為非正則內(nèi)點(diǎn)，其最接近的邊界點(diǎn)為Q點(diǎn)，則有：ukj=u(Q)=我Q)，(k,j)￡A(2)(ii)線性插值。這種方案是通過(guò)用同一條網(wǎng)格線上與點(diǎn)尸相鄰的邊界點(diǎn)R與內(nèi)點(diǎn)T作線性插值得到非正則內(nèi)點(diǎn)P(k,j)處u值的近似。由點(diǎn)R與T的線性插值確定u(P)的近似值為u=JL-我R)+du(T)k,jh+dh+d其中d二Rp|，h=|PT|，其截?cái)嗾`差為O(h2)。由式(13)所給出的差分格式稱為五點(diǎn)菱形格式，實(shí)際計(jì)算時(shí)經(jīng)常取h=T，此時(shí)五點(diǎn)菱形格式可化為23k+1,j+uk_1,j+uk,j+1+uk,j—1—4uk，尸f,j (14)簡(jiǎn)記為(15)1Au(15)+u―1,j+uk+u―1,j+uk,j+1+uk,j—1—4uk,j其中0u「uk+1,j求解差分方程組最常用的方法是同步迭代法，同步迭代法是最簡(jiǎn)單的迭代方式。除邊界節(jié)點(diǎn)外，區(qū)域內(nèi)節(jié)點(diǎn)的初始值是任意取定的。例1用五點(diǎn)菱形格式求解Laplace方程第一邊值問(wèn)題-363-

中(x,y)I =ig[(i+x)2+y2]▼ (x,y)e「 1其中△={(X,y)I0<X,y<1}。取h=T=-o(x,y)￡(x,y)￡△r=3A圖1網(wǎng)格劃分圖解節(jié)點(diǎn)編號(hào)為(k,j),k=0,1,2,3,j=0,1,2,3。網(wǎng)格中有四個(gè)內(nèi)點(diǎn)，均為正則內(nèi)點(diǎn)。由五點(diǎn)菱形格式，得方程組：(u+u+u+u—4u)=0(16)(u +u+u+u—4u)=0(16)(u1,3+u2,20,2—4u21)=0(u2,3+u3,2+u2,1+u1,2代入邊界條件u,u1,0F4，，，，<,u,1—40,1,u0,210,u3,1,u3,2的值3,2—41818ui,解非齊次線性方程組求得u1=0.6348,計(jì)算的Matlab程序如下:爭(zhēng)+u1,0 0,18u+u288=—皿+u2,好2,182,2ff1,3<u2,3(16)式可以化成/0,28+uf3,2fu12=1.06,u21=0.7985,u22=1.1698(17)clc,clearf1=@(x)2*log(1+x);f2=@(x)10g((1+x)J2+1);f3=@(y)1og(1+y.A2);f4=@(y)1og(4+y.A2);u=zeros(4);m=4;n=4;h=1/3;u(1,1:m)=feva1(f3,0:h:(m-1)*h)';u(n,1:m)=feva1(f4,0:h:(m-1)*h)';u(1:n,1)=feval(f1,0:h:(n-1)*h);u(1:n,m)=feva1(f2,0:h:(n-1)*h);b=-[u(2,1)+u(1,2);u(4,2)+u(3,1);u(2,4)+u(1,3);u(3,4)+u(4,3)];a=[-4110;1-401;10-41;011-4];-364-x=a\b實(shí)際上，可以使用同步迭代法、異步迭代法、逐次超松弛迭代法等方法求方程組(16)的解。同步迭代法的迭代格式為1rU(i+1)= [U(i)+U(i)+U(i)+U(i)]k,j 4k-1,j k+1,j k,j-1 k,j+1異步迭代法的迭代格式為1r 1U(i+1)=[U(i+1)+U(i)+U(i+1)+U(i)]k,j 4k-1,j k+1,j k,j-1 k,j+1由于在異步迭代法中有一半是用了迭代的新值，所以可以預(yù)料異步迭代法的收斂速度比同步迭代法的收斂速度要快一倍左右。下面我們用逐次超松弛迭代法求例1的數(shù)值解。程序如下(以下所有的程序放在一個(gè)文件中)：functionsol=mainf1=@(x)2*log(1+x);f2=@(x)log((1+x).A2+1);f3=@(y)10g(1+y.A2);f4=@(y)10g(4+y.A2);a=1；b=1;h=1/3;to1=0.00001;max1=1000;so1=dirich(f1,f2,f3,f4,a,b,h,to1,max1)functionU=dirich(f1,f2,f3,f4,a,b,h,to1,max1)%Input-f1,f2,f3,f4areboundaryfunctionsinputasstrings% -aandbrightendpointsof[0,a]and[0,b]% -hstepsize% -to1istheto1erance%Output-Uso1utionmatrix;%Iff1,f2,f3andf4areM-fi1efunctionsca11U=dirich(@f1,@f2,@f3,@f4,a,b,h,to1,max1).%iff1,f2,f3andf4areanonymousfunctionsca11U=dirich(f1,f2,f3,f4,a,b,h,to1,max1).%Initia1izeparametersandUn=fix(a/h)+1;m=fix(b/h)+1;ave=(a*(feva1(f1,0)+feva1(f2,0))...+b*(feva1(f3,0)+feva1(f4,0)))/(2*a+2*b);U=ave*ones(n,m);%BoundaryconditionsU(1,1:m)=feva1(f3,0:h:(m-1)*h)';U(n,1:m)=feva1(f4,0:h:(m-1)*h)';U(1:n,1)=feva1(f1,0:h:(n-1)*h);U(1:n,m)=feva1(f2,0:h:(n-1)*h);%逐次超松弛迭代法(SOR)參數(shù)w=4/(2+sqrt(4-(cos(pi/(n-1))+cos(pi/(m-1)))A2));%Refineapproximationsandsweepoperatorthroughoutthegrid-365-

err=1;cnt=0;while((err>tol)&(cnt<=max1))err=0;forj=2:m-1fori=2:n-1relx=w*(U(i,j+1)+U(i,j-1)+U(i+1,j)+U(i-1,j)-4*U(i,j))/4;U(i,j)=U(i,j)+relx;err=abs(relx);endendcnt=cnt+1;endU=flipud(U');上述程序中的函數(shù)dirich（f1,f2,f3,f4,a,b,h,tol,max1）可以求解如下類型的問(wèn)題:d2ud2u + =0,(x,y)ea={(%,y):0<x<a,0<y<b}d.x2 dy2而且滿足條件u(x,0)=f(x)，u(x,b)=f(x)，0<x<au(0,y)=f3(y)，u(a,y)=f4(y)，0<y<b上述程序中步長(zhǎng)h和a,b之間的關(guān)系為，存在整數(shù)n和m，使得a=nh，b=mh。當(dāng)h=T時(shí)，利用點(diǎn)（k,j）,（k土1,j土1）構(gòu)造的差分格式白(u+1川+u"1+u_1,j+1+u-1,j-1-4uk廣f,j稱為五點(diǎn)矩形格式，簡(jiǎn)記為_(kāi)u-fh2□k,j Jk,j其中口ukj=uk+1,j+1+uk+1,j-1+uk-1,j+1+uk-1,j-1-4uk,j。2.2'拋物型方程的差分解法以一維熱傳導(dǎo)方程(5)(a>0)TOC\o"1-5"\h\z\o"CurrentDocument"du d2u 八(a>0)——-a =0dt dt2為基本模型討論適用于拋物型方程定解問(wèn)題的幾種差分格式。首先對(duì)xt平面進(jìn)行網(wǎng)格剖分。分別取h,T為x方向與t方向的步長(zhǎng)，用兩族平行直線x=xk=kh(k=0,±1比2,L)，t=tj=jt(j=0,1,2JL)，將xt平面剖分成矩形網(wǎng)格，節(jié)點(diǎn)為(xk,tj)(k=0,±1比2,L,j=0,1,2L)。為簡(jiǎn)便起見(jiàn)，記(k,j)=(xk,yj)，u(k,j)= u(xk,yj)，Q=3(xk)， g1廣 g1(tj)， g2廣g2(tj)，入］上="(tj)，入2jj，。2.2.1微分方程的差分近似-366-

TOC\o"1-5"\h\z一.、 du d2u在網(wǎng)格內(nèi)點(diǎn)(k,j)處，對(duì)不分別采用向前、向后及中心差商公式，對(duì)一采用二dt dx2階中心差商公式，一維熱傳導(dǎo)方程(5)可分別表示為u(k,j+1)-u(k,j)u(k+1,j)-2u(k,j)±u(k-1,j)=O(t+h2)T -a h2u(k+1,j)-2u(k,j)+u(k-1,j)=O(T+h2)-a h2 -u(k,j+1)-u(k,j-1)u(k+1,j)-2u(k,j)+u(k-1,j)=O(T+h2)2T -a h2由此得到一維熱傳導(dǎo)方程的不同的差分近似u-uu_2u+uk,j+1 k,j-a k+1,j k/k-1,j=QT h2u-uu-2u+uk,j k,j-1-ak+1,jk,j k-1,j_qT h2u-uu-2u+uk,j+1 k,j-1-a k+1,j k/ k-1,j=Q2T h2初、邊值條件的處理為用差分方程求解定解問(wèn)題(6)，(7)等，還需對(duì)定解條件進(jìn)行離散化。對(duì)初始條件及第一類邊界條件，可直接得到u,q=uu,q=u(X，Q)=Qu —■u(Q,t)—gu —u(l,t)—gln,jT ]其中n―-,m——。hT(k=Q,±1,L或k=Q,1JL,n)(j=Q,1L,m-1)對(duì)第二、三類邊界條件則需用差商近似。下面介紹兩種較簡(jiǎn)單的處理方法。du(i)在左邊界(x=Q)處用向前差商近似偏導(dǎo)數(shù)一，在右邊界(x=l)處用向后差dx商近似偏導(dǎo)數(shù)生，即dxdudx(j—Q,1L,m(j—Q,1L,m)u(n,j)-u(n-6+O(h)l6dx(n,j)即得邊界條件(8)的差分近似為-367-

熱-u.1,j 0,j-九u二p(j=0,1L,m)(25)今h1(j=0,1L,m)(25)J-u、n-4 n-1,j+兒u=g.h 2jn,j 2jdu(ii)用中心差商近似上，即dxdud.xu(1,j)dud.xu(1,j)—u(—1,j)+O(h2)(0,j)dudxu(n+1,j)—u(n-1,j)+O(h2)(j=0,1L,m)(n,j)則得邊界條件的差分近似為和-u1,j -1,j -Xugo(j=0,1L,m)(26).2h 1j0,(j=0,1L,m)(26).n+Lj n-1,j+Xugg. 2h 2jn,j 2j這樣處理邊界條件，誤差的階數(shù)提高了，但式（26）中出現(xiàn)定解區(qū)域外的節(jié)點(diǎn)（-1,j）和（n+1,j），這就需要將解拓展到定解區(qū)域外?？梢酝ㄟ^(guò)用內(nèi)節(jié)點(diǎn)上的u值插值求出u-1,j和un+1j，也可以假定熱傳導(dǎo)方程（5）在邊界上也成立，將差分方程擴(kuò)展到邊界節(jié)點(diǎn)上，由,此消去u-1j和un+1.。幾種常用的差分格式下面我們以熱傳導(dǎo)方程的初邊值問(wèn)題⑺為例給出幾種常用的差分格式。⑴古典顯式格式aT為便于計(jì)算，令r=——，式（20）改寫(xiě)成以下形式h2晨,j+1gT+1,j+（1-2r）uk：T-1,j將式（20）與（23），（24）結(jié)合，我們得到求解問(wèn)題（7）的一種差分格式.gru...+（1-2r）u +ru... （k=1,2JL,n-1,j=0,1J_,m-1）.k,j+1 k+1,j ' ）k,j k-1,j第二3 （kg1,2L,n） （27）.k,0k"0,jgg1j，“n，廣g2j （jg1,2L,m）由于第0層（jg0）上節(jié)點(diǎn)處的u值已知（u=中），由式（25）即可算出u在第一層k,0k（jg1）上節(jié)點(diǎn)處的近似值u 。重復(fù)使用式（25），可以逐層計(jì)算出各層節(jié)點(diǎn)的近似值，k,1因此此差分格式稱為古典顯示格式。又因式中只出現(xiàn)相鄰兩個(gè)時(shí)間層的節(jié)點(diǎn)，故此式是二層顯式格式。（ii）古典隱式格式將（21）整理并與式（23），（24）聯(lián)立，得差分格式如下-368-

* =u+r(u—2u +u )(k=1,2J_,n—1,j=0,1L,m-1)余k,j+1 k,j k+1,j+1 k,j+1 k—1,j+1第二3 (k=0,1L,n) (28)▲k,0 *k3=g,u.=g (j=0,1L,m)▼0,j1jn,j2j其中r=上。雖然第0層上的u值仍為已知，但不能由式(28)直接計(jì)算以上各層節(jié)h2點(diǎn)上的值u ，必須通過(guò)解下列線性方程組k,jT+2r，，—r，，，h2點(diǎn)上的值u ，必須通過(guò)解下列線性方程組k,jT+2r，，—r，，，<—r1+2rO—rO—r—r/丫u+rgn—2,j+1888/Tu81,j+1—r1+2r888=88，，，，1,j1j+1u2,jMn—1,j+1un—2,j+rgn—1,j/888882j+1才能由uk3嚴(yán)算對(duì)，川，故此差分方程稱為古典隱式格式。此方程組是三對(duì)角方程組，(iii)杜福特一弗蘭克爾(DoFort—Frankel)格式DoFort—Frankel格式是三層顯式格式，它是由式u—uk,j+12ru—uk,j+12rk,j-1—a^k+J k,j+1h2k,j—1 k―j與初始條件及第一類邊界條件結(jié)合得到的。具體形式如下:2r1—22r*k,0▲*k,0▲▲u▲0,j▼(u +u)+u1+2r k+1,j k-1,j 1+2r k,j-1=Wk=g,u=g1j n,j 2j,(k=1,2,L,n—1,j=1,2,L,m—1)(k=0,1L,n) (29)(j=0,1L,m)用這種格式求解時(shí)，除了第0層上的值uk0由初值條件得到，必須先用二層格式求出第1層上的值uk1，然后再按格式(29)逐層計(jì)算uk(j=2,3,L,m)。TOC\o"1-5"\h\z例2用古典顯示格式求初邊值問(wèn)題。 "蜀u d2u▲——— =0 0<t<3,0<x<3￡t dx2*(x,0)=x2 0<x<3▲(0,t)=0, u(3,t)=9 0<t<3▼的數(shù)值解，取h=1,r=0.5?！?ar解：這里：a=1,r=—=0.5,我x)=x2,g(t)=0,g(t)=9。h2 1 2由格式-369-

* =ru +(1-2r)u余k* =ru +(1-2r)u余k,j+1 k+1,j k,jk,0可得到0,j+ru ,(k=1,2,L,n-1;j=0,1L,m-1)(k=0,1L,n)(j=0,1L,m)熱.=0.5(u余物k,050,j=X2k=0,k+1,j+"k-1,ju=93,j(k=1,2j=0,1,L,5)(k=0,1,2,3)(j=0,1,L,6)將初值u代入上式，即可算出：k,0u11=0.5x(u20+u00)=0.5x(4+0)=2u=0.5x(u+u)=0.5x(9+1)=5將邊界條件u01=0,u31=9及上述結(jié)果代入又可求得u12=2.5，u=5.5如此逐層計(jì)算，得全部節(jié)點(diǎn)上的數(shù)值解為u13=2.75，u=5.75，…2.2.4求解拋物形方程的Matlab程序（1）前向差分法求解熱傳導(dǎo)方程的MATLAB程序設(shè)u(X,0)=f(x)，其中，0<X<a，而且，u(0,t)=q,u(a,t)=c2，其中0<t<b，求解區(qū)間R={(x,t):0<x<a,0<t<b}內(nèi)u(x,t)=c2u(x,t)的近似解。t XXfunctionU=forwdif(f,c1,c2,a,b,c,n,m);%Input-f=u(x,0)asastring'f'% - c1=u(0,t)andc2=u(a,t)% - a and brightendpointsof [0,a]and [0,b]% - c the constantintheheat equation% - n and mnumberofgridpointsover[0,a] and[0,b]%Output-Usolutionmatrix;%InitializeparametersandUh=a/(n-1);k=b/(m-1);r=c人2*k/h人2;s=1-2*r;U=zeros(n,m);%BoundaryconditionsU(1,1:m)=c1;U(n,1:m)=c2;%GeneratefirstrowU(2:n-1,1)=feval(f,h:h:(n-2)*h)';%GenerateremainingrowsofUforj=2:m-370-

fori=2:n-1U(i,j)=s*U(i,j-1)+r*(U(i-1,j-1)+U(i+1,j-1));endendU=U';利用上述函數(shù)，我們可以編寫(xiě)如下程序計(jì)算例2的數(shù)值解。c=1;a=3;b=3;c1=0;c2=9;n=4;m=7;f=@(x)x.A2;sol=forwdif(f,c1,c2,a,b,c,n,m)（2）Crank-Nicholson求解熱傳導(dǎo)方程的MATLAB程序設(shè)u(x,0)=f(x)，其中，0<x<a，而且，u(0,t)=c1,u(a,t)=c2，其中0<t<b，求解區(qū)間R={(x,t):0<x<a,0<t<b}內(nèi)u(x,t)=c2u(x,t)的近似解。t xxfunctionU=crnich(f,c1,c2,a,b,c,n,m)%Input-f=u(x,0)% - c1=u(0,t)andc2=u(a,t)% - a and brightendpointsof [0,a]and [0,b]% - c the constantintheheat equation% - n and mnumberofgridpointsover[0,a] and[0,b]%Output-Usolutionmatrix;%Iffisan

%IffisanM-filefunctioncallU=crnich(@f,c1,c2,a,b,c,n,m).anonymousfunctioncallU=crnich(f,c1,c2,a,b,c,n,m).%Iffisan

%Iffisan%InitializeparametersandUh=a/(n-1);k=b/(m-1);r=c人2*k/h人2;s1=2+2/r;s2=2/r-2;U=zeros(n,m);%BoundaryconditionsU(1,1:m)=c1;U(n,1:m)=c2;%GeneratefirstrowU(2:n-1,1)=f(h:h:(n-2)*h)';%Formthediagonalandoff-diagonalelemntsofAand%theconstantvectorBandsolvetridiagonalsystemAX=BVd(1,1:n)=s1*ones(1,n);Vd(1)=1;Vd(n)=1;Va=-ones(1,n-1);Va(n-1)=0;Vc=-ones(1,n-1);Vc(1)=0;Vb(1)=c1;Vb(n)=c2;forj=2:mfori=2:n-1Vb(i)=U(i-1,j-1)+U(i+1,j-1)+s2*U(i,j-1);endX=trisys(Va,Vd,Vc,Vb);-371-

U(1:n,j)=X';endU=U';functionX=trisys(A,D,C,B)%Input-Aisthesubdiagonalofthecoefficientmatrix% -Disthemaindiagonalofthecoefficientmatrix% -Cisthesuperdiagonalofthecoefficientmatrix% -Bistheconstantvectorofthelinearsystem%Output-XisthesolutionvectorN=length(B);fork=2:Nmult=A(k-1)/D(k-1);D(k)=D(k)-mult*C(k-1);B(k)=B(k)-mult*B(k-1);endX(N)=B(N)/D(N);fork=N-1:-1:1X(k)=(B(k)-C(k)*X(k+1))/D(k);end2.3雙曲型方程的差分解法一階雙曲型方程的差分格式考慮一階雙曲型方程的初值問(wèn)題：跑U 出=0t>0,-8<X<+8R7+achc (30)名(X,0)=中(X) -8<X<+8將X-1平面剖分成矩形網(wǎng)格，取X方向步長(zhǎng)為h，t方向步長(zhǎng)為C，網(wǎng)格線為x=x=kh(k=0,±1,±2,L),t=t=戶(j=0,1,2L)為簡(jiǎn)便，記k j(k,j)=(xk,t),u(k,j)=u(xk,t),Q=我x)以不同的差商近似偏導(dǎo)數(shù)，可以得到方程的不同的差分近似u-uu-uk,j+*—kj+ak+1,jh-kj—0u—uu—uk,j+1 心+ak-ji k—1,j=0T hu—uk+1,j2hk—1,j截?cái)嗾`差分別為O(T+h)，O(T+h)與O(T+h2)。結(jié)合離散化的初始條件，可以得到幾種簡(jiǎn)單的差分格式-372-

TOC\o"1-5"\h\z鄢=u -ar(u -u)\o"CurrentDocument"；k* k,j k+i,j k,j(k=0,±1,±2,L,j=0,1,2,L) (34)&二中k,0 k鄢=u -ar(u -u )；kj+1 Q k,j k-1,j (k=0,±1,±2,L,j=0,1,2,L) (35)二Wk,0 k球=uar(u -u) , ,.k,j+1 k,j-2 k+1,j k-1,j(k=0,±1,±2,L,j=0,1,2L) (36)領(lǐng)=中k,0 kT其中：r=-。如果已知第j層節(jié)點(diǎn)上的值u，按上面三種格式就可求出第j+1層h k,j的值u 。因此，這三種格式都是顯式格式。k,j+1au auTOC\o"1-5"\h\z如果對(duì)外采用向后差商，絲采用向前差商，則方程可化成at axu（k，/）二u（k，/-1） u（k+Lj）小u（k，/）+0（^+hx_o （37）T+ah+ua十）一相應(yīng)的差分格式為：和=u-ar（u -u）.k,j+1 k,j k+1,j+1 k,j+1 （k=0,±1,±2,L,j=0,1,2,L） （38）▼=Wk,0 k此差分格式是一種隱式格式，必須通過(guò)解方程組才能由第j層節(jié)點(diǎn)上的值u，求出第k,jj+1層節(jié)點(diǎn)上的值ukj+1。t>0,-8<X<+8磔t>0,-8<X<+8例3對(duì)初值問(wèn)題：甲+書(shū)X▼(x,0)=我x)4x<0其中中(其中中(x)=".2x=0，用差分格式求其數(shù)值解u(j=1,2,3,4)，取r=-=。k,j h2x>0解：記xrkh(k=0,±1,±2,L)，由初始條件:4k=-1,-2,L女①=我x)=k=0TOC\o"1-5"\h\zkk .2▼0k=1,2,L按差分格式：4u =u-ar(u -u).k,j+1 k,j k+1,j k,j(k=0,±1,±2,L,j=0,1,2L)▼=Wk,0 k-373-

，…… 3 1 ，……計(jì)算公式為：u =-u--u 。計(jì)算結(jié)果略。k,j+12k，j2k+1，j如果用差分格式：.=u,一ar（u—u,一） , ,.k,j+1 k,j k，j kTj （k=0,±1,±2,L,j=0,1,2,L）二W▼k,0k求解，計(jì)算公式為：-1^u—（u+u）k,j+12k,jk-1,j計(jì)算結(jié)果略。與準(zhǔn)確解：4x<t■u（x,t）——x—t.2，x>0比較知，按前一個(gè)差分格式所求得的數(shù)值解不收斂到初值問(wèn)題的解，而后一個(gè)差分格式的解收斂到原問(wèn)題的解。2.3.2波動(dòng)方程的差分格式對(duì)二階波動(dòng)方程（10）TOC\o"1-5"\h\z\o"CurrentDocument"d2u d2u—a2St2 d.x2\o"CurrentDocument"du Su如果令2=萬(wàn)"2—不,則方程（10）可化成一階線性雙曲型方程組（39）記V—（v,v）T記V—（v,v）T，則方程組（39）可表成矩陣形式12（40）加Ra2/SV?（40）況4 0戶x Adx矩陣A有兩個(gè)不同的特征值入-土a，故存在非奇異矩陣P，使得PAP-1方程組（31）可化成作變換0—PV—(W,W)T，方程組（31）可化成12（41）方程組（41）由兩個(gè)獨(dú)立的一階雙曲型方程聯(lián)立而成。因此我們可以把二階波動(dòng)方程化成一階方程組進(jìn)行求解。下面給出如下波動(dòng)方程和邊界條件的差分格式-374-

u(x,t)=c2u(x,y)tt xx*(0,t)=0,u(a,t)=0X(x,0)=f(t) 0<x<a*(x,0)=g(x) 0<x<at將矩形R={(x,t):0<x<a,0<t<b}劃分成(n—1)x(m-1)個(gè)小矩形，別為：軟=h,?=k，形成一個(gè)網(wǎng)格。把方程(42)離散化成差分方程u—2u+uu—2u+ui,j+1 i,ji,j-1 k2 =c2i+1,jh夕jTj長(zhǎng)寬分為方便起見(jiàn)，可將r=ck/h代入上式，可得：長(zhǎng)寬分u —2u+u =r2(u—2u+u) (45)i,j+1 i,j i,j—1 i+1,j i,j i—1,j設(shè)行j和j—1的近似值已知，可用上式求網(wǎng)格的行j+1u =(2—2r2)u+r2(u「.+uQ-u「 (46)i,j+1 i,j i+1,ji—1,j i,j-1用上式時(shí)，必須注意，如果計(jì)算的某個(gè)階段帶來(lái)的誤差最終會(huì)越來(lái)越小，則方法是穩(wěn)定的。為了保證上式的穩(wěn)定性，必然使r=ck/h<1。還存在其他一些差分方程方法，稱為隱格式法，它們更難實(shí)現(xiàn)，但對(duì)r無(wú)限制。2.3.3差分方法求解波動(dòng)方程的MATLAB程序求解區(qū)間R={(x,t):0<x<a,0<t<b}，以(43)為初邊值條件的波動(dòng)方程的差分方法程序。%**********************************************************functionU=finedif(f,g,a,b,c,n,m)%Input-f=u(x,0)asastring'f'% - g=ut(x,0)asa string'g'% - a and brightendpointsof [0,a]and [0,b]% - c the constant inthewave equation% - n and mnumber ofgridpointsover[0,a] and[0,b]%Output-Usolutionmatrix;%IffandgareM-filefunctionscallU=finedif(@f,@g,a,b,c,n,m).%iffandgareanonymousfunctionscallU=finedif(f,g,a,b,c,n,m).%InitializeparametersandUh=a/(n-1);k=b/(m-1);r=c*k/h;r2=r△2;r22=r△2/2;si=1-r人2;s2=2-2*r2U=zeros(n,m);%Computfirstandsecondrowsfori=2:n-1U(i,1)=feval(f,h*(i-1));U(i,2)=s1*feval(f,h*(i-1))+k*feval(g,h*(i-1))...

+r22*(feval(f,h*i)+feval(f,h*(i-2)));end-375-%ComputeremainingrowsofUforj=3:m,fori=2:(n-1),U(i,j)=s2*U(i,j-1)+r2*(U(i-1,j-1)+U(i+1,j-1))-U(i,j-2);endendU=U';Cvt* K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*K,t*2個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)個(gè)一維狀態(tài)空間的偏微分方程的MATLAB解法Matlab有專門的指令pdepe解一維的拋物型或橢圓型方程組的初邊值問(wèn)題，使用指令pdepe要寫(xiě)三個(gè)輔助函數(shù)文件，即微分方程函數(shù)文件，初始條件函數(shù)文件，邊界條件函數(shù)文件。pdepe的用法Matlab提供了一個(gè)指令pdepe，用以解以下的一維偏微分方程式TOC\o"1-5"\h\z/,包、包 s 包 包、c(x,t,u, ) -(xmf(x,t,u, ))+s(x,t,u, ) (47)SxSt-x-mSx Sx Sx其中u可以為向量，時(shí)間介于10<t<t之間，而位置x則介于有限區(qū)域［a,b］之間。m值表示問(wèn)題的對(duì)稱性，其可為01或2分別表示平板(slab)，圓柱(cylindrical)或球體(spherical)的情形。式中f(x,t,u,SuSx)為通量項(xiàng)(flux)，而s(x,t,u,證8)為源項(xiàng)(source)oc(x,t,u,Su；Sx)為偏微分方程的系數(shù)對(duì)角矩陣。若某一對(duì)角線元素為0,則表示該偏微分方程為橢圓型偏微分方程，若為正值(不為0),則為拋物型偏微分方程。請(qǐng)注意c的對(duì)角線元素一定不全為0。偏微分方程初始值可表示為TOC\o"1-5"\h\zu(x,t)=v(x) (48)而邊界條件為 0 0p(x,t,u)+q(x,t)f(x,t,u,Su/Sx)=0,x=a,t0<t<t (49)pb(x,t,u)+qb(x,t)f(x,t,u,SuSx)-0,x=b,t0<t<t解上述初邊 (50)值條件的偏微分方程的MATLAB命令pdepe的用法如下：sol=pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)其中：m為問(wèn)題之對(duì)稱參數(shù)；xmesh為空間變量x的網(wǎng)格點(diǎn)(mesh)位置向量，即xmesh=［x0,x/L,xN］，其中x0=a(起點(diǎn))，xN-b(終點(diǎn))。tspan為時(shí)間變量t的向量，即tspan=［t0,t,L,tM］，其中t°為起始時(shí)間，tM=t為終點(diǎn)時(shí)間。mpdefun為使用者提供的pde函數(shù)文件。其函數(shù)格式如下：lc,f,s］=pdefun(x,t,u,dudx)亦即，使用者僅需提供偏微分方程中的系數(shù)向量。c,f和s均為歹U(column)向量，而向量c即為矩陣c的對(duì)角線元素。icfun提供解u的起始值，其格式為u-icfun(x)，值得注意的是u為列向量。-376-bcfun為使用者提供的邊界條件函數(shù)，格式如下：(pa,qa,pb,qb]=bcfun(xa,ua,xb,ub,t)其中ua和ub分別表示左邊界(xa=a)和右邊界(xb=b)u的近似解。輸出變量中，pa和qa分別表示左邊界pQ和qb的列向量，而pb和qb則為右邊界pb和qb的列向量。 a "sol為輸出的解，sol為三維矩陣，sol(j,k,i)是u的第i個(gè)分量在t=tspan(j)，x=xmesh(k)的值。options為求解器的相關(guān)解法參數(shù)。詳細(xì)說(shuō)明參見(jiàn)odeset的使用方法。注：MATLABPDE求解器pdepe的算法，主要是將原來(lái)的橢圓型和拋物型偏微分方程轉(zhuǎn)化為一組常微分方程。此轉(zhuǎn)換的過(guò)程是基于使用者所指定的mesh點(diǎn)，以二階空間離散化(spatialdiscretization技術(shù)為之(KeelandBerzins,1990)然后以ode15s的指令求解。采用ode15s的ode解法，主要是因?yàn)樵陔x散化的過(guò)程中，橢圓型偏微分方程被轉(zhuǎn)化為一組代數(shù)方程，而拋物型的偏微分方程則被轉(zhuǎn)化為一組聯(lián)立的微分方程。因而，原偏微分方程被離散化后，變成一組同時(shí)伴有微分方程與代數(shù)方程的微分代數(shù)方程組，故以ode15s便可順利求解。x的取點(diǎn)(mesh)位置對(duì)解的精確度影響很大，若pdepe求解器給出“…h(huán)asdifficultyfindingconsistentinitialconsidition”的訊息時(shí)，使用者可進(jìn)一步將mesh點(diǎn)取密一點(diǎn)，即增加mesh點(diǎn)數(shù)。另外，若狀態(tài)u在某些特定點(diǎn)上有較快速的變動(dòng)時(shí)，亦需將此處的點(diǎn)取密集些，以增加精確度。值得注意的是pdepe并不會(huì)自動(dòng)做xmesh的自動(dòng)取點(diǎn)，使用者必須觀察解的特性，自行作取點(diǎn)的操作。一般而言，所取的點(diǎn)數(shù)至少需大于3以上。tspan的選取主要是基于使用者對(duì)那些特定時(shí)間的狀態(tài)有興趣而選定。而間距(stepsize)的控制由程序自動(dòng)完成。若要獲得特定位置及時(shí)間下的解，可配合以pdeval命令。使用格式如下：luout,duoutdx一=pdeval(m,xmesh,ui,xout)其中m代表問(wèn)題的對(duì)稱性。m=0表示平板；m=1表示圓柱體；m=2表示球體。其意義同pdepe中的自變量m。xmesh為使用者在pdepe中所指定的輸出點(diǎn)位置向量。xmesh=[x0,x「L,xN]oui即sol(j,:,i)o也就是說(shuō)其為pdepe輸出中第i個(gè)輸出ui在各點(diǎn)位置】xmesh處，時(shí)間固定為t=tspan(j)下的解。jxout為所欲內(nèi)插輸出點(diǎn)位置向量。此為使用者重新指定的位置向量。uout為基于所指定位置xout，固定時(shí)間tf下的相對(duì)應(yīng)輸出。duoutdx為相對(duì)應(yīng)的dudx輸出值。以下將以數(shù)個(gè)例子，詳細(xì)說(shuō)明pdepe的用法。/2求解一維偏微分方程試解偏微分方程2du du82而二布2其中0<x<1,且滿足以下初邊值條件-377-

⑴初值條件：u(X,0)=sin(兀x)(ii)邊界條件BC1：u(0,t)=0；BC2：兀et/u(17)=0。+ax注：本問(wèn)題的解析解為u(x,t)=e-tsin(Kx)解下面將敘述求解的步驟與過(guò)程。當(dāng)完成以下各步驟后，可進(jìn)一步將其匯總為主程序ex20_1.m，然后求解。步驟1將欲求解的偏微分方程改寫(xiě)成如下的標(biāo)準(zhǔn)形式。au aoau口八冗2—o—OF0

at_axoaxo此即au)=一ax=0和m=0。步驟2編寫(xiě)偏微分方程的系數(shù)向量函數(shù)。function[c,f,s]=ex20_1pdefun(x,t,u,dudx)c=pi人2;f=dudx;s=0;步驟3編寫(xiě)起始值條件。functionu0=ex20_1ic(x)u0=sin(pi*x);步驟4編寫(xiě)邊界條件。在編寫(xiě)之前，先將邊界條件改寫(xiě)成標(biāo)準(zhǔn)形式，如式(49)，函數(shù)，即和因而，(50)，找出相對(duì)應(yīng)的函數(shù)pa(?),qa(.)和pb(?),qb。，然后寫(xiě)出函數(shù)，即和因而，aBC1：u(0,t)+0-a(0,t)=0,x=0

axauBC2：兀e-1+1? (1t)=0,x=1axp=u(0,t),q=0p=Ke-1,q=1邊界條件函數(shù)可編寫(xiě)成function[pa,qa,pb,qb]=ex20_1bc(xa,ua,xb,ub,t)pa=ua;qa=0;pb=pi*exp(-t);qb=1;步驟5取點(diǎn)。例如x=linspace(0,1,20);%x取20點(diǎn)t=linspace(0,2,5);%時(shí)間取5點(diǎn)輸出步驟6利用pdepe求解。m=0;%依步驟1之結(jié)果-378-sol=pdepe(m,@ex20_1pdefun,@ex20_1ic,@ex20_1bc,x,t);%這里sol實(shí)際上是二維矩陣 — —步驟7顯示結(jié)果。surf(x,t,sol)title('pde數(shù)值解')，xlabel('位置')，ylabel('時(shí)間')，zlabel('u')若要顯示特定點(diǎn)上的解，可進(jìn)一步指定x或t的位置，以便繪圖。例如，欲了解時(shí)間為2(終點(diǎn))時(shí)，各位置下的解，可輸入以下指令(利用pdeval指令)：figure(2);%繪成圖2M=length(t);%取終點(diǎn)時(shí)間的下標(biāo)xout=linspace(0,1,100);%輸出點(diǎn)位置[uout,dudx]=pdeval(m,x,sol(M,:),xout);plot(xout,uout);%繪圖title('時(shí)間為2時(shí)，各位置下的解')，xlabel('x'),ylabel('u')綜合以上各步驟，可寫(xiě)成一個(gè)程序求解例4。其程序(全部程序放在一個(gè)文件中)如下：functionsol=ex20_1%************************************%求解一維熱傳導(dǎo)偏微分方程的一個(gè)綜合函數(shù)程序%************************************m=0;x=linspace(0,1,20);%xmesht=linspace(0,2,20);%tspan%************%以pdepe求解軍%************sol=pdepe(m,@ex20_1pdefun,@ex20_1ic,@ex20_1bc,x,t);%這里sol實(shí)際上是二維矩陣 — —%************%繪圖輸出%************figure(1),surf(x,t,sol)title('pde數(shù)值解')，xlabel('位置x'),ylabel('時(shí)間t'),zlabel('數(shù)值解u')%*************%與解析解做比較%*************figure(2)surf(x,t,exp(-t)'*sin(pi*x));title('解析解')，xlabel('位置x'),ylabel('時(shí)間t'),zlabel('數(shù)值解u')%*****************%t=tf=2時(shí)各位置之解%*****************figure(3)M=length(t);%取終點(diǎn)時(shí)間的下標(biāo)xout=linspace(0,1,100);%輸出點(diǎn)位置[uout,dudx]=pdeval(m,x,sol(M,:),xout);plot(xout,uout);%繪圖title('時(shí)間為2時(shí)，各位置下的解')，xlabel('x'),ylabel('u')2******************-379-

%偏微分方程函數(shù)%******************function[c,f,s]=ex20_1pdefun(x,t,u,dudx)c=pi人2;f=dudx;s=0;%******************%初始條件函數(shù)%******************functionu0=ex20_1ic(x)u0=sin(pi*x);%******************%邊界條件函數(shù)%******************function[pa,qa,pb,qb]=ex20_1bc(xa,ua,xb,ub,t)pa=ua;qa=0;pb=pi*exp(-t);qb=1;例5試解以下聯(lián)立的偏微分方程系統(tǒng)i2) d2uj2)=0.170—^+F(uj2)TOC\o"1-5"\h\zd.x2 1其中F(u一u)=exp(5.73(u一u))一exp(-11.46(u一u))，且0<x<1和t>0。此聯(lián)立偏微分方程系統(tǒng)滿足以下初邊值條件。 1 2⑴初值條件u(x,0)=1，u(x,0)=0(ii)邊值條件 1 2\o"CurrentDocument" 1(0,t)=0，u(0,t)=0dx 2uj1t)=1^u^(1t)=0解步驟1：改寫(xiě)偏微分方程為標(biāo)準(zhǔn)形式玉024'"/皆*"2 /xMF(u「u2匕鏟3t<2r2x0.170吃了<F(u1-u2)/運(yùn)算表示兩個(gè)向量的對(duì)應(yīng)元素相乘。因此FF(u-FF(u-u)/s=' 1 28<F(u[u2)fc=%,f=' 3ux8<f 0.170-^.w< 3xf和m=0。另外，左邊界條件（x=0處）。寫(xiě)成-380-

X0/X/

，

<

2X0/X/

，

<

2X0.024*，0.170<學(xué)X0DSX_,絲7鮮s%/X0/Pa='<2F同理，右邊界條件同理，右邊界條件（X=1處）為石.024X-1<X-1<<10*，，0.170<步驟2：編寫(xiě)偏微分方程的系數(shù)向量函數(shù)。function[c,f,s]=ex20_2pdefun(x,t,u,dudx)c=[11]';f=[0.0240.170],.*dudx;y=u(1)-u(2);F=exp(5.73*y)-exp(-11.47*y);s=[-FF]';步驟3：編寫(xiě)初始條件函數(shù)functionu0=ex20_2ic(x)u0=[10]';一步驟4：編寫(xiě)邊界條件函數(shù)function[pa,qa,pb,qb]=ex20_2bc(xa,ua,xb,ub,t)pa=[0ua(2)]';qa=[10]';pb=[ub(1)-10]';qb=[01]';步驟5：取點(diǎn)。由于此問(wèn)題的端點(diǎn)均受邊界條件的限制，且時(shí)間t很小時(shí)狀態(tài)的變動(dòng)很大(由多次求解后的經(jīng)驗(yàn)得知），如，x=[00.005解后的經(jīng)驗(yàn)得知），如，x=[00.005t=[00.0050.010.050.10.20.50.70.90.950.990.9951];0.010.050.10.511.52];以上幾個(gè)主要步驟編寫(xiě)完成后，事實(shí)上就可直接完成主程序來(lái)求解。此問(wèn)題的求解程序如下：functionsol=ex20_2^**************************************^*%求解一維偏微分方程組的一個(gè)綜合函數(shù)程序^*^**************************************^*m=0;x=[00.0050.010.050.10.20.50.70.90.950.990.9951];t=[00.0050.010.050.10.511.52];^*^************************************^*%利用pdepe求解^*^************************************^*-381-

sol=pdepe(m,@ex20_2pdefun,@ex20_2ic,@ex20_2bc,x,t);u1=sol(:,:,1);%第一個(gè)狀態(tài)之?dāng)?shù)值解輸出u2=sol(:,:,2);%第二個(gè)狀態(tài)之?dāng)?shù)值解輸出^************************************^*%繪圖輸出^*^************************************^*figure(1),surf(x,t,u1)title('u1之?dāng)?shù)值解'),xlabel('x'),ylabel('t')%figure(2),surf(x,t,u2)title('u2之?dāng)?shù)值解'),xlabel('x'),ylabel('t')^*^**************************************^*%pdefun函數(shù)^*^**************************************^*function[c,f,s]=ex20_2pdefun(x,t,u,dudx)c=[11],;f=[0.0240.170]'.*dudx;y=u(1)-u(2);F=exp(5.73*y)-exp(-11.47*y);s=[-FF]';^*^***************************************^*%初始條件函數(shù)^*^***************************************^*functionu0=ex20_2ic(x)u0=[10]';^*^***************************************^*%邊界條件函數(shù)^*^***************************************^*function[pa,qa,pb,qb]=ex20_2bc(xa,ua,xb,ub,t)pa=[0ua(2)],;qa=[10]';pb=[ub(1)-10],;qb=[01]';3.3應(yīng)用實(shí)例9.611=0T9.611=0T二卜DT CD2T 」DT□二卜——+0.59^—+ DL □Dr2 rDrCD2fCD2f 1+0.94^-+□Dr20.047=0T其中T為溫度(℃其中T為溫度(℃),f為反應(yīng)率，L為軸向距離，r為徑向距離。此系統(tǒng)的邊界條件為r=0：T=0,f=0,寥=f=0Dr DrDTr=2：—0.65一=112(T—273)，

Dr當(dāng)L=0時(shí)，T(r,0)=125，f(r,0)=0。將原方程改寫(xiě)成如式(47)的標(biāo)準(zhǔn)式生二0Dr-382-因此?%

c=4/和m=+因此?%

c=4/和m=+1(圓柱)。/L-=r.i嗎巧下 包/,0.59r\o"CurrentDocument"， 8: 乘八1初611&094rf.T<0.047/\o"CurrentDocument"4 aL8丫 at/~0.59瓦8?1B61Uf= 8s=， 8094f8T<0.047f< al8另外，左邊界條件(r=0處)寫(xiě)成?17D◎Pa=&f=<18f同理右邊界條件(r=2)可寫(xiě)成112(T—273yD.65/0.59/112(T—273y D.65/0.59/p=8 8,q=8 8b< 0 8與< 1 /根據(jù)以上的分析，可編寫(xiě)MATLAB程序求解此PDE問(wèn)題，其求解程序如下:functionsol=main20_6m=1;%********************%取點(diǎn)%********************r=linspace(0,2,20);L=linspace(0,5,40);%***********************%利用pdepe求解%***********************sol=pdepe(m,@ex20_6pdefun,@ex20_6ic,@ex20_6bc,r,L);T=sol(:,:,1);%溫度f(wàn)=sol(:,:,2);%反應(yīng)率%***********************%繪圖輸出%***********************figure(1),surf(L,r,T')title('temp'),xlabel('L'),ylabel('r'),zlabel('temp(0C)’)%figure(2),surf(L,r,f')title('reactionrate'),xlabel('L'),ylabel('r'),zlabel('reactionrate')-383-

%*************************************************%偏微分方程函數(shù)%*************************************************function[c,f,s]=ex20_6pdefun(r,L,u,DuDr)T=u(1);f=u(2);c=[11]';f=[0.590.94]'.*DuDr;s=[9.611;0.047]/T;%**********************************%初始條件函數(shù)%**********************************functionu0=ex20_6ic(x)u0=[1250]';一%**********************************%邊界條件函數(shù)%**********************************function[pa,qa,pb,qb]=ex20_6bc(ra,ua,rb,ub,L)pa=[00]';qa=[11]';pb=[112*(ub(1)-273)0]';qb=[0.65/0.591]';例7擴(kuò)散系統(tǒng)之濃度分布管中儲(chǔ)放靜止液體B，高度為L(zhǎng)=10cm,放置于充滿A氣體的環(huán)境中。假設(shè)與B液體接觸面之濃度為CA0=0.01mol!m3，且此濃度不隨時(shí)間改變而改變，即在操作時(shí)間內(nèi)（h=10天）維持定值。氣體A在液體B中之?dāng)U散系數(shù)為D=2x10-9m2/s。試決定以下兩種情況下，氣體A溶于液體B中的流通量（flux）。ABA與B不發(fā)生反應(yīng)；A與B發(fā)生反應(yīng)解（a）因氣體A與液體B不發(fā)生反應(yīng)，故其擴(kuò)散現(xiàn)象的質(zhì)量平衡方程如下：ac s2cfa=dabaz2依題意，其初始及邊界條件為I.C. CA(z,0)=0,z>0一八 一 八 ac 八 八B.C.C(0,t)=C,t>0；號(hào)=0,t>0A A0 azz=Lac在獲得濃度分布后，即可用Fick定律acabazNabazz=z=0計(jì)算流通量。與標(biāo)準(zhǔn)式（47）比較，可得C=1,f=Da『Ca2z,s=0，和m=0。另外，經(jīng)與式（49）和（50）比較后得知，左邊界及右邊界條件的系數(shù)分別為左邊界(z=0)：p=CA(0,t)-CA0,q=0。TOC\o"1-5"\h\z\o"CurrentDocument""" 10 "\o"CurrentDocument"右邊界(z=L)：p=0,q= 。b bDAB-384-

（b）在氣體A與液體B會(huì)發(fā)生一次反應(yīng)的情況下，其質(zhì)量平衡方程需改寫(xiě)為其中k=2x10-7。而起始及邊界條件同上。與標(biāo)準(zhǔn)式（47）比較，可得m=0，C=1，f=DABaCA/az，和s=kCA。而邊界條件的系數(shù)同（a）。 的A A利用以上的處理結(jié)果，可編寫(xiě)MATLAB程序如下：functionex207%****************文************%****************文************%擴(kuò)散系統(tǒng)之濃度分布%*****************************clear,clcglobalDABkCA0%******************************%給定數(shù)據(jù)2******************************CA0=0.01;L=0.1;DAB=2e-9;k=2e-7;h=10*24*3600;%*******************************%取點(diǎn)2*******************************t=linspace(0,h,100);z=linspace(0,L,10);%*******************************%case(a)%*******************************m=0;7bc,z,t);sol=pdepe(m,@ex20_7pdefuna,@ex20_7ic,@ex20CA=sol;7bc,z,t);fori=1:length(t)[CA_i,dCAdz_i]=pdeval(m,z,CA(i,:),0);NAz(i)=-dCAdzi*DAB;endfigure(1),subplot(211)surf(z,t/(24*3600),CA),title('case(a)')xlabel('length(m)'),ylabel('time(day)'),zlabel('conc.(mol/m人3)')subplot(212)plot(t/(24*3600),NAz'*24*3600)xlabel('time(day)'),ylabel('flux(mol/m人2.day)')%************************************%c