Galerkin法解常微分方程边值问题

Galerkin法解常微分方程边值问题

2020-11-01
| 科学计算 | | julia , Galerkin , 微分方程 , 数值计算 | Comment 评论 
# 奇怪：这个库必须放在最前面才能一次加载成功
using Plots
gr()

# 解决Colab不显示输出数学公式的问题
using Markdown: MD, LaTeX
function latex(expr)
    expr |> sympy.latex |> LaTeX
end;

选择三角基函数

\[ \begin{aligned} & u''+u=-x \qquad 0 < x < 1 \\ &u(0)=u(1)=0 \end{aligned} \]

用微分算符表示之：

\[ \begin{aligned} &\boxed{Lu=f} \\ \\ &Lu=-\frac{d^2 u}{dx^2}-u \\ &f = x \end{aligned} \]

根据虚功原理，改写成变分等式：

\[ \begin{aligned} &\boxed{a(u,v)=(f,v)} \qquad \forall v \\ \\ &a(u,v)=\int^1_0{\left(\frac{d u}{dx}\frac{d v}{dx}-u v\right)dx} \\ &(f,v) = \int^1_0{x vdx} \end{aligned} \]

选择一族基函数（比如选择三角函数）：

\[ \varphi_i(x)=\sin(i \pi x) \qquad i=1,\dots,n \]

于是可用 \(u_n(x)\) 来逼近 \(u(x)\) :

\[ u_n(x)=\sum^n_{i=1}{c_i \varphi_i(x)} \]

\(v(x)\) 不妨选择 \(\varphi_j(x)\quad j=1,\dots,n\) , 于是得到一个线性方程组：

\[ \sum^n_{i=1}{a(\varphi_i,\varphi_j) c_j}=(f,\varphi_j)\qquad j=1,\dots,n \]

解出系数 \(c_j\) 即求解完毕。

用Julia实现如下：

using SciPy
using LinearAlgebra

# a(φ_i,φ_j)的被积函数
F(x,i,j) = i*j*pi^2*cos(i*pi*x)*cos(j*pi*x)-sin(i*pi*x)*sin(j*pi*x)
# (f,φ_j)的被积函数
f(x,j) = x*sin(j*pi*x)

# 自变量序列
x = [(j-1)*0.01 for j in 1:101];

function NSolve(n)
    
    # a(φ_i,φ_j) 组成的矩阵
    A = [SciPy.integrate.quad(x->F(x,i,j), 0, 1)[1] for i in 1:n, j in 1:n];

    # (f,φ_j)组成的数组
    b = [SciPy.integrate.quad(x->f(x,j), 0, 1)[1] for j in 1:n];

    # 解出系数c_j
    c = inv(A)*b;
    
    # 返回近似解
    [dot(c,[sin(i*pi*x[j]) for i in 1:n]) for j in 1:101]
end;

依次求 \(n=1,2,3,4\) 下的近似解：

plot(x,NSolve(1),color="red",label="u1",xlabel="x",ylabel="u(x)")
plot!(x,NSolve(2),color="green",label="u2")
plot!(x,NSolve(3),color="blue",label="u3")
plot!(x,NSolve(4),color="black",label="u4")

0155.svg

此图表明，随着n的变大，曲线变换越小，解就应该越精确。

作为对照，我们还可解出这个方程的解析解：

using SymPy

@vars t
y = SymFunction("y")
diffeq = diff(y(t),t,2)+y(t)⩵-t

# 通解
ex = dsolve(diffeq, y(t))
ex |> latex |> MD
\[ y{\left(t \right)} = C_{1} \sin{\left(t \right)} + C_{2} \cos{\left(t \right)} - t \] 
# 根据边界条件确定积分常数
eqs = [ex(t=>0,y(0)=>0),ex(t=>1,y(1)=>0)]
C = solve(eqs,[Sym("C1"),Sym("C2")])

Dict{Any,Any} with 2 entries:
  C1 => 1/sin(1)
  C2 => 0

# 解析解
ex(C) |> latex |> MD
\[ y{\left(t \right)} = - t + \frac{\sin{\left(t \right)}}{\sin{\left(1 \right)}} \] 
plot(x,NSolve(4),color="black",label="u4",xlabel="x",ylabel="u(x)")
plot!(rhs(ex)(C),xlims = (0, 1), ylims = (0, 0.08),color="red",label="U")

0156.svg

这个图表明 \(n=4\) 的近似解已经很接近解析解。

