偏微分方程求解

1、關(guān)于偏微分方程求解第一張，PPT共四十六頁，創(chuàng)作于2022年6月基礎(chǔ)知識預(yù)習(xí)微分方程的求解包含：常微分方程的求解（上節(jié)課已經(jīng)講過）這里不再贅述。：偏微分方程的求解（本次教學(xué)內(nèi)容）第二張，PPT共四十六頁，創(chuàng)作于2022年6月偏微分方程概念偏微分方程（Partial Differential Equation，簡稱PDE）指含有未知函數(shù)及其偏導(dǎo)數(shù)的方程。描述自變量、未知函數(shù)及其偏導(dǎo)數(shù)之間的關(guān)系。偏微分方程分為線性偏微分方程式與非線性偏微分方程式，常常有幾個解而且涉及額外的邊界條件。常微分方程：在微分方程中，若自變量的個數(shù)只有一個的微分方程。偏微分方程：自變量的個數(shù)有兩個或兩個以上的微分方程。第三

2、張，PPT共四十六頁，創(chuàng)作于2022年6月求解偏微分方程的方法求解偏微分方程的數(shù)值方法：1. 有限元法（Finite Element Method, FEM）- hp-FEM2. 有限體積法（Finite Volume Method, FVM）3. 有限差分法（Finite Difference Method, FDM）。其它：廣義有限元法（Generalized Finite Element Method, FFEM）、擴(kuò)展有限元法（eXtended Finite Element Method, XFEM）、無網(wǎng)格有限元法（Meshfree Finite Element Method）、離散

3、迦遼金有限元法（Discontinuous Galerkin Finite Element Method, DGFEM）等。第四張，PPT共四十六頁，創(chuàng)作于2022年6月MATLAB解偏微分方程MATLAB提供了兩種方法解決PDE 問題：pdepe()函數(shù)，它可以求解一般的PDEs，具有較大的通用性，但只支持命令行形式調(diào)用。PDE 工具箱，可以求解特殊PDE 問題，PDEtool 有較大的局限性，比如只能求解二階PDE 問題，并且不能解決偏微分方程組，但是它提供了GUI界面，從繁雜的編程中解脫出來了，同時還可以通過File-Save As直接生成M代碼使用pdeval()直接計算某個點(diǎn)的函數(shù)值

4、？第五張，PPT共四十六頁，創(chuàng)作于2022年6月一般偏微分方程組(PDEs)的MATLAB求解直接求解一般偏微分方程(組)，它的調(diào)用格式為sol=pdepe(m,pdefun,pdeic,pdebc,x,t)問題描述函數(shù)初值條件邊界條件輸出參數(shù)自變量參數(shù)第六張，PPT共四十六頁，創(chuàng)作于2022年6月【輸入?yún)?shù)】（1）pdefun：是PDE 的問題描述函數(shù)，它必須換成下面的標(biāo)準(zhǔn)形式PDE 就可以編寫下面的入口函數(shù)c,f,s=pdefun(x,t,u,du)m,x,t就是對應(yīng)于(式1)中相關(guān)參數(shù)和自變量，du是u的一階導(dǎo)數(shù)，由給定的輸入變量即可表示出出c,f,s這三個函數(shù)第七張，PPT共四十六頁，

5、創(chuàng)作于2022年6月【輸入?yún)?shù)】（2）pdeic：是PDE 的初值條件，必須化為下面的形式我們使用下面的簡單的函數(shù)來描述為u0=pdeic(x)第八張，PPT共四十六頁，創(chuàng)作于2022年6月【輸入?yún)?shù)】（3）pdebc：是PDE的邊界條件描述函數(shù)，必須先化為下面的形式于是邊值條件可以編寫下面函數(shù)描述為pa,qa,pb,qb=pdebc(x,t,u,du)其中a 表示下邊界，b 表示下邊界第九張，PPT共四十六頁，創(chuàng)作于2022年6月【輸入?yún)?shù)】（4）m：就是對應(yīng)于(式1)中相關(guān)參數(shù)x,t：就是對應(yīng)于(式1)中自變量第十張，PPT共四十六頁，創(chuàng)作于2022年6月【輸出參數(shù)】sol：是一個三維數(shù)組

