Mathematica 12 為偏微分方程（PDE）的符號和數(shù)值求解提供了強大的功能。本文將重點介紹版本12中全新推出的基于有限元方法（FEM）的非線性PDE求解器。首先簡要回顧用于求解 PDE 的 Wolfram 語言基本語法，包括如何指定狄利克雷和諾伊曼邊界條件；隨后我們將通過一個具體的非線性問題，說明 Mathematica 12的 FEM 求解過程。最后，我們將展示一些物理和化學(xué)實例，如Gray-Scott模型和與時間相關(guān)的納維-斯托克斯方程。更多信息可以在 Wolfram 語言教程"有限元編程"中找到，本文大部分內(nèi)容都以此為基礎(chǔ)（教程鏈接見文末）。

 

1. 前言

Wolfram Research 的旗艦產(chǎn)品 Mathematica 通過5000多種內(nèi)置函數(shù)驅(qū)動著 Wolfram 語言。在作為數(shù)學(xué)建模和分析基礎(chǔ)的常/偏微分方程領(lǐng)域，Mathematica 12 具有功能強大的求解器來對其進行符號或數(shù)值求解。最近，基于有限元法的數(shù)值求解函數(shù)得到顯著增強，并有望求解任意區(qū)域上的PDE并獲得特征值/特征函數(shù)。在此，我們將著重介紹 FEM 在最新版本12中對非線性偏微分方程的求解，并通過實例介紹在實際問題中的應(yīng)用流程。使用有限元方法求解非線性 PDE 的詳細過程和代碼信息向公眾開放，請參見Wolfram 語言教程"有限元編程"。

 

2. 微分方程的數(shù)值求解過程

 

在 Wolfram 語言中，對微分方程進行數(shù)值求解的函數(shù)有兩個：NDSolve 和 NDSolveValue。兩者僅在輸出格式上有細微差異，內(nèi)部處理則完全一致。因此，在下面的文本表述中，將用簡短的 "NDSolve"，而將便于輸出的"NDSolveValue" 用于代碼范例中。若要在 Mathematica 上使用 FEM，首先加載軟件包：

 

 

然后，我們要做的就是給 NDSolve 提供一個 PDE、區(qū)域和初始及邊界條件。以在單位圓上的泊松方程 –∇²u = 1 為例，如果以在 x>=0 上 u=0 作為邊界條件：

 

 

所得出解的圖形為：

 

2.1 輸入表達式

目前，在 NDSolve 中適用于有限元法的偏微分方程式必須具有以下形式：

 

 

此處，待求解的因變量 u 在 Rⁿ上為一維函數(shù)時，m、d、a、f 為標(biāo)量，α、γ 和 β 為 n 維向量，c 為 n*n 矩陣。此外，c、a、f、α、γ 和 β 可依賴 t、x ∈ Rⁿ、u、∇u 等，但 m 和 d 僅對 x 有依賴性。對于多個因變量 u ∈ R^d建立聯(lián)立方程式時，方程式 (1) 中，γ 和 f 為 d 維向量，其他系數(shù)是向量分量的矩陣。但是，當(dāng)進行系數(shù)矩陣或向量 u 的運算時，微分運算符會被矩陣/向量的各個分量內(nèi)的運算所繼承，例如，∇(u₁, …, u_d)^T = (∇u₁, …, ∇u_d)^T 和

 

等等。

從自然科學(xué)到工程應(yīng)用的多數(shù) PDE 都是 (1) 的特例。例如，波動方程是除 m 和 c 以外的系數(shù)都為 0 時的情況，且表示不可壓縮流體的速度以及壓力場 u = (v, p)^T ∈ R⁴ 的 Navier–Stokes 方程式

 

可用方程式 (1) 的形式表示。

 

 

下面，我們考慮的問題將暫時與時間無關(guān)，并處理與空間維數(shù)有關(guān)的有限元法．與時間有關(guān)的問題將在第 3 節(jié)末尾作簡要說明，并且在 4.3 和 4.4 節(jié)中給出范例。

重要的是，為了判斷 FEM 是否適用，提供給 NDSolve 的方程應(yīng)視為方程式 (1) 的形式 ("Coefficient Form")，包括各個系數(shù)的依賴性。舉一個簡單的例子，
 

該式相當(dāng)于方程式(1) 中 c = –∇u, f = –4，并且將其他系數(shù)設(shè)置為 0 的情況。但作為提供給 NDSolve 的 PDE 進行輸入時，

 