当 \(n=10\) 时，近似解和解析解在图中无法区分:

plot(x,NSolve(10),color="black",label="u10",xlabel="x",ylabel="u(x)")
plot!(rhs(ex)(C),xlims = (0, 1), ylims = (0, 0.08),color="red",label="U")

0157.svg

选择多项式基函数

\[ \begin{aligned} & u''+u=x^2 \qquad 0 < x < 1 \\ &u(0)=0,\quad u(1)=1 \end{aligned} \]

先做变换 \(u\rightarrow u+x\) ，得到0边界条件，解完再修正即可：

\[ \begin{aligned} & u''+u=x^2-x \qquad 0 < x < 1 \\ &u(0)=0,\quad u(1)=0 \end{aligned} \]

改写成变分等式：

\[ \begin{aligned} &\boxed{a(u,v)=(f,v)} \qquad \forall v \\ \\ &a(u,v)=\int^1_0{\left(\frac{d u}{dx}\frac{d v}{dx}-u v\right)dx} \\ &(f,v) = \int^1_0{(x-x^2) vdx} \end{aligned} \]

选择一族基函数（比如选择多项式函数,还必须注意满足边界条件）：

\[ \varphi_i(x)=x(1-x)x^{i-1} \qquad i=1,\dots,n \]

求解和前面类似线性方程组即可。

用Julia实现如下：

# a(φ_i,φ_j)的被积函数
F(x,i,j) =(i-(i+1)*x)*(j-(j+1)*x)*x^(i+j-2)-x^2*(1-x)^2*x^(i+j-2)
# (f,φ_j)的被积函数
f(x,j) = (x-x^2)*x*(1-x)*x^(j-1)

function NSolve(n)
    
    # a(φ_i,φ_j) 组成的矩阵
    A = [SciPy.integrate.quad(x->F(x,i,j), 0, 1)[1] for i in 1:n, j in 1:n];

    # (f,φ_j)组成的数组
    b = [SciPy.integrate.quad(x->f(x,j), 0, 1)[1] for j in 1:n];

    # 解出系数c_j
    c = inv(A)*b;
    
    # 返回近似解(注意：u+x才是目标解)
    [dot(c,[x[j]*(1-x[j])*x[j]^(i-1) for i in 1:n]) for j in 1:101]+x
end;

作为对比，同时解出精确解：

diffeq = diff(y(t),t,2)+y(t)⩵t^2

# 通解
ex = dsolve(diffeq, y(t))

# 根据边界条件确定积分常数
eqs = [ex(t=>0,y(0)=>0),ex(t=>1,y(1)=>1)]
C = solve(eqs,[Sym("C1"),Sym("C2")])

# 解析解
ex(C) |> latex |> MD
\[ y{\left(t \right)} = t^{2} + \frac{2 \left(1 - \cos{\left(1 \right)}\right) \sin{\left(t \right)}}{\sin{\left(1 \right)}} + 2 \cos{\left(t \right)} - 2 \] 
plot(rhs(ex)(C),lims=(-0.2,1.2),color="red",lw=2,label="U",xlabel="x",ylabel="u(x)")
plot!(x,NSolve(10),color="black",lw=1,label="u")

0158.svg

谱方法

\[ \begin{aligned} & -u''+u=2x\sin x-2\cos x \qquad 0 \le x \le 2\pi \\ &u(0)=u(2\pi)=0 \end{aligned} \]

设：

\[ \begin{aligned} & Lu=-u''+u\\ & f=2x\sin x-2\cos x \end{aligned} \] 
f(x)=2*x*sin(x)-2*cos(x);   # 输入可以是数值也可是符号

按前面类似得方法（有点小不同），可写出：

\[ \begin{aligned} &\boxed{a(u,v)=(f,v)} \qquad \forall v \\ \\ &a(u,v)=\int^{2\pi}_0{\left(\frac{d u}{dx}\frac{d \bar{v}}{dx}+u \bar{v}\right)dx} \overset{\Delta}{=}\int^{2\pi}_0{\alpha dx}\\ &(f,v) = \int^{2\pi}_0{f(x) \bar{v}dx} \end{aligned} \] 
@vars x
α(u,v)=diff(u,x)*diff(conj(v),x)+u*conj(v);

选择一族基函数：

\[ \varphi_k(x)=e^{i k x}-1 \qquad k=\pm 1,\dots,\pm n \]

那么解可用下式逼近：