6、，sol(:,:,i)表示ui的解，換句話說uk對應(yīng)x(i)和t(j)時的解為sol(i,j,k)第十一張，PPT共四十六頁，創(chuàng)作于2022年6月實(shí)例講解（題目）例：第十二張，PPT共四十六頁，創(chuàng)作于2022年6月初值條件邊界條件第十三張，PPT共四十六頁，創(chuàng)作于2022年6月實(shí)例講解（解法）【解】第一步根據(jù)（1）對照給出的偏微分方程，則原方程可以改寫為第十四張，PPT共四十六頁，創(chuàng)作于2022年6月輸入?yún)?shù)（1）目標(biāo)PDE函數(shù)% 目標(biāo)PDE函數(shù)function c,f,s=pdefun (x,t,u,du)c=1;1;f=0.024*du(1);0.17*du(2);temp=u(1)-u(

7、2);s=-1;1.*(exp(5.73*temp)-exp(-11.46*temp);第十五張，PPT共四十六頁，創(chuàng)作于2022年6月輸入?yún)?shù)（2）初值條件初值條件改寫為% 初值條件函數(shù)function u0=pdeic(x)u0=1;0;第十六張，PPT共四十六頁，創(chuàng)作于2022年6月輸入?yún)?shù)（3）邊界條件邊界條件改寫為% 邊界條件函數(shù)function pa,qa,pb,qb=pdebc(xa,ua,xb,ub,t)%a表示左邊界，b表示右邊界pa=0;ua(2);qa=1;0;pb=ub(1)-1;0;qb=0;1;第十七張，PPT共四十六頁，創(chuàng)作于2022年6月（4）主調(diào)函數(shù)clcx=

8、0:0.05:1;t=0:0.05:2;m=0;sol=pdepe(m,pdefun,pdeic,pdebc,x,t);figure(numbertitle,off,name,PDE Demoby Matlabsky)%創(chuàng)建個窗口，窗口名字是name后邊的名字NumberTitle,off是關(guān)掉默認(rèn)顯示名字。subplot(211)surf(x,t,sol(:,:,1)%sol(:,:,i)表示ui的解title(The Solution of u_1)xlabel(X)ylabel(T)zlabel(U)subplot(212)surf(x,t,sol(:,:,2)%sol(:,:,i)表示

9、ui的解title(The Solution of u_2)xlabel(X)ylabel(T)zlabel(U)第十八張，PPT共四十六頁，創(chuàng)作于2022年6月第十九張，PPT共四十六頁，創(chuàng)作于2022年6月PDEtool求解特殊PDE問題MATLAB的偏微分工具箱(PDE toolbox)可以比較規(guī)范的求解各種常見的二階偏微分方程(特殊二階的PDE)第二十張，PPT共四十六頁，創(chuàng)作于2022年6月典型偏微分方程的描述 第二十一張，PPT共四十六頁，創(chuàng)作于2022年6月 第二十二張，PPT共四十六頁，創(chuàng)作于2022年6月 （3）雙曲線型偏微分方程的一般形式 第二十三張，PPT共四十六頁，創(chuàng)作

10、于2022年6月 （4）特征值型偏微分方程的一般形式，注 意它是（1）的變形，不能算獨(dú)立的一類 第二十四張，PPT共四十六頁，創(chuàng)作于2022年6月MATLAB 采用有限元的方法求解各種PDEMATLAB 為我們提供一個pdetool (在command window 中鍵輸pdetool打開)的交互界面，可以求解二元偏微分u(x1,x2)(注意只能求解二元)。方程的參數(shù)由a、c、d和f確定，求解域由圖形確定，求解域確定好后，需要對求解域進(jìn)行柵格化(這個是自動)。第二十五張，PPT共四十六頁，創(chuàng)作于2022年6月 第二十六張，PPT共四十六頁，創(chuàng)作于2022年6月 第二十七張，PPT共四十六頁，

11、創(chuàng)作于2022年6月偏微分方程邊界條件的描述Dirichlet(狄利克萊)條件Neumann(紐曼)條件第二十八張，PPT共四十六頁，創(chuàng)作于2022年6月 第二十九張，PPT共四十六頁，創(chuàng)作于2022年6月 第三十張，PPT共四十六頁，創(chuàng)作于2022年6月求解實(shí)例 第三十一張，PPT共四十六頁，創(chuàng)作于2022年6月【解】由給定的PDE，可以得出d=1,c=1,a=2,f=10第三十二張，PPT共四十六頁，創(chuàng)作于2022年6月step1:點(diǎn)擊工具欄的【PDE】按鈕，如下輸入PDE的參數(shù)，注意選擇Hyperbolic 第三十三張，PPT共四十六頁，創(chuàng)作于2022年6月step2:繪制求解域?qū)ψ鴺?biāo)軸