PDE 則在 NDSolve 開始處理前被計算，結(jié)果 2u´ (x)u´´ (x) 被視為方程式 (1) 的第一項．由于這并非是方程式 (1) 中 Coefficient Form 的形式，不能用 FEM 求解（u´´(x) 被視為 u´(x) 的系數(shù)，造成系數(shù)依賴于二階導(dǎo)數(shù)函數(shù)的結(jié)果）。為保持方程式 (1) 的形式不變提供給 NDSolve，需要使用 Inactive 或者 Inactivate，

 

如上所示，保留對 ∇ 的計算即可。

 

2.2 指定區(qū)域

您可以指定任何維數(shù)的任何區(qū)域。如果是像上述泊松方程示例那樣的簡單的圖形，則可以通過 Disk 和 Polygon 的組合來創(chuàng)建；如果它是由等式和不等式表示的區(qū)域，則可以使用 ParametricRegion、 ImplicitRegion 等。此外，還可以通過 ImageMesh 將由照片生成的圖像轉(zhuǎn)換成可供  NDSolve  使用的區(qū)域數(shù)據(jù)。

2.3 指定邊界條件

直接在邊界 ∂Ω 上給出函數(shù)值的狄里克雷邊界條件：

 

如上所示，只需用 PDE 指定即可。此處 bc 是給出邊界值的函數(shù)，predicate是 f(x, y)=bc 需要滿足的邊界。如果僅將 predicate 設(shè)置為 True，則將指定整個∂Ω。

廣義諾伊曼邊界條件（洛平邊界條件）由 NeumannValue 指定。洛平邊界條件以下列形式定義了垂直向外穿透邊界的通量分量：

為 ∂Ω 上的外向法線（單位）向量，右側(cè)的 g–qu 是由用戶給定的值。但請注意，NeumannValue 與 DirichletCondition 的指定方法不同。這是因為在有限元逼近中，PDE 乘以測試函數(shù) ? 并積分到區(qū)域 Ω 中以獲得弱形式。在等式（1）的第一項 ? 上積分，項則變?yōu)椋?/span>

 

在邊界 ∂Ω 上積分的被積函數(shù)剛好與在洛平邊界條件應(yīng)指定的值相對應(yīng)。因此，通過用 g–qu 的積分代替此項，NDSolve 則可正確處理該邊界條件。

例如，在區(qū)域 u(x,y) = 0, x ≥ 1/2 中，對于單位圓邊界 x ≤ 0，指定

·∇u = xy²的狄利克雷條件和諾伊曼條件，求解拉普拉斯方程 –∇²u = 0：

 

（在這種情況下，PDE 被明確識別為 Coefficient Form，因此不必將微分運算符設(shè)置為 Inactive，但也可以使用 Inactive[Div], Inactive[Grad]。）

Div 前的負號是為了讓等式(1) 中的 c 成為 1，從而讓諾伊曼條件

·c∇u =  ·∇u = xy² 可直接輸入 NeumannValue。繪制解的圖形：

 

需要注意的是，由于 NeumannValue 是根據(jù)方程(1) 的 PDE 以方程 (4) 的形式確定的，因此可能需要"手動"調(diào)整諾伊曼條件。例如，對于泊松方程  –∇²u + 1/5 = 0，同時指定諾伊曼條件·∇u = xy²(x ≥ 1/2) 和狄利克雷條件 u(x,y) = 0(x ≤ 0)，如果用于 NDSolve 的 PDE 本身為 –5 ∇²u + 1 = 0，則 NDSolve 的公式 (1) 中認定 c = 5，與 ·c∇u  對應(yīng)的 NeumannValue 必須為 5xy²?？雌饋硭坪鹾苊黠@，但需要注意的是，與狄利克雷條件不同，諾伊曼（洛平）條件并非獨立于 PDE 而指定。–∇²u + 1/5 = 0 和 –5 ∇²u + 1 = 0 的情況如下所示。

 

如果輸入式為 –∇²u + 1/5 == 0，則：

 

如果輸入式為 –5 ∇²u + 1 == 0，則：

對于洛平條件 3u + ∇u = xy² 的情況也是如此，對于 NDSolve 內(nèi)的 PDE，如果左側(cè)為–∇²u+1/5，則右側(cè)為 5*NeumannValue[1/2*x*y2 – 3/2*u[x, y], x≥1/2]。

3. Wolfram 語言中的非線性有限元法
 

