有限元方法軟件scikit-fem的兩個(gè)例子

2023-07-08 10:07 作者:quantumtopology 0人讀過 | 我要投稿

最近在使用有限元方法。有限元軟件有很多，最出名的應(yīng)該是fenics，它功能非常多，功能包裝得也很好，使用方便。唯一的缺點(diǎn)是太大了，如果只是寫個(gè)小軟件去安裝它就顯得非常臃腫。想象一下，你使用外部庫比你自身的程序大的多，而這個(gè)外部庫還是不常用的。fenics在wsl上還不太好安裝。因?yàn)閣indows下認(rèn)為兩個(gè)名字相同但大小寫不同的文件是同一文件，但安裝fenics時(shí)需要安裝一個(gè)和已有庫同名的軟件包，所以就安裝不上。后來找來找去找到一個(gè)名為scikit-fem的包。它很小，看它官網(wǎng)上幫助文檔發(fā)現(xiàn)也能解很多問題。但缺點(diǎn)是這幫助文檔寫的太簡略了！有好多代碼和函數(shù)需要猜它什么意思。這了記錄兩個(gè)典型例子，可能會(huì)對使用scikit-fem的人有所幫助。

一維波動(dòng)方程

該例子是ex44.py中的例子。對于方程 $u_%7Btt%7D%20%3D%20c%5E2u_%7Bxx%7D$ ，它需要轉(zhuǎn)化成一個(gè)一階方程組求解。所以這個(gè)例子對于求解方程組也很有幫助。

$%5Cbegin%7Bcases%7D%0A%5Cfrac%7B%5Cpartial%20u_1%7D%7B%5Cpartial%20t%7D%20%3D%20-u_2%20%5C%5C%0A%5Cfrac%7B%5Cpartial%20u_2%7D%7B%5Cpartial%20t%7D%20%3D%20-c%5E2%20%5Cfrac%7B%5Cpartial%5E2%20u_1%7D%7B%5Cpartial%20x%5E2%7D%0A%5Cend%7Bcases%7D$

使用Crank-Nicolson差分格式

$%5Cbegin%7Bcases%7D%0Au_1%5E%7Bk%2B1%7D%20-%20u_1%5Ek%20%26%3D%20-0.5%5CDelta%20t(u_2%5E%7Bk%2B1%7D%2Bu_2%5Ek)%20%5C%5C%0Au_2%5E%7Bk%2B1%7D%20-%20u_2%5Ek%20%26%3D%20-0.5c%5E2%5CDelta%20t%20%5Cleft(%20%5Cfrac%7B%5Cpartial%5E2%20u_1%5E%7Bk%2B1%7D%7D%7B%5Cpartial%20x%5E2%7D%20%20%20%2B%20%5Cfrac%7B%5Cpartial%5E2%20u_1%5E%7Bk%7D%7D%7B%5Cpartial%20x%5E2%7D%20%20%20%5Cright)%0A%5Cend%7Bcases%7D$

把其中第二個(gè)方程轉(zhuǎn)化成變分形式（第一個(gè)方程不用轉(zhuǎn)化，因?yàn)闆]有導(dǎo)數(shù)項(xiàng)），并整理得到

$%5Cbegin%7Bcases%7D%0Au_1%5E%7Bk%2B1%7D%2B0.5%5CDelta%20t%20u_2%5E%7Bk%2B1%7D%20%26%3D%20u_1%5E%7Bk%7D%20-%200.5%5CDelta%20t%20u_2%5Ek%20%5C%5C%0A%5Cint%20u_2%5E%7Bk%2B1%7D%20v%20-%200.5c%5E2%5CDelta%20t%5Cint%20%5Cfrac%7B%5Cpartial%20u_1%5E%7Bk%2B1%7D%7D%7B%5Cpartial%20x%7D%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20x%7D%20%26%3D%20%5Cint%20u_2%5Ek%20v%20%2B%200.5c%5E2%5CDelta%20t%20%5Cint%20%5Cfrac%7B%5Cpartial%20u_1%5E%7Bk%7D%7D%7B%5Cpartial%20x%7D%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20x%7D%0A%5Cend%7Bcases%7D$

現(xiàn)在看如下代碼塊

可以看到矩陣A和B的形式分別為

$%5Cbegin%7Bpmatrix%7D%0A%26I%20%260.5%5Cmathrm%7Bd%7DtI%20%5C%5C%0A%26-c%5E2L%200.5%5Cmathrm%7Bd%7Dt%20%26M%0A%5Cend%7Bpmatrix%7D%2C%20%5Cquad%0A%5Cbegin%7Bpmatrix%7D%0A%26I%20%26-0.5%5Cmathrm%7Bd%7DtI%20%5C%5C%0A%26c%5E2L%200.5%5Cmathrm%7Bd%7Dt%20%26M%0A%5Cend%7Bpmatrix%7D$

分別對應(yīng)上面推導(dǎo)公式的右端和左端。時(shí)間t^k時(shí)刻的u乘以B，得到Bu^k，然后求解Au^{k+1} = Bu^k，得到u^{k+1}。