\[ u_n(x)=\sum^n_{k=-n \\ k\ne 0}{c_k \varphi_k(x)}=\sum^n_{k=-n \\ k\ne 0}{c_k (e^{i k x}-1)}=\boxed{\sum^n_{k=-n}{c_k e^{i k x}}} \\ c_0=-\sum^n_{k=-n \\ k\ne 0}{c_k} \]

这就是所谓“谱方法”的命名由来。

进而：

\[ a(u_n,\varphi_j)=\sum^n_{k=-n}{c_k a(e^{ikx},e^{ijx})}-\sum^n_{k=-n}{c_k a(e^{ikx},1)} \]

求出 \(a(e^{ikx},e^{ijx})\) 分别在 \(k=j,k\ne j\) 情况下的值：

@vars x real=true k j integer=true

# 当k=j时
SymPy.integrate(α(exp(IM*k*x),exp(IM*j*x))(k=>j) |> simplify, (x, 0, 2*PI)) |> latex |> MD
\[ 2 \pi \left(j^{2} + 1\right) \] 
# 当k!=j时
SymPy.integrate(α(exp(IM*k*x),exp(IM*j*x)) |> simplify, (x, 0, 2*PI)) |> latex |> MD
\[ 0 \]

这意味着：

\[ a(u_n,\varphi_j)=c_j a(e^{ijx},e^{ijx})-c_0 a(1,1)=\boxed{2\pi(j^2+1)c_j-2\pi c_0} \]

此外还可计算出 \((f,\varphi_j)\) ：

fφ = SymPy.integrate(f(x)*conj(exp(IM*j*x)-1), (x, 0, 2*PI)) .|> simplify |> latex |> MD
\[ \begin{cases} \pi \left(1 + 2 i \pi\right) & \text{for}\: j = -1 \\\pi \left(1 - 2 i \pi\right) & \text{for}\: j = 1 \\\frac{4 \pi j^{2}}{j^{2} - 1} & \text{otherwise} \end{cases} \]

最后联立方程组

\[ \begin{aligned} &a(u_n,\varphi_j)=(f,\varphi_j) \qquad j=\pm 1,\dots,\pm n \\ &\sum^n_{k=-n}{c_k}=0 \end{aligned} \]

求出系数 \(c_j\) ：

# 自变量序列
x=2*pi*[(i-1)*0.01 for i in 1:101];

function NSolve(n)
    
    # 给矩阵分配空间
    A = zeros(2*n+1,2*n+1);
    b = (1.0+0.0*im)*zeros(2*n+1);

    # (f,φ_j)
    b[1] = pi+2*im*pi^2;
    b[2] = pi-2*im*pi^2;
    for j in 2:n
        b[2*j-1] = 4*pi*j^2/(j^2-1);
        b[2*j] = b[2j-1];
    end
    # a(u_n,φ_j)
    for j in 1:n
        A[2*j-1,2*j-1] = 2*pi*(j^2+1);
        A[2*j,2*j] = A[2*j-1,2*j-1];
    end
    A[1:2*n,2*n+1] = -2*pi*ones(2*n);
    # c_k
    A[2*n+1,:] = ones(2*n+1);
    
    # 解出系数c_j
    c = inv(A)*b
    
    # 返回近似解
    real([sum(c[2*k-1]*(exp(-im*k*x[i])-1)+c[2*k]*(exp(im*k*x[i])-1) for k in 1:n) for i in 1:101])
end;

作为对比，同时解出精确解：

@vars t
y = SymFunction("y")

diffeq = -diff(y(t),t,2)+y(t)⩵2*t*sin(t)-2*cos(t)

# 通解
ex = dsolve(diffeq, y(t))

# 根据边界条件确定积分常数
eqs = [ex(t=>0,y(0)=>0),ex(t=>2*PI,y(2*PI)=>0)]
C = solve(eqs,[Sym("C1"),Sym("C2")])

# 解析解
ex(C) |> latex |> MD
\[ y{\left(t \right)} = t \sin{\left(t \right)} \]

先绘制精确解：

plot(rhs(ex(C)),color="black",xlims=(-0.5, 6.5),label="U",xlabel="x",ylabel="u(x)")

0159.svg

叠加上近似解 \(u_4(x)\) ：

plot!(x,NSolve(4),color="green",label="u4")

0160.svg

叠加上近似解 \(u_8(x)\) ：

plot!(x,NSolve(8),color="red",label="u8")

0161.svg

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