title: "有限元法解常微分方程边值问题"
date: 2020-11-02T12:25:49+08:00
draft: true
categories: ["科学计算"]
tags: ["julia","有限元法","微分方程","数值计算"]

有限元法可归结为如下几个步骤：

1. 转换成变分问题（应该会用到边界条件）

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)$构成基底：

![[0166.jpg) ]]

$$ \boxed{\xi = \frac{x-x_{i-1}}{h_i}\qquad x_{i-1} < x < x_i\qquad 0<\xi<1} $$

$$ \varphi_0(x)=\left{\begin{aligned}&1-\xi, &x_0 < x < x_1 \\ &0, & 其它\end{aligned}\right . $$

$$ \varphi_i(x)=\left{\begin{aligned}&\xi, &x_{i-1} < x < x_i \\ &1-\xi,& x_i < x < x_{i+1} \\ &0, & 其它\end{aligned}\right . \qquad i=1,2,\dots,n-1 $$

$$ \varphi_n(x)=\left{\begin{aligned}&\xi, &x_{n-1} < x < x_n \\ &0, & 其它\end{aligned}\right . $$

不难验证，前面的线性插值函数$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

$$ y{\left(t \right)} = \sin{\left(\frac{\pi t}{2} \right)} $$

绘制近似解和精确解

# 精确解
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")

![[0162.svg]]

《微分方程数值解法（第4版）-李荣华＆刘播-高等教育出版社-2009》