2. 非線性方程

?非線性方程轉(zhuǎn)化成變分形式會(huì)得到非線性泛函。此時(shí)把用基組展開的試探函數(shù)帶入到變分形式中將得到關(guān)于展開系數(shù)的非線性方程。此時(shí)用牛頓迭代法去求解即可。在使用scikit-fem時(shí)，我們需要把變分形式線性化交給程序求解。一個(gè)泛函 $F%5Bu%5D$ 泛函導(dǎo)數(shù)的定義為

$%5Clim_%7Bt%5Crightarrow%200%7D%20%5Cfrac%7BF%5Bu%2Bt%5Cdelta%20u%5D%20-%20F%5Bu%5D%7D%7Bt%7D%3D%5Cfrac%7B%5Cdelta%20F%7D%7B%5Cdelta%20u%7D%5Cdelta%20u$

引入?yún)?shù)t是為了能和普通函數(shù)求導(dǎo)一樣求出泛函導(dǎo)數(shù)，其實(shí)是不必要的。這里選用的例子是ex10.py。對于極小曲面問題（幫助文檔里寫錯(cuò)了，掉了根號）

$%5Cnabla%5Ccdot%5Cleft(%20%5Cfrac%7B%5Cnabla%20u%7D%7B%5Csqrt%7B1%20%2B%20%7C%7C%20%5Cnabla%20u%20%7C%7C%5E2%7D%7D%5Cright)%20%3D%200$

它的變分形式為

$F(u%2Cv)%20%3D%20%5Cint%20%5Cfrac%7B%5Cnabla%20u%5Ccdot%20%5Cnabla%20v%7D%7B%5Csqrt%7B1%20%2B%20%7C%7C%5Cnabla%20u%7C%7C%5E2%7D%7D%20%5Cmathrm%7Bd%7D%5COmega%20%20%3D%200$

為了求解泛函導(dǎo)數(shù)，把如下泛函在t = 0處展開到一階（為了方便省略積分號）

$%5Cfrac%7B%5Cnabla%20(u%20%2B%20t%5Cdelta%20u)%20%5Cnabla%20v%7D%7B%5Csqrt%7B1%20%2B%20%7C%7C%5Cnabla%20(u%20%2B%20t%5Cdelta%20u)%7C%7C%5E2%7D%7D%20%3D%20%5Cfrac%7B%5Cnabla%20u%20%5Cnabla%20v%7D%7B%5Csqrt%7B1%20%2B%20%7C%7C%5Cnabla%20u%7C%7C%5E2%7D%7D%20%2B%20J%5Ccdot%20t$

$J%3D%20%5Cfrac%7B%5Cnabla%20(u%20%2B%20%5Cdelta%20u)%20%5Cnabla%20v%7D%7B%5Csqrt%7B1%20%2B%20%7C%7C%5Cnabla%20(u%20%2B%20t%5Cdelta%20u)%7C%7C%5E2%7D%7D%20%0A-%20%5Cfrac%7B%5Cnabla%20(u%20%2B%20t%5Cdelta%20u)%20%5Cnabla%20v%7D%7B(1%20%2B%20%7C%7C%5Cnabla%20(u%20%2B%20t%5Cdelta%20u)%7C%7C%5E2)%5E%7B3%2F2%7D%7D%20%5Cfrac%7B1%7D%7B2%7D(2%5Cnabla%20u%5Ccdot%20%5Cnabla%20%5Cdelta%20u%2B2t%20%5Cnabla%20%5Cdelta%20u%5Ccdot%20%5Cnabla%20%5Cdelta%20u)$

注意這里使用了

$%7C%7C%5Cnabla(u%20%2B%20t%5Cdelta%20u)%7C%7C%5E2%20%3D%20%5Cnabla(u%20%2B%20t%5Cdelta%20u)%5Ccdot%20%5Cnabla(u%20%2B%20t%5Cdelta%20u)%5C$

當(dāng)t趨于0時(shí)

$J(t%5Crightarrow%200)%3D%20%5Cfrac%7B%5Cnabla%20(u%20%2B%20%5Cdelta%20u)%20%5Cnabla%20v%7D%7B%5Csqrt%7B1%20%2B%20%7C%7C%5Cnabla%20u%7C%7C%5E2%7D%7D%20%0A-%20%5Cfrac%7B%5Cnabla%20u%20%5Cnabla%20v%7D%7B(1%20%2B%20%7C%7C%5Cnabla%20u%7C%7C%5E2)%5E%7B3%2F2%7D%7D%20%5Cfrac%7B1%7D%7B2%7D(2%5Cnabla%20u%5Ccdot%20%5Cnabla%20%5Cdelta%20u)$

看ex10.py中的代碼塊

這里w = p['prev']就是前一步的解，即公式中的u。代碼中的u則是上式中的 $%5Cdelta%20u$ ， $u%5E%7Bk%2B1%7D%20%3D%20u%5Ek%20%2B%20%5Cdelta%20u$ ?？梢钥吹絁acobian恰好是和推導(dǎo)得到的泛函導(dǎo)數(shù)相對應(yīng)。