12、的操作可以在【Options】主菜單中操作，包括設(shè)置網(wǎng)格、坐標(biāo)系范圍等(1)【Options】-Axis Limits設(shè)置如下第三十四張，PPT共四十六頁，創(chuàng)作于2022年6月 第三十五張，PPT共四十六頁，創(chuàng)作于2022年6月 (2)點(diǎn)擊工具欄上的第三個按鈕【繪制橢圓】，任意繪制一個橢圓，雙擊橢圓，設(shè)置如下 第三十六張，PPT共四十六頁，創(chuàng)作于2022年6月重復(fù)上面的操作，參數(shù)如下 第三十七張，PPT共四十六頁，創(chuàng)作于2022年6月得到 第三十八張，PPT共四十六頁，創(chuàng)作于2022年6月 (3)在set formula 中如下輸入，“+”表示求并集，“-”表示求差集，注意沒有直接求交接的操作

13、符第三十九張，PPT共四十六頁，創(chuàng)作于2022年6月step3:邊界條件和初值條件初值條件可以通過【Solve】-【Parameters】設(shè)置邊值條件設(shè)置如下(1)點(diǎn)擊工具欄的第6 個按鈕【區(qū)域邊界】，顯示如下第四十張，PPT共四十六頁，創(chuàng)作于2022年6月 (2)【Boundary】-【Remove All Subdomain Borders】移除所有子域的邊界，將得到所有子域合并成一個求解域(3) 【Boundary】-【Secify Boundary Conditons】設(shè)置邊界如下，注意我們這里只有Dirichlet條件第四十一張，PPT共四十六頁，創(chuàng)作于2022年6月step4:生成

14、使用有限元方法求解方程所需的柵格點(diǎn)擊工具欄的第8/9 個按鈕，對求解域生成柵格，多次點(diǎn)擊可以在原來基礎(chǔ)上繼續(xù)細(xì)化柵格，直到自己覺得滿意為止，當(dāng)然可以通過【Mesh】主菜單進(jìn)行精確控制第四十二張，PPT共四十六頁，創(chuàng)作于2022年6月 step5:求解方程點(diǎn)解工具欄的第10 個按鈕“=”【求解方程】 step6:求解結(jié)果繪圖 點(diǎn)擊第11 個按鈕【繪制圖形】，里面的選項很豐富，可以繪制等高線等好多，甚至播放動畫，具體大家可以自己慢慢摸索第四十三張，PPT共四十六頁，創(chuàng)作于2022年6月動畫播放設(shè)置：(1)【Solve】-【Parameters】設(shè)置合適的時間向量Time(2)【Plot】-【Parameters】選中【Animation】，點(diǎn)擊后面的【Options】，設(shè)置播放速度和次數(shù)，比如6fp