上次專欄寫了scikit-fem的兩個(gè)例子。但隨著后續(xù)使用發(fā)現(xiàn)，這個(gè)軟件包用起來也相當(dāng)不方便，最大的問題還是出在文檔不完整。我嘗試求解一個(gè)向量值的方程（這個(gè)專欄）就會(huì)講，但嘗試了好久也沒搞清楚在scikit-fem中要怎么實(shí)現(xiàn)。最后在github上問作者，他寫的代碼用到的是這個(gè)包中還處于試驗(yàn)階段的模塊NonlinearForm，還要依賴外部軟件包autograd，于是我瞬間就不想用了??。后來想想，其實(shí)一維方程解起來并不麻煩。有限元方法麻煩的是在高維以及高階基函數(shù)的情形，此時(shí)需要規(guī)劃產(chǎn)生的線性方程組矩陣元的排列順序。因?yàn)槲矣玫降闹挥幸痪S情形，那為什么不自己寫呢。

1. 簡(jiǎn)單的擴(kuò)散方程

   下面我們求解方程

   

   初始條件為u=1。

使用向后歐拉公式：

兩邊同時(shí)乘以$v$，再積分得到變分形式

一般trial function和test function取成一樣的基函數(shù)，注意只有當(dāng)test function 時(shí)，，所以這一項(xiàng)只出現(xiàn)在第一個(gè)方程里。

我們可以先計(jì)算出 和的矩陣然后求解。注意到，因?yàn)橛疫吔鐥l件是不需要求解的，所以未知數(shù)的個(gè)數(shù)是n個(gè)，從C0到C(n-1)。

我們使用一維的線性元。假設(shè)把區(qū)間分解成。則相應(yīng)的基函數(shù)為

i=1,2,...,n-1

在使用這個(gè)基組時(shí)，每個(gè)小區(qū)間$[x_{i-1}, x_i]$上有兩個(gè)基函數(shù)：和。把它們進(jìn)行如下的線性組合

所以把基函數(shù)按照函數(shù)值線性組合得到的就是線性插值的結(jié)果。

下面是完整的代碼：

2. 非線性的擴(kuò)散方程組

邊界條件為

同樣使用向后歐拉公式，得到方程的變分形式。然后把解用基函數(shù)展開。這里只給出展開后的結(jié)果。

這是關(guān)于和的非線性方程，要使用牛頓法求解。但這里我們不寫出它的Jacobian矩陣表達(dá)式，因?yàn)樘÷?！煩！了！。直接把它?dāng)作一個(gè)非線性方程交給scipy的newton_krylov函數(shù)求解即可。完整代碼如下：