為了應(yīng)用 FEM，需要在目標(biāo)區(qū)域中生成網(wǎng)格，在此不做深入討論。對于有興趣的讀者，在此介紹幾個簡單工具供了解：Triangle 和 TetGen（鏈接見文末）。前者用于生成二維區(qū)域，后者用于生成三維區(qū)域。Triangle 可進行 Delaunay 三角剖分、約束 Delaunay 三角剖分(constrained Delaunay triangulation) 和限定Delaunay三角剖分(conforming Delaunay triangulation)，TetGen 可以在三維區(qū)域中生成四面體網(wǎng)格，例如約束 Delaunay 四面體化和邊界符合Delaunay 網(wǎng)格。Wolfram 語言將在需要時自動使用，但您也可以對其進行自定義靈活使用。有關(guān) Triangle 的詳情請參照 [4.1] 節(jié)，TetGen 的相關(guān)說明請參見 [4.2] 節(jié)。

 

在線性 PDE 的情況下，聯(lián)立線性方程組是從 PDE 的弱形式到離散化來求解的，但這也用于求解非線性 PDE。以下為基本流程：

 

1. 在成為種子的候選解附近線性化非線性PDE

2. 對線性化方程進行離散化求解

3. 如果種子和所獲得的解的差異在允許的誤差內(nèi)，則結(jié)束

4. 使用獲得的解作為新種子，返回到第1步的線性化工作

 

也就是說，它遵循的過程與用 Newton-Raphson 方法求解非線性代數(shù)方程式的過程相同。在上面提到的 Wolfram 文檔教程有詳細介紹，以下為簡單介紹。

首先，如果我們刪除與公式(1) 的時間導(dǎo)數(shù)相關(guān)的部分，則有

 

若將，

 

 

則變?yōu)橐韵潞唵涡问剑?/span>

 

盡管將非線性 PDE 進行線性化，與求 1 個變量的非線性方程組的數(shù)值解相同，將任意函數(shù) u₀ 作為種子，由此漸進逼近使 ∇²·Γ (u^∗) – F (u^∗) = 0 的真正解 u^∗，令 u^∗與 u₀ 之間的差為 r = u^∗ – u₀，則得出：

 

將 Γ、F 在 u₀ 周圍通過泰勒展開，忽略 O(r²)  的高階項，則得出：

 

 

對 Γ ' 或 F ' 的微分求導(dǎo)，可通過計算 ∂Γ/∂∇u、∂Γ/∂u、∂F/∂∇u、∂F/∂u 等得出。當(dāng)它們在 u₀ 處求值時，等式(9) 成為每個離散點（節(jié)點）上 u 的聯(lián)立線性方程。在這里，通過同時聯(lián)立的初始條件和邊界條件，從而形成一個封閉的聯(lián)立方程并且得出 r。

種子 u₀ 默認為 u(x) = 0, ∀ x ∈ Ω，是 NDSolve 的一個選項，例如，可指定為 InitialSeeding→{u[x,y]==x+Exp[-Abs[y]]} 考慮到線性化的漸近解可能導(dǎo)致意想不到的局部解，提供與問題相符的好種子也是有用的。另外，從等式(13)計算殘差 r 時，左側(cè)出現(xiàn)的雅可比矩陣 ∇·Γ '(u₀) – F '(u₀) 的計算量很大，這極大地影響了整體計算時間。因此，在 Wolfram 語言中，當(dāng)應(yīng)用非線性 FEM 時，將使用仿射協(xié)變牛頓法（Affine Covariant Newton）代替 Newton-Raphson 法，并且在允許的范圍內(nèi)可以重復(fù)使用上一步中的雅可比法。從而顯著減少雅可比的計算次數(shù)。

對于時間相關(guān)的積分，可以通過離散化空間維度以獲得方程組（矩陣），然后將其作為關(guān)于時間的常微分方程，從而應(yīng)用各種計算方法。直線數(shù)值法 (Method of lines)等[5]，在某些情況下，也可以將 FEM 應(yīng)用于時間維度。

 

4. 實例應(yīng)用

4.1 非線性磁導(dǎo)率下的磁場分布

電流周圍會產(chǎn)生磁場 。當(dāng)電流以類似電動機的配置流過線圈時，會產(chǎn)生什么樣的磁場分布呢？對于非線性材料，其構(gòu)成定子和轉(zhuǎn)子的鐵磁材料的磁導(dǎo)率取決于磁場，我們將特別對此情況進行計算?；灸Ｐ凸绞谴艌雠c矢量電勢之間的關(guān)系 以及安培定律 。將它們組合在一起，并使用庫侖計 ，

 

假設(shè)電流僅在 z 方向上具有分量，為了簡化問題，假設(shè)矢量電勢的 x 和 y 分量為常數(shù)，并假設(shè)磁導(dǎo)率在 z 方向上也為常數(shù)。由此，在等式(10)中只有 z 分量是有效的，它是標(biāo)量  u = A_z 的 PDE：

 

 

