手动实现用有限元方法求解一维方程
上次专栏写了scikit-fem的两个例子。但随着后续使用发现，这个软件包用起来也相当不方便，最大的问题还是出在文档不完整。我尝试求解一个向量值的方程（这个专栏）就会讲，但尝试了好久也没搞清楚在scikit-fem中要怎么实现。最后在github上问作者，他写的代码用到的是这个包中还处于试验阶段的模块NonlinearForm，还要依赖外部软件包autograd，于是我瞬间就不想用了😂。后来想想，其实一维方程解起来并不麻烦。有限元方法麻烦的是在高维以及高阶基函数的情形，此时需要规划产生的线性方程组矩阵元的排列顺序。因为我用到的只有一维情形，那为什么不自己写呢。
简单的扩散方程
下面我们求解方程
初始条件为u=1。
使用向后欧拉公式：
两边同时乘以$v$，再积分得到变分形式
一般trial function和test function取成一样的基函数，注意只有当test function 时，
，所以
这一项只出现在第一个方程里。
我们可以先计算出 和
的矩阵然后求解。注意到，因为右边界条件是不需要求解的，所以未知数的个数是n个，从C0到C(n-1)。
我们使用一维的线性元。假设把区间分解成。则相应的基函数为
i=1,2,...,n-1
在使用这个基组时，每个小区间$[x_{i-1}, x_i]$上有两个基函数：和
。把它们进行如下的线性组合
所以把基函数按照函数值线性组合得到的就是线性插值的结果。
下面是完整的代码：
2. 非线性的扩散方程组
边界条件为
同样使用向后欧拉公式，得到方程的变分形式。然后把解用基函数展开。这里只给出展开后的结果。
这是关于和
的非线性方程，要使用牛顿法求解。但这里我们不写出它的Jacobian矩阵表达式，因为太！麻！烦！了！。直接把它当作一个非线性方程交给scipy的newton_krylov函数求解即可。完整代码如下：
