python求解常微分方程之Galerkin method：权函数

权函数的概念

在伽辽金法中，权函数是一组满足边界条件的函数，它们用于将微分方程转化为一组代数方程。权函数的选择对近似解的精度和稳定性有重要影响。

对结果的影响

权函数的选择直接影响近似解的精度。如果权函数能够很好地逼近真实解，那么近似解的精度就会很高。此外，权函数的正交性可以简化求解过程，提高计算效率。

怎么选择以提高精度

1. 满足边界条件：权函数必须满足与微分方程相同的边界条件。
2. 正交性：选择正交函数作为权函数可以简化求解过程。
3. 完备性：权函数集合应该能够覆盖整个解空间。
4. 光滑性：权函数应该足够光滑，以便于计算其导数。

举例说明

对于给定的微分方程 $\frac{d^{2}u}{d{x}^2}-u=-x$，边界条件为 $u(0)=0$ 和 $u(1)=0$，我们可以选择以下权函数：

1. 单个权函数：
   $$w_1(x)=x(1-x)$$
   这是最基本的选择，适用于初步近似。

2. 多个权函数：
   为了提高精度，我们可以选择四个权函数，它们构成一组基，可以更全面地逼近解：
   $$w_1(x)=x(1-x),w_2(x)=x^2(1-x),w_3(x)=x^3(1-x),w_4(x)=x^4(1-x)$$
   这些权函数都是关于 $x$ 的多项式，它们满足边界条件 $w(0)=0$ 和 $w(1)=0$，并且它们在区间$[0,1]$ 上是线性独立的。

设置多个权函数

当我们增加权函数的数量时，我们实际上是在增加近似解的自由度。每个权函数对应一个待定系数，我们需要为每个权函数设置一个积分方程，并求解这些方程以找到待定系数。对于四个权函数，我们设置四个积分方程：

$${\int }_0^1{w}_i(x)\left(\frac{d^{2}u}{d{x}^2}-u+x\right)dx=0,i=1,2,3,4$$

将近似解$$u(x)=C_{1}x(1-x)+C_2{x}^2(1-x)+C_3{x}^3(1-x)+C_4{x}^4(1-x)$$代入上述积分方程，并求解得到的线性方程组以找到 $C_1,C_2,C_3,C_4$的值。

问题描述

使用伽辽金法求解二阶微分方程：
$$\frac{d^{2}u}{d{x}^2}-u=-x,0<x<1$$
边界条件为：
$$u(0)=0,u(1)=0$$
权函数设为：
$$w_1=x(1-x),w_2=x^2(1-x)$$
近似解设为：
$$u(x)=C_{1}x(1-x)+C_2{x}^2(1-x)$$

求解步骤(二个权函数)

1. 计算近似解的导数：
   $$u(x)=C_{1}x(1-x)+C_2{x}^2(1-x)$$
   $$u^{\prime }(x)=C_1(1-2x)+C_2(2x-3x^2)$$
   $$u^{\prime \prime }(x)=-2C_1+2C_2(1-4x)$$

2. 将近似解及其导数代入微分方程：
   $$\frac{d^{2}u}{d{x}^2}-u=(-2C_1+2C_2-4C_{2}x+2C_2{x}^2)-(C_{1}x+C_1{x}^2-C_2{x}^2-C_2{x}^3)+x$$
   $$=-3C_1+2C_2-4C_{2}x+2C_2{x}^2-C_{1}x-C_1{x}^2+C_2{x}^2+C_2{x}^3+x$$
   $$=(-3C_1+2C_2)+(-4C_2-C_1)x+(2C_2-C_1)x^2+C_2{x}^3+x$$

3. 定义残差函数：
   $$R(x)=(-3C_1+2C_2)+(-4C_2-C_1)x+(2C_2-C_1)x^2+C_2{x}^3+x$$

4. 最小化残值：
   $${\int }_0^1{w}_1(x)R(x)dx=0$$
   $${\int }_0^1{w}_2(x)R(x)dx=0$$

5. 求解线性方程组：
   通过求解上述两个积分方程，我们可以得到关于 $C_1$ 和$C_2$ 的线性方程组。这可以通过数值积分和线性代数方法来完成。

6. 得到近似解：
   将求得的 $C_1$ 和 $C_2$ 代入近似解的表达式中，得到最终的近似解。

最终答案

$$\boxed{u(x)=C_{1}x(1-x)+C_2{x}^2(1-x)}$$

其中 $C_1$ 和$C_2$是通过求解线性方程组得到的系数。

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt

# 定义符号变量
x, C1, C2 = sp.symbols('x C1 C2')

# 定义近似解（一个权函数）
u_single = C1 * x * (1 - x)

# 计算微分方程的残差（一个权函数）
R_single = sp.diff(u_single, x, 2) - u_single + x

# 定义权函数（一个权函数）
w_single = x * (1 - x)

# 计算积分（一个权函数）
I_single = sp.integrate(w_single * R_single, (x, 0, 1))

# 解方程求C1（一个权函数）
res_single = sp.solve(I_single, C1)
C1_single_value = res_single[0]

# 定义近似解（两个权函数）
u_double = C1 * x * (1 - x) + C2 * x**2 * (1 - x)

# 计算微分方程的残差（两个权函数）
R_double = sp.diff(u_double, x, 2) - u_double + x

# 定义权函数（两个权函数）
w1 = x * (1 - x)
w2 = x**2 * (1 - x)

# 计算积分（两个权函数）
I1 = sp.integrate(w1 * R_double, (x, 0, 1))
I2 = sp.integrate(w2 * R_double, (x, 0, 1))

# 解方程求C1和C2（两个权函数）
res_double = sp.solve((I1, I2), (C1, C2))
C1_double_value, C2_double_value = res_double[C1], res_double[C2]

# 计算解析解
def u_analytical(x):
    return x - (np.exp(x) - np.exp(-x)) / (np.exp(1) - np.exp(-1))

# 将解析解和近似解转换为可计算的形式
x_vals = np.linspace(0, 1, 100)
u_analytical_vals = u_analytical(x_vals)
u_single_approx_vals = [float(u_single.subs(C1, C1_single_value).subs(x, xv)) for xv in x_vals]
u_double_approx_vals = [float(u_double.subs({C1: C1_double_value, C2: C2_double_value}).subs(x, xv)) for xv in x_vals]

# 计算残差
R_single_vals = [float(R_single.subs(C1, C1_single_value).subs(x, xv)) for xv in x_vals]
R_double_vals = [float(R_double.subs({C1: C1_double_value, C2: C2_double_value}).subs(x, xv)) for xv in x_vals]

# 绘制函数值对比图
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.plot(x_vals, u_analytical_vals, label='Analytical Solution', color='blue')
plt.plot(x_vals, u_single_approx_vals, label='Single Weight Function Approximation', linestyle='--', color='red')
plt.plot(x_vals, u_double_approx_vals, label='Double Weight Functions Approximation', linestyle='-.', color='green')
plt.title('Function Value Comparison')
plt.xlabel('x')
plt.ylabel('u(x)')
plt.legend()

# 绘制残差对比图
plt.subplot(1, 2, 2)
plt.plot(x_vals, R_single_vals, label='Residuals Single Weight', color='red')
plt.plot(x_vals, R_double_vals, label='Residuals Double Weight', color='green')
plt.title('Residual Comparison')
plt.xlabel('x')
plt.ylabel('R(x)')
plt.legend()

plt.tight_layout()
plt.show()