對于磁導(dǎo)率 μ(B)，使用根據(jù)以下測量數(shù)據(jù)擬合的方程。

 

下圖顯示了電動機的橫截面示意圖，假設(shè)電流在黃色和橙色部分中沿垂直于屏幕的方向流動，則通過非線性 PDE 公式(11)計算磁場的強度分布。讓我們計算一下。

設(shè)定磁導(dǎo)率并指定電流元件規(guī)格。在 NDSolve 中施加狄里克雷邊界條件，即電動機定子外部的磁場為零。

 

將獲得的結(jié)果與電機結(jié)構(gòu)線框一起顯示：

 

4.2 驅(qū)動空腔問題

以下是描述穩(wěn)態(tài)下不可壓縮流體的 Navier-Stokes 方程：

 

 

上面第一個方程第二項中的對流項（在此，ρ = 1 且外力場為零）基本上是非線性的。此處，由于 u 是向量，如果是二維，則第一個方程式由兩個方程式 u_x 和 u_y 組成，微分算子∇作用于該方程式（請參見下面的代碼）。讓我們計算二維空腔中的速度場。狄里克雷條件是填充在空腔中的流體的上側(cè)以恒定速度 u_x 驅(qū)動到右側(cè)，剩余側(cè)的通量為零，并且圖形左下角的壓力為零。Wolfram 語言代碼如下：

 

可視化獲得的速度場：

壓力分布如下：

4.3 Gray-Scott 模型

由于化學(xué)反應(yīng)和物質(zhì)擴散而導(dǎo)致的多種物質(zhì)的濃度變化被稱為反應(yīng)擴散系統(tǒng)。以下以此為模型計算非線性聯(lián)立 PDE（Gray-Scott 模型）的示例。從外部將作為原料的化學(xué)物質(zhì) U 連續(xù)地引入填充有另一種物質(zhì) V 的反應(yīng)容器中，并進行自催化反應(yīng)。

 

 

之后，U 變成最終產(chǎn)品 P，然后將 P 排放到系統(tǒng)外部。U 和 V 濃度 u 和 v 隨時間變化所描述的就是這個模型：

 

 

D_u 和 D_v 是各自的擴散系數(shù)，F 是物質(zhì) U 的補充率，K 是定義反應(yīng)速度 V→P 的參數(shù)。從輸入表達式到可視化（動畫）的代碼如下所示：

4.4 繞圓柱體流動

通過將隨時間變化的項添加到公式(12)，來描述隨時間變化的不可壓縮流體流動。

 

對于流體在兩個無限寬的平行板間流動的狀態(tài)，讓我們計算一個無限長的圓柱體垂直于流向位于該空間中時的流速分布。未知數(shù)是速度(u,v)和壓力 p，其中垂直于板和圓柱的平面為 xy 平面。Wolfram 語言代碼如下。

Navier-Stokes 方程式：

 

設(shè)置入口處水池的大小和速度分布。定義 rampFunction，該函數(shù)可提供平滑的速度變化，以使速度在特定時間不會從零變?yōu)榉橇?。由于流域的大小和流體速度，此處的雷諾數(shù)約為 200。

 

邊界條件和初始條件：

 

速度分布從 t  = 0 到 10 的變化是由 NDSolve 在監(jiān)視 t 的同時計算的。在普通 PC (3.1GHz Intel Core i5，內(nèi)存16GB) 上大約需要運算 6 分鐘。

根據(jù)每個點的速度絕對值進行著色并創(chuàng)建動畫。

可以確認以下為該區(qū)域生成的網(wǎng)格：

 

5. 結(jié)束語

Mathematica 12（Wolfram語言 12）極大地擴展了有限元方法的應(yīng)用范圍，使得包括 Navier-Stokes 方程在內(nèi)的許多非線性偏微分方程的求解變?yōu)榭赡?。由?Wolfram 語言在符號計算方面的優(yōu)勢，無論 PDE 形式如何，都可以在保證求解的高效性和統(tǒng)一性的同時，保證其高度通用性。有關(guān) FEM 的內(nèi)部處理的詳細信息已經(jīng)發(fā)布。與此同時，Mathematica 教程也對彈性分析、聲學(xué)分析和熱/振動傳導(dǎo)分析等在許多領(lǐng)域的應(yīng)用實例有詳細說明，以供大家參考。

 

參考教程和工具：

 

有限元編程：https://reference.wolfram.com/language/FEMDocumentation/tutorial/FiniteElementProgramming.html

 

網(wǎng)格生成和三角剖分：

https://www.cs.cmu.edu/~quake/triangle.html

http://wias-berlin.de/software/index.jsp?id=TetGen