【正文】 3. 求熱傳導(dǎo)方程初邊值問題 的數(shù)值解，并將計(jì)算結(jié)果與準(zhǔn)確解比較。 習(xí) 題 1. 求二維拉普拉斯方程在邊界條件下的數(shù)值解。 4）設(shè)置“plot菜單下的“parameters…”選項(xiàng)如下：選擇 Height(3D plot) 和 Animation兩項(xiàng)。 其中的變量 x，y 是 MATLAB 可接受的內(nèi)置變量。 ix=find(x.^2+y.^2)。 3）然后編寫函數(shù) fun1(x,y)如下： function f=fun1(x,y)。我們這里就 省略了。 例 11 使用 PDETOOL 重新求例 8 的數(shù)值解。 （1）在邊界條件輸入框中，可以使用如下變量： 二維坐標(biāo) x 和 y，邊界線段長度參數(shù) s（s 是以箭頭的方向沿邊界線段從 0 增加到 1），外法向矢量的分量 nx 和 ny（如果需要邊界的切線方向，可以通過 tx=ny 和 ty=nx 表示），解 u。 3. 如 果 求拋 物 型 或 雙曲 型 方 程 的數(shù) 值 解 ， 還需 要 通 過 “solve ” 菜 單下 的 “parameters…”選項(xiàng)設(shè)置初值條件。參數(shù)設(shè)置見圖 10，顯示結(jié)果見圖 11。 注意： 1. MATLAB 會(huì)以圖形的方式展示結(jié)果，使用者亦可點(diǎn)選 plot 下之“parameters”功能，選擇適當(dāng)?shù)姆绞斤@示圖形及數(shù)據(jù)。使用者亦可進(jìn)一步利用將網(wǎng)點(diǎn)取得密一 點(diǎn)(refine mesh)，見圖 8。對話框輸入情況見圖 7。其中在弧形部份與標(biāo)準(zhǔn)式知，g=5 且q=0。 步驟 4：選取 Boundary 功能，以輸入邊界條件。 步驟 3：選取 PDE 功能項(xiàng)，以輸入 PDE 方程的系數(shù)及類型。同時(shí)，如欲進(jìn)一步修改各圖形對象之大小及位置數(shù)據(jù)，可在該圖上雙擊鼠標(biāo)左鍵，然后在對象對話框上輸入數(shù)據(jù)。其余圖形 C1，R1 和 C2 可選取適當(dāng)對象并類似地畫出，以形成如圖 4 的圖形區(qū)域。 步驟 2：利用 Draw 功能，畫出問題之幾何圖形。 步驟 1：在命令窗口中鍵入 pdetool 以進(jìn)入 GUI(graphical user interface)界面。 s方程式的求解為例，詳細(xì)說明 pdetool 的用法。 （3）最后，在 Plot 模式下，顯示答案。 在定義 PDE 問題之后，可依以下兩個(gè)步驟求解 （1）在 mesh 模式下，產(chǎn)生 mesh 點(diǎn)，以便將原問題離散化。 （2）利用 boundary 模式，指定邊界條件。 ：放大縮小功能，便于圖形繪制及顯示。 ：在指定 PDE 系統(tǒng)，邊界條件及區(qū)域后，按此鈕即開始解題。 ：產(chǎn)生圖形區(qū)域內(nèi)離散化的網(wǎng)點(diǎn)。在此功能選定后，使用者可在任一圖形邊界上按住鼠 標(biāo)左鍵雙擊，然后在對話框中輸入邊界條件。 ：用以繪制多邊型等不規(guī)則區(qū)域，欲關(guān)閉此功能需按鼠標(biāo)右鍵。 ：以中心點(diǎn)向外的方式繪制橢圓或圓。 ：由周圍界線的方式繪制橢圓或圓形區(qū)域。 ：從中心點(diǎn)至某一角邊的方式繪制矩形或正方形。現(xiàn)將這些按鈕的主要功能敘述如下： 前五個(gè)按鈕為 PDE 系統(tǒng)之邊界范圍繪制功能，由左至右之用法為： ：以對角繪制矩形或正方形。 在說明此解法工具之前，先介紹此 PDE 圖形界面的菜單下方的功能圖標(biāo)(icon)按鈕。 （2）產(chǎn)生離散化之點(diǎn)，并將原 PDE 方程式離散化。) 偏微分方程的 pdetool 解法 圖形界面解法簡介 對于一般的區(qū)域，任意邊界條件的偏微分方程，我們可以利用 MATLAB 中 pdetool提供的偏微分方程用戶圖形界面解法。 subplot(2,2,1),pdemesh(p,e,t)。 u=assempde(b,p,e,t,c,a,f)。 [p,e,t]=initmesh(g)。f=39。 c=1。 b=39。 解 它的精確解為 下面求它的數(shù)值解，編寫程序如下： g=39。 mv(j)=getframe。 uu=hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d)。 %(4)在時(shí)間段為0到5的31個(gè)點(diǎn)上求解 n=31。 u0=atan(cos(pi*x))。 %(3)定義初始條件 x=p(1,:)39。d=1。a=0。squareb339。squareg39。1上滿足Neumann條件 = 0 解 這里是雙曲型方程，其中c = 1, a = 0, f = 0, d = 1。 end movie(mv,10) 例9 求解正方形區(qū)域{(x, y) |≤ x, y ≤ 1}上的波方程 初始條件為u(0) = arctan(cos(πx))， 3sin(πx) exp(cos(πy ))，邊界條件為x =177。 %(5)動(dòng)畫圖示結(jié)果 for j=1:nframe pdesurf(p,t,u1(:,j))。 tlist=linspace(0,nframe)。 ix=find(sqrt(p(1,:).^2+p(2,:).^2))。 %（2）產(chǎn)生初始的三角形網(wǎng)格 [p,e,t]=initmesh(g)。f=0。 %邊界上為零條件 c=1。 %定義正方形區(qū)域 b=39。編寫程序如下： %(1)問題定義 g=39。 pdesurf(p,t,u) 例8 求解正方形區(qū)域{(x, y) |≤ x, y ≤ 1}上的熱傳導(dǎo)方程初始條件為 邊界條件為Dirichlet條件u = 0。Tol39。 [p,e,t]=refinemesh(g,p,e,t)。 rtol=1e3。 c=39。 b=39。 解 這是橢圓型方程，其中，編寫程序如下：g=39。 subplot(2,2,2),pdesurf(p,t,u) subplot(2,2,3),pdesurf(p,t,uexact39。 error=[error err]。 err=norm(uexact39。 u=assempde(b,p,e,t,c,a,f)。 err=1。 %（2）產(chǎn)生初始的三角形網(wǎng)格 [p,e,t]=initmesh(g)。a=0。circleb139。circleg39。 求解偏微分方程 例 6 求解泊松方程 ，求解區(qū)域?yàn)閱挝粓A盤，邊界條件為在圓盤邊界上u = 0。 （ii）Neumann 條件： () + qu = g on這里為區(qū)域的單位外法線，h, r, q, g是定義在上的復(fù)值函數(shù)。 （v）非線性橢圓偏微分方程 其中c, a, f 可以是u的函數(shù)。 （ii）拋物型偏微分方程 其中c,a, f ,d 可以依賴于時(shí)間t。 下面我們討論的方程是定義在平面上的有界區(qū)域Ω上，區(qū)域的邊界記作。 167。 pr=0。 %************************************************ % 邊界條件函數(shù) %************************************************ function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t) global DAB k CA0 pl=CAlCA0。 s=k*CA。 %********************************************* % case (a) %********************************************* function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz) global DAB k CA0 c=1。 f=DAB*dCAdz。flux (mol/m^)39。time (day)39。) subplot(212) plot(t/(24*3600),NAz39。) zlabel(39。) ylabel(39。) xlabel(39。 end % figure(2) subplot(211) surf(z,t/(24*3600),CA) title(39。 for i=1:length(t) [CA_i,dCAdz_i]=pdeval(m,z,CA(i,:),0)。 sol=pdepe(m,ex20_3_2pdefunb,ex20_3_2ic,ex20_3_2bc,z,t)。flux (mol/m^)39。time (day)39。) subplot(212) plot(t/(24*3600),NAz39。) zlabel(39。) ylabel(39。) xlabel(39。 end figure(1) subplot(211) surf(z,t/(24*3600),CA) title(39。 for i=1:length(t) [CA_i,dCAdz_i]=pdeval(m,z,CA(i,:),0)。 sol=pdepe(m,ex20_3_2pdefuna,ex20_3_2ic,ex20_3_2bc,z,t)。 z=linspace(0,L,10)。 h=10*24*3600。 DAB=2e9。 利用以上的處理結(jié)果，可編寫 MATLAB 參考程序如下： function ex20_3_2 %***************************** % 擴(kuò)散系統(tǒng)之濃度分布 %***************************** clear clc global DAB k CA0 %****************************** % 給定數(shù)據(jù) %****************************** CA0=。另外，經(jīng)與式（37）比較后得知，左邊界及右邊界條件之系數(shù)分別為 左邊界( z = 0 )： 右邊界( z = L )： (b)與標(biāo)準(zhǔn)式（35）比較，可得m = 0，C = 1，和 ?，F(xiàn)就各狀況的處理過程簡述如下。 在獲得濃度分布后，即可以 Fick’s law計(jì)算流通量。 (a) A 與 B 不發(fā)生反應(yīng)； (b) A 與 B 發(fā)生以下之反應(yīng)