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

加权余量法-WRM

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

求解思路：

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

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

Galerkin Method 简介

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

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

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

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

权函数为近似解\(\bar u\)所含待定系数的导数，即：\(w_i=\frac {d \bar u}{dC_i}\)

举例说明

已知微分方程：

\(\left\{\begin{matrix} \frac {d^2u}{dx^2} -u =-x, &0<x<1\\u(0)=0 , &u(1)=0\end{matrix}\right.\)

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

Galerkin Method 具体过程

1. 假设近似解：\(\bar u = C_1x(1-x)\)，显然满足边界条件，\(C_1\)为待定系数
2. 将\(\bar u\)代入原微分方程中，那么余量为：
   \(R=\frac {d^2 \bar u}{dx^2} -\bar u + x=C_1*(-2) - C_1x(1-x) +x\)
3. 选择权函数\(w=\frac{d \bar u}{dC_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)*(C_1x^2+(1-C_1)x-2C_1)dx=0\)，解得\(C_1=\frac{5}{22}\)
5. 所以近似解：\(\bar 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)\)拟合较为良好