有限元法解常微分方程边值问题

有限元法解常微分方程边值问题

2020-11-02
| 科学计算 | | julia , 有限元法 , 微分方程 , 数值计算 | Comment 评论

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

  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)列向量= [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》