微分方程数值解法3-Galerkin Method 伽辽金法

加权余量法-WRM

微分方程数值求解方法中，有一类为加权余余量法（weighted residual methods，WRM）

求解思路：

• 选近似解，要求满足边界条件
• 求余量
• 选权函数，并确定权函数的系数。根据权函数选择方法的不同，WRM有许多分支

【参考】：刘金原，《计算物理学》，2012.6

Galerkin Method 简介

微分方程数值求解方法中，有一类为加权余余量法（weighted residual methods，WRM），Galerkin Method 伽辽金法属于加权余量法的一种

余量：是指用方程近似解代入方程后，得到近似解与真实解的偏差，即残差或误差

近似解通常采用带有任意系数的一组线性无关函数表示

加权余量法的思想就是如何确定这组系数，使余量最小

权函数为近似解$\overline{u}$所含待定系数的导数，即：$w_i=\frac{d\overline{u}}{d{C}_i}$

举例说明

已知微分方程：

$\{\begin{array}{cc}\frac{d^{2}u}{d{x}^2}-u=-x,& 0<x<1\\ u(0)=0,& u(1)=0\end{array}$

那么方程解析解：$u=x-\frac{e^x-e^{-x}}{e-e^{-1}}$

Galerkin Method 具体过程

1. 假设近似解：$\overline{u}=C_{1}x(1-x)$，显然满足边界条件，$C_1$为待定系数
2. 将$\overline{u}$代入原微分方程中，那么余量为：
   $R=\frac{d^2\overline{u}}{d{x}^2}-\overline{u}+x=C_1\ast (-2)-C_{1}x(1-x)+x$
3. 选择权函数$w=\frac{d\overline{u}}{d{C}_i}$使得余量$R$在区间(0,1)内的加权积分为0，即
   $I={\int }_0^{1}wRdx={\int }_0^{1}w(C_1(-2-x(1-x))+x)dx=0$
   这里的$w=x(1-x)$
4. 得到$I={\int }_0^{1}x(1-x)\ast (C_1{x}^2+(1-C_1)x-2C_1)dx=0$，解得$C_1=\frac{5}{22}$
5. 所以近似解：$\overline{u}=\frac{5}{22}x(1-x)$

由上可知，使用近似解为插值多项式（幂函数），可以很好使得权函数、余量以及加权余量均为多项式，容易得到加权函数的积分表达式，进而求解待定系数

代码实现

import sympy
from sympy.abc import *
from sympy import symbols

x, C1 = symbols("x,C1")
u = 'C1*x*(1-x)' # 字符串
u = sympy.sympify(u) # 字符串转换为计算式，假设原方程的解为该方程u，因为原方程初值u(0)=0,u(1)=0

w = sympy.diff(u, C1) # u对C1求导，w=du/dC1，w为权函数，这里w=x*(1 - x)
print(w)

R = sympy.diff(u, x, 2)- u + x # u对x二次求导，R=d2u/dx2-u+x，即需要求解的方程为R=0
print(R)

# 选择权函数w使得余量在区间【0，1】内加权积分为0，即integrate【0，1】wR dx=0,均为幂函数
I = sympy.integrate(w*R,(x,0,1))
print('______',I,'______')

res = sympy.solve(I,C1) # res为C1的值，积分式I=0时CI的值
print(res) # 那么原方程R=0的近似解u=5/22*x*(1-x)

# R=0（R=d2u/dx2-u+x，u(0)=0,u(1)=0）的精确解u=x-(e^x-e^{-x})/(e-e^{-1});绘制图形对比

import numpy as np
import matplotlib.pyplot as plt
plt.rc('font',family='Alibaba PuHuiTi') # 运行matplotlib_font.py，查询内置的字体，含中文的字库即可

x = np.linspace(0,1,40) # x = np.linspace(-1,2,40)
u=x-(np.exp(x)-np.exp(-x))/(np.exp(1)-np.exp(-1))  # 解析解
u_error = 5/22*x*(1-x)  # 伽辽金法近似解

# fig = plt.figure()
# ax = fig.add_subplot(1,1,1)

plt.plot(x,u,'g',label='解析解') # 解析解-绿色
plt.plot(x,u_error,'y',label='近似解') # 近似解-黄色
plt.legend()
plt.show()

解析解与近似解对比

放大绘图区间，区间$(0,1)$拟合较为良好