有限元法解常微分方程边值问题
评论
有限元法可归结为如下几个步骤：
- 转换成变分问题（应该会用到边界条件）
2）对解域进行刨分（可以是不均匀）
3）构造基函数（本篇采用基于线性插值的基函数）
4）推导出有限元方程
5）求解有限元方程
6）收敛性和误差估计
# 奇怪：这个库必须放在最前面才能一次加载成功
using Plots
gr()
# 解决Colab不显示输出数学公式的问题
using Markdown: MD, LaTeX
function latex(expr)
expr |> sympy.latex |> LaTeX
end;
考虑两点边值问题：
\[ \begin{aligned} &Lu=-\frac{d}{dx}\left(p\frac{du}{dx}\right)+qu=f ,\qquad a < x < b \\ & u(a)=0,\frac{d}{dx}u(b)=0 \end{aligned} \]对应的变分方程(用到两个边界条件)：
\[ \begin{aligned}&a(u,v)=(f,v),\qquad \forall v \\ & a(u,v)=\int^b_a{\left(p\frac{du}{dx}\frac{dv}{dx}+quv\right)dx}\end{aligned} \]区间剖分（可以不均匀）：
\[ a = x_0 < x_1 < \dots < x_n=b \\ h_i=x_i-x_{i-1}, \qquad h = \mathrm{max} \{h_i\} \]与此对应的近似解序列 \(\{u_i\}\) （待求）为：
\[ u_0=0, u_1,u_2,\dots,u_n \]通过这 \(n+1\) 个点的值可进行线性插值得到近似解 \(u_h(x)\) ：
\[ \begin{aligned}u_h(x)=\frac{x_i-x}{h_i}u_{i-1}+\frac{x-x_{i-1}}{h_i}u_i,\qquad x_{i-1} < x < x_i \end{aligned} \]可引入线性无关的“山形函数”序列 \(\varphi_i(x)\) 构成基底：
不难验证，前面的线性插值函数 \(u_h(x)\) ，可用这个序列 \(\varphi_i\) 为基底展开(下面考虑了 \(u_0=0\) )：
\[ u_h(x)=\sum^n_{i=1}{u_i \varphi_i(x)} \]由图可知（亦可简单验算）：
\[ \varphi_i(x)\varphi_j(x)=0, \qquad |i-j|\ge 2 \\ \frac{d\varphi_i}{dx}\frac{d\varphi_j}{dx}=0, \qquad |i-j|\ge 2 \]所以 \(a(\varphi_i,\varphi_j)\) 的非零值只可能有：
\[ \begin{aligned}&a(\varphi_{1},\varphi_1),a(\varphi_{2},\varphi_1) \\ &a(\varphi_{j-1},\varphi_j),a(\varphi_{j},\varphi_j),a(\varphi_{j+1},\varphi_j),\qquad j=2,\dots,n-1\\ & a(\varphi_{n-1},\varphi_n),a(\varphi_{n},\varphi_{n}) \end{aligned} \]可分别算出这三类非零值：
\[ a(\varphi_{j-1},\varphi_j)=\int^1_0{\left[-h_j^{-1}p(x_{j-1}+h_j\xi)+h_j q(x_{j-1}+h_j\xi)(1-\xi)\xi\right]d\xi} \] \[ \begin{aligned}a(\varphi_j,\varphi_j)&=\int^1_0{\left[h_j^{-1}p(x_{j-1}+h_j\xi)+h_j q(x_{j-1}+h_j\xi)\xi^2\right]d\xi} \\ & +\int^1_0{\left[h_{j+1}^{-1}p(x_{j}+h_{j+1}\xi)+h_{j+1} q(x_{j} +h_{j+1}\xi)(1-\xi)^2\right]d\xi}\end{aligned} \] \[ a(\varphi_{j+1},\varphi_j)=\int^1_0{\left[-h_{j+1}^{-1}p(x_{j}+h_{j+1}\xi)+h_{j+1} q(x_{j}+h_{j+1}\xi)(1-\xi)\xi\right]d\xi} \] \[ a(\varphi_n,\varphi_n)=\int^1_0{\left[h_n^{-1}p(x_{n-1}+h_n\xi)+h_n q(x_{n-1}+h_n\xi)\xi^2\right]d\xi} \]using SciPy
#
# “山形”基函数
#
function φ(t,i,x,h,n)
if i > 1 && t>=x[i-1] && t<=x[i]
return (t-x[i-1])/h[i-1]
elseif i < n+2 && t>x[i] && t<=x[i+1]
return 1-(t-x[i])/h[i]
else
return 0
end
end;
#
# a(φ_i,φ_j)
#
# 假设已有p(t),q(t)
#
function aa(i,j,x,h,n)
if abs(i-j) >= 2
return 0;
end;
if i==j-1
α(ξ)=-p(x[j-1]+h[j-1]*ξ)/h[j-1]+h[j-1]*q(x[j-1]+h[j-1]*ξ)*(1-ξ)*ξ;
return SciPy.integrate.quad(α, 0, 1)[1];
end;
if i==j && j<n+1
β(ξ)=p(x[j-1]+h[j-1]*ξ)/h[j-1]+h[j-1]*q(x[j-1]+h[j-1]*ξ)*ξ^2+p(x[j]+h[j]*ξ)/h[j]+h[j]*q(x[j]+h[j]*ξ)*(1-ξ)^2;
return SciPy.integrate.quad(β, 0, 1)[1];
end;
if i==j && j==n+1
βn(ξ)=p(x[j-1]+h[j-1]*ξ)/h[j-1]+h[j-1]*q(x[j-1]+h[j-1]*ξ)*ξ^2;
return SciPy.integrate.quad(βn, 0, 1)[1];
end;
if i==j+1
γ(ξ)=-p(x[j]+h[j]*ξ)/h[j]+h[j]*q(x[j]+h[j]*ξ)*(1-ξ)*ξ;
return SciPy.integrate.quad(γ, 0, 1)[1];
end;
end;
正式求解的julia代码：
using Random
# 参数选定
a=0; b=1
f(t)=(pi^2/2)sin(pi*t/2);
p(t) = 1;
q(t) = pi^2/4;
function NSolve(n)
# 随机生成的不均匀的区间刨分
h = [0.5*rand()+0.5 for i in 1:n];
h = (b-a)*h/sum(h);
x = zeros(n+1);
x[1]=a; x[n+1]=b;
x[2:n]=[sum(h[1:i]) for i in 1:(n-1)];
h[n]=x[n+1]-x[n];
# 待求近似解
u = zeros(n+1);
# 计算a(φ_i,φ_j)矩阵
A = [aa(i,j,x,h,n) for j in 2:n+1, i in 2:n+1];
# 计算(f,φ_j)列向量
fφ = [SciPy.integrate.quad(t->f(t)*φ(t,j,x,h,n),
x[j-1],
(j==n+1 ? b : x[j+1])
)[1] for j in 2:n+1];
# 解出插值端点值
u[2:n+1] = inv(A)*fφ;
# 近似解序列
(x,u)
end;
做为对比，也可解出精确解：
using SymPy
@vars t
y = SymFunction("y")
diffeq = -diff(y(t),t,2)+(PI^2/4)y(t)⩵(PI^2/2)*sin(t*PI/2)
# 通解
ex = dsolve(diffeq, y(t))
# 根据边界条件确定积分常数
ex1 = rhs(ex)
ex2 = diff(ex1,t)
eqs = [ex1(t=>0)⩵0,ex2(t=>1)⩵0]
C = solve(eqs,[Sym("C1"),Sym("C2")])
# 解析解
ex = ex(C)
ex |> latex |> MD
绘制近似解和精确解
# 精确解
plot(rhs(ex),color="black",xlims=(-0.1, 1.1),label="U",xlabel="x",ylabel="u(x)")
# 近似解 n=4,8
plot!(NSolve(4),color="green",label="u4")
plot!(NSolve(8),color="red",label="u8")
《微分方程数值解法（第4版）-李荣华＆刘播-高等教育出版社-2009》