常微分方程（十二）边值问题
Chen Kai BOSS
  2019-06-03 10:45:00 2019-06-03 10:45      5.7k 字  28 分钟  

初值问题告诉我们起点的全部信息，而边值问题则在两端各给出部分信息——就像知道出发点和目的地，但不知道中间的完整路径。这类问题在工程中极为常见：梁的弯曲、热传导的稳态、薛定谔方程的本征态问题。本章将深入探讨边值问题的理论和数值求解方法。

什么是边值问题？

从初值问题到边值问题

打靶法示意图

回顾初值问题（ IVP）：

所有条件都在同一点 给出。

边值问题（ BVP）：

条件在不同点给出：左边界和右边界。

为什么边值问题更难？

对于初值问题，存在唯一性定理（ Picard 定理）保证解的存在唯一性。但边值问题可能：

  1. 无解：条件不兼容
  2. 有无穷多解：条件不足以确定
  3. 有唯一解：这是我们希望的情况

示例 1：考虑 - BVP1：, - 解：，任意都满足，无穷多解

  • BVP2：, - 需要，矛盾，无解

  • BVP3：, - 唯一解

边界条件的类型

第一类（ Dirichlet）：给定函数值

第二类（ Neumann）：给定导数值

第三类（ Robin/混合）：给定函数值和导数的线性组合

线性边值问题

标准形式

有限差分网格

二阶线性 BVP：

解的结构

线性 BVP 的解可以表示为：

其中 是特解，, 是齐次方程的基本解。

叠加原理方法

方法：

  1. 解两个初值问题：
    • 问题 A：，， - 问题 B：，，$y'(a) = 1解为y = y_A + c y_B由边界条件y(b) = $ 确定：

（要求）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp, solve_bvp
from scipy.optimize import brentq, fsolve

def superposition_method(p, q, r, a, b, alpha, beta):
    """
    线性 BVP 的叠加原理方法
    y'' + p(x)y' + q(x)y = r(x)
    y(a) = alpha, y(b) = beta
    """
    # 转换为一阶系统
    def system_A(x, Y):
        y, dy = Y
        return [dy, r(x) - p(x)*dy - q(x)*y]
    
    def system_B(x, Y):
        y, dy = Y
        return [dy, -p(x)*dy - q(x)*y]
    
    # 解问题 A
    sol_A = solve_ivp(system_A, [a, b], [alpha, 0], dense_output=True)
    
    # 解问题 B
    sol_B = solve_ivp(system_B, [a, b], [0, 1], dense_output=True)
    
    # 计算常数 c
    yA_b = sol_A.sol(b)[0]
    yB_b = sol_B.sol(b)[0]
    
    if abs(yB_b) < 1e-10:
        raise ValueError("方法失效： y_B(b) ≈ 0")
    
    c = (beta - yA_b) / yB_b
    
    # 构造解
    def solution(x):
        return sol_A.sol(x)[0] + c * sol_B.sol(x)[0]
    
    return solution

# 示例： y'' - y = -x, y(0) = 0, y(1) = 1
p = lambda x: 0
q = lambda x: -1
r = lambda x: -x
a, b = 0, 1
alpha, beta = 0, 1

sol = superposition_method(p, q, r, a, b, alpha, beta)

# 精确解： y = x + A*sinh(x) + B*cosh(x)
# y(0) = 0 => B = 0
# y(1) = 1 => 1 + A*sinh(1) = 1 => A = 0
# 实际上精确解是 y = x (可以验证)
# 等等，让我重新算： y'' = y - x
# 特解 y_p = x（验证： 0 = x - x ✓）
# 齐次解 y_h = A*e^x + B*e^{-x}
# y(0) = A + B = 0
# y(1) = 1 + A*e + B*e^{-1} = 1
# 从第一个： B = -A
# 代入第二个： A*e - A*e^{-1} = 0 => A = 0
# 所以 y = x

x_plot = np.linspace(a, b, 100)
y_num = np.array([sol(x) for x in x_plot])
y_exact = x_plot  # 精确解是 y = x

plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.plot(x_plot, y_num, 'b-', linewidth=2, label='数值解')
plt.plot(x_plot, y_exact, 'r--', linewidth=2, label='精确解')
plt.xlabel('x')
plt.ylabel('y')
plt.title('叠加原理方法')
plt.legend()
plt.grid(True, alpha=0.3)

plt.subplot(122)
plt.plot(x_plot, abs(y_num - y_exact), 'k-')
plt.xlabel('x')
plt.ylabel('|误差|')
plt.title('误差')
plt.grid(True, alpha=0.3)
plt.yscale('log')

plt.tight_layout()
plt.show()

打靶法

基本思想

格林函数的物理意义

打靶法（ Shooting Method）将 BVP 转化为 IVP：

  1. 猜测缺失的初始条件
  2. 解初值问题到右边界
  3. 检查是否满足右边界条件
  4. 如果不满足，调整猜测值，重复

就像射击：调整枪的角度（初始斜率）直到击中目标（满足右边界条件）。

非线性 BVP 的打靶法

考虑：

我们不知道。

算法：

  1. 定义残差函数 - 其中 是初值问题, 的解

  2. 求解（用 Newton 法或二分法）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
def shooting_method(f, a, b, alpha, beta, s_guess, tol=1e-8, max_iter=50):
    """
    非线性 BVP 的打靶法
    y'' = f(x, y, y')
    y(a) = alpha, y(b) = beta
    
    返回: x 数组, y 数组
    """
    def residual(s):
        """给定初始斜率 s，计算边界残差"""
        def system(x, Y):
            y, dy = Y
            return [dy, f(x, y, dy)]
        
        sol = solve_ivp(system, [a, b], [alpha, s], dense_output=True)
        return sol.sol(b)[0] - beta
    
    # 用 Brent 方法或割线法求解
    try:
        # 先尝试二分法（需要找到包围根的区间）
        s_low, s_high = -10, 10
        while residual(s_low) * residual(s_high) > 0:
            s_low *= 2
            s_high *= 2
        s_opt = brentq(residual, s_low, s_high, xtol=tol)
    except:
        # 如果二分法失败，用 Newton 法
        s_opt = fsolve(residual, s_guess, full_output=False)[0]
    
    # 用最优 s 计算最终解
    def system(x, Y):
        y, dy = Y
        return [dy, f(x, y, dy)]
    
    sol = solve_ivp(system, [a, b], [alpha, s_opt], dense_output=True)
    
    x_out = np.linspace(a, b, 100)
    y_out = sol.sol(x_out)[0]
    
    return x_out, y_out, s_opt

# 示例：非线性 BVP
# y'' = -e^y, y(0) = 0, y(1) = 0
# 这是 Bratu 问题的一个形式
f = lambda x, y, dy: -np.exp(y)
a, b = 0, 1
alpha, beta = 0, 0

x, y, s = shooting_method(f, a, b, alpha, beta, s_guess=0)

plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.plot(x, y, 'b-', linewidth=2)
plt.xlabel('x')
plt.ylabel('y')
plt.title(f'打靶法解：初始斜率 s = {s:.6f}')
plt.grid(True, alpha=0.3)

# 验证边界条件
print(f"y(0) = {y[0]:.6f} (应为 {alpha})")
print(f"y(1) = {y[-1]:.6f} (应为 {beta})")

# 打靶过程的可视化
s_values = np.linspace(-2, 2, 9)
plt.subplot(122)
for s_try in s_values:
    def system(x, Y):
        return [Y[1], -np.exp(Y[0])]
    sol = solve_ivp(system, [a, b], [alpha, s_try], dense_output=True)
    x_plot = np.linspace(a, b, 50)
    y_plot = sol.sol(x_plot)[0]
    plt.plot(x_plot, y_plot, '-', alpha=0.5, label=f's={s_try:.1f}')

plt.axhline(y=beta, color='r', linestyle='--', label='目标 y(1)=0')
plt.xlabel('x')
plt.ylabel('y')
plt.title('打靶过程：不同初始斜率的轨迹')
plt.legend(fontsize=8)
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('shooting_method.png', dpi=150, bbox_inches='tight')
plt.show()

打靶法的敏感性

打靶法有时会很敏感：解对初始斜率的依赖可能是指数级的。考虑：

通解是。在 处， 而。

微小的初始误差会被放大约 倍！

多重打靶法

思想：将区间分成多段，每段独立打靶，然后在内部节点匹配。

这减少了误差的指数放大。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
def multiple_shooting(f, a, b, alpha, beta, n_segments=4):
    """
    多重打靶法
    将[a,b]分成 n_segments 段
    """
    # 节点
    nodes = np.linspace(a, b, n_segments + 1)
    h = nodes[1] - nodes[0]
    
    # 初始猜测：线性插值
    y_guess = np.linspace(alpha, beta, n_segments + 1)
    s_guess = np.ones(n_segments) * (beta - alpha) / (b - a)
    
    def residuals(params):
        """
        params = [y_1, ..., y_{n-1}, s_0, s_1, ..., s_{n-1}]
        y_0 = alpha, y_n = beta 是固定的
        """
        y_interior = params[:n_segments-1]
        s_all = params[n_segments-1:]
        
        y_nodes = np.concatenate([[alpha], y_interior, [beta]])
        
        res = []
        
        for i in range(n_segments):
            x_left, x_right = nodes[i], nodes[i+1]
            y_left, s_left = y_nodes[i], s_all[i]
            
            def system(x, Y):
                return [Y[1], f(x, Y[0], Y[1])]
            
            sol = solve_ivp(system, [x_left, x_right], [y_left, s_left])
            y_right_computed = sol.y[0, -1]
            s_right_computed = sol.y[1, -1]
            
            # 连续性条件
            res.append(y_right_computed - y_nodes[i+1])
            if i < n_segments - 1:
                res.append(s_right_computed - s_all[i+1])  # 斜率连续
        
        return res
    
    # 初始参数
    initial_params = np.concatenate([y_guess[1:-1], s_guess])
    
    # 求解
    solution = fsolve(residuals, initial_params)
    
    # 重构解
    y_interior = solution[:n_segments-1]
    s_all = solution[n_segments-1:]
    y_nodes = np.concatenate([[alpha], y_interior, [beta]])
    
    # 在每段上求解并合并
    x_full, y_full = [], []
    for i in range(n_segments):
        x_left, x_right = nodes[i], nodes[i+1]
        
        def system(x, Y):
            return [Y[1], f(x, Y[0], Y[1])]
        
        sol = solve_ivp(system, [x_left, x_right], [y_nodes[i], s_all[i]],
                       dense_output=True)
        
        x_seg = np.linspace(x_left, x_right, 25)
        y_seg = sol.sol(x_seg)[0]
        
        x_full.extend(x_seg[:-1] if i < n_segments-1 else x_seg)
        y_full.extend(y_seg[:-1] if i < n_segments-1 else y_seg)
    
    return np.array(x_full), np.array(y_full)

有限差分法

基本思想

Sturm-Liouville问题

有限差分法直接离散化微分方程，将 BVP 转化为代数方程组。

离散化导数：（中心差分）

线性 BVP 的有限差分

考虑：

在节点（，）上离散化：

整理得：

矩阵形式

这是一个三对角线性系统：

其中： - - - -（需要修正边界项）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
def finite_difference_linear(p, q, r, a, b, alpha, beta, N):
    """
    线性 BVP 的有限差分法
    y'' + p(x)y' + q(x)y = r(x)
    y(a) = alpha, y(b) = beta
    """
    h = (b - a) / N
    x = np.linspace(a, b, N + 1)
    
    # 内部节点
    x_interior = x[1:-1]
    n = len(x_interior)
    
    # 构造三对角矩阵
    A = np.zeros((n, n))
    b_vec = np.zeros(n)
    
    for i in range(n):
        xi = x_interior[i]
        pi, qi, ri = p(xi), q(xi), r(xi)
        
        # 对角线
        A[i, i] = -2 + h**2 * qi
        
        # 下对角线
        if i > 0:
            A[i, i-1] = 1 - h * pi / 2
        
        # 上对角线
        if i < n - 1:
            A[i, i+1] = 1 + h * pi / 2
        
        # 右端项
        b_vec[i] = h**2 * ri
        
        # 边界修正
        if i == 0:
            b_vec[i] -= (1 - h * pi / 2) * alpha
        if i == n - 1:
            b_vec[i] -= (1 + h * pi / 2) * beta
    
    # 求解线性系统
    y_interior = np.linalg.solve(A, b_vec)
    
    # 添加边界值
    y = np.concatenate([[alpha], y_interior, [beta]])
    
    return x, y

# 示例： y'' - y = -x, y(0) = 0, y(1) = 1
p = lambda x: 0
q = lambda x: -1
r = lambda x: -x
a, b = 0, 1
alpha, beta = 0, 1

# 不同网格密度
N_values = [10, 20, 40, 80]

plt.figure(figsize=(14, 5))

plt.subplot(121)
for N in N_values:
    x, y = finite_difference_linear(p, q, r, a, b, alpha, beta, N)
    plt.plot(x, y, 'o-', markersize=3, label=f'N = {N}')

# 精确解
# y'' = y - x 的精确解
# y = x 满足 y'' = 0, y - x = -x + x = 0 ??? 
# 让我重新验证： y = x, y' = 1, y'' = 0
# y'' - y = 0 - x = -x ✓
x_exact = np.linspace(a, b, 200)
y_exact = x_exact
plt.plot(x_exact, y_exact, 'k-', linewidth=2, label='精确解')

plt.xlabel('x')
plt.ylabel('y')
plt.title('有限差分法')
plt.legend()
plt.grid(True, alpha=0.3)

# 收敛性分析
plt.subplot(122)
errors = []
for N in [10, 20, 40, 80, 160, 320]:
    x, y = finite_difference_linear(p, q, r, a, b, alpha, beta, N)
    y_exact_at_nodes = x
    error = np.max(np.abs(y - y_exact_at_nodes))
    errors.append(error)

h_values = (b - a) / np.array([10, 20, 40, 80, 160, 320])
plt.loglog(h_values, errors, 'bo-', label='最大误差')
plt.loglog(h_values, h_values**2, 'r--', label='O(h ²) 参考')
plt.xlabel('步长 h')
plt.ylabel('最大误差')
plt.title('收敛阶验证')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('finite_difference.png', dpi=150, bbox_inches='tight')
plt.show()

非线性 BVP 的有限差分

对于非线性方程，离散化后得到非线性代数方程组，需要用 Newton 法求解。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
def finite_difference_nonlinear(f, df_dy, df_ddy, a, b, alpha, beta, N, 
                                 tol=1e-10, max_iter=100):
    """
    非线性 BVP 的有限差分法（ Newton 迭代）
    y'' = f(x, y, y')
    y(a) = alpha, y(b) = beta
    """
    h = (b - a) / N
    x = np.linspace(a, b, N + 1)
    n = N - 1  # 内部节点数
    
    # 初始猜测：线性插值
    y = np.linspace(alpha, beta, N + 1)
    
    for iteration in range(max_iter):
        # 构造残差向量
        F = np.zeros(n)
        
        for i in range(1, N):
            xi = x[i]
            yi = y[i]
            dyi = (y[i+1] - y[i-1]) / (2*h)
            d2yi = (y[i+1] - 2*y[i] + y[i-1]) / h**2
            
            F[i-1] = d2yi - f(xi, yi, dyi)
        
        # 检查收敛
        if np.max(np.abs(F)) < tol:
            break
        
        # 构造 Jacobian 矩阵
        J = np.zeros((n, n))
        
        for i in range(1, N):
            xi = x[i]
            yi = y[i]
            dyi = (y[i+1] - y[i-1]) / (2*h)
            
            # J[i-1, j-1] = ∂ F[i-1]/∂ y[j]
            # F[i-1] = (y[i+1] - 2y[i] + y[i-1])/h ² - f(x[i], y[i], dy[i])
            
            # 对 y[i-1]的偏导
            if i > 1:
                J[i-1, i-2] = 1/h**2 + df_ddy(xi, yi, dyi) / (2*h)
            
            # 对 y[i]的偏导
            J[i-1, i-1] = -2/h**2 - df_dy(xi, yi, dyi)
            
            # 对 y[i+1]的偏导
            if i < N - 1:
                J[i-1, i] = 1/h**2 - df_ddy(xi, yi, dyi) / (2*h)
        
        # Newton 更新
        dy = np.linalg.solve(J, -F)
        y[1:-1] += dy
    
    return x, y

# 示例： y'' = -e^y, y(0) = 0, y(1) = 0 (Bratu 问题)
f = lambda x, y, dy: -np.exp(y)
df_dy = lambda x, y, dy: -np.exp(y)
df_ddy = lambda x, y, dy: 0

x_fd, y_fd = finite_difference_nonlinear(f, df_dy, df_ddy, 0, 1, 0, 0, N=50)

# 与打靶法比较
x_sh, y_sh, _ = shooting_method(f, 0, 1, 0, 0, s_guess=0)

plt.figure(figsize=(10, 5))
plt.plot(x_fd, y_fd, 'b.-', label='有限差分')
plt.plot(x_sh, y_sh, 'r--', label='打靶法')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Bratu 问题：两种方法比较')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

配置法与样条方法

Galerkin 方法

特征函数的正交性

思想：将解展开为基函数的线性组合，然后通过残差正交性确定系数。

设，其中 是基函数。

Galerkin 条件：残差与每个基函数正交：

有限元方法简介

有限元方法是配置法的重要特例：

  1. 将区间分成单元
  2. 在每个单元上用简单的多项式（如线性或二次）
  3. 要求函数在单元边界连续
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
def fem_linear_bvp(p, q, r, a, b, alpha, beta, N):
    """
    简单有限元方法（线性元）
    -y'' + p(x)y' + q(x)y = r(x)
    """
    h = (b - a) / N
    x = np.linspace(a, b, N + 1)
    n = N - 1
    
    # 刚度矩阵和质量矩阵的组装
    K = np.zeros((n, n))  # 刚度矩阵（来自 y''项）
    M = np.zeros((n, n))  # 质量矩阵（来自 y 项）
    C = np.zeros((n, n))  # 对流矩阵（来自 y'项）
    F = np.zeros(n)       # 载荷向量
    
    for i in range(n):
        xi = x[i + 1]  # 第 i 个内部节点
        
        # 对角线元素（来自左右两个单元的贡献）
        K[i, i] = 2 / h
        M[i, i] = 2 * h / 3 * q(xi)
        
        # 非对角线元素
        if i > 0:
            K[i, i-1] = -1 / h
            M[i, i-1] = h / 6 * q(x[i])
        if i < n - 1:
            K[i, i+1] = -1 / h
            M[i, i+1] = h / 6 * q(x[i+2])
        
        # 载荷向量（简单近似）
        F[i] = h * r(xi)
    
    # 边界修正
    F[0] += alpha / h
    F[-1] += beta / h
    
    # 组装总矩阵
    A = K + M
    
    # 求解
    y_interior = np.linalg.solve(A, F)
    y = np.concatenate([[alpha], y_interior, [beta]])
    
    return x, y

特征值问题

Sturm-Liouville 问题

标准形式：

这里 是特征值， 是特征函数。

有限差分求特征值

离散化后得到广义特征值问题。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
def sturm_liouville_eigenvalues(p, q, w, a, b, N, num_eig=5):
    """
    Sturm-Liouville 特征值问题
    -(p(x)y')' + q(x)y = λ w(x)y
    y(a) = y(b) = 0
    """
    h = (b - a) / N
    x = np.linspace(a, b, N + 1)
    n = N - 1
    
    # 构造矩阵 A（来自-(py')' + qy）
    A = np.zeros((n, n))
    B = np.zeros((n, n))  # 质量矩阵（来自 wy）
    
    for i in range(n):
        xi = x[i + 1]
        
        # 主对角线
        p_plus = p(xi + h/2)
        p_minus = p(xi - h/2)
        A[i, i] = (p_plus + p_minus) / h**2 + q(xi)
        B[i, i] = w(xi)
        
        # 次对角线
        if i > 0:
            A[i, i-1] = -p_minus / h**2
        if i < n - 1:
            A[i, i+1] = -p_plus / h**2
    
    # 求解广义特征值问题
    from scipy.linalg import eigh
    eigenvalues, eigenvectors = eigh(A, B)
    
    # 重构特征函数（添加边界值 0）
    eigenfunctions = []
    for j in range(num_eig):
        ef = np.concatenate([[0], eigenvectors[:, j], [0]])
        # 归一化
        ef = ef / np.sqrt(np.sum(ef**2 * h))
        eigenfunctions.append(ef)
    
    return x, eigenvalues[:num_eig], eigenfunctions

# 示例：简单的 Sturm-Liouville 问题
# -y'' = λ y, y(0) = y(π) = 0
# 精确特征值：λ_n = n ², 特征函数： sin(nx)
p = lambda x: 1
q = lambda x: 0
w = lambda x: 1
a, b = 0, np.pi
N = 100

x, eigenvalues, eigenfunctions = sturm_liouville_eigenvalues(p, q, w, a, b, N, num_eig=6)

# 与精确值比较
exact_eigenvalues = np.array([1, 4, 9, 16, 25, 36])

print("Sturm-Liouville 特征值：")
print("数值解\t\t 精确解\t\t 相对误差")
for i, (num, exact) in enumerate(zip(eigenvalues, exact_eigenvalues)):
    rel_err = abs(num - exact) / exact
    print(f"{num:.6f}\t{exact}\t\t{rel_err:.2e}")

# 绘制特征函数
fig, axes = plt.subplots(2, 3, figsize=(14, 8))
axes = axes.flatten()

for i in range(6):
    axes[i].plot(x, eigenfunctions[i], 'b-', linewidth=2, label='数值')
    exact_ef = np.sin((i+1) * x)
    exact_ef = exact_ef / np.sqrt(np.sum(exact_ef**2) * (b-a)/len(x))
    axes[i].plot(x, exact_ef, 'r--', linewidth=2, label='精确')
    axes[i].set_title(f'λ_{i+1} = {eigenvalues[i]:.4f}')
    axes[i].set_xlabel('x')
    axes[i].legend()
    axes[i].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('eigenvalues.png', dpi=150, bbox_inches='tight')
plt.show()

实际应用

梁的弯曲

简支梁的挠度 满足：

其中 是杨氏模量， 是截面惯性矩， 是载荷分布。

边界条件（简支）： 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
def beam_deflection(EI, q, L, N):
    """
    简支梁的弯曲
    EI*y'''' = q(x)
    y(0) = y(L) = 0, y''(0) = y''(L) = 0
    """
    h = L / N
    x = np.linspace(0, L, N + 1)
    n = N + 1
    
    # 四阶导数的有限差分：
    # y''''(x_i) ≈ (y_{i-2} - 4y_{i-1} + 6y_i - 4y_{i+1} + y_{i+2}) / h^4
    
    # 构造方程组
    # 对于四阶问题，需要更复杂的边界处理
    # 这里用简化的方法
    
    A = np.zeros((n, n))
    b = np.zeros(n)
    
    for i in range(n):
        if i == 0 or i == n - 1:
            # y = 0 边界
            A[i, i] = 1
            b[i] = 0
        elif i == 1:
            # y'' = 0 在 x=0 附近
            # y''(x_1) ≈ (y_0 - 2y_1 + y_2)/h ² = 0
            # 但 y_0 = 0，所以 -2y_1 + y_2 = 0
            A[i, 0] = 1
            A[i, 1] = -2
            A[i, 2] = 1
            b[i] = 0
        elif i == n - 2:
            # y'' = 0 在 x=L 附近
            A[i, n-3] = 1
            A[i, n-2] = -2
            A[i, n-1] = 1
            b[i] = 0
        else:
            # 内部点： EI*y'''' = q
            A[i, i-2] = EI / h**4
            A[i, i-1] = -4 * EI / h**4
            A[i, i] = 6 * EI / h**4
            A[i, i+1] = -4 * EI / h**4
            A[i, i+2] = EI / h**4
            b[i] = q(x[i])
    
    y = np.linalg.solve(A, b)
    return x, y

# 均布载荷
EI = 1
L = 1
q = lambda x: 1  # 均布载荷

x, y = beam_deflection(EI, q, L, N=50)

# 精确解： y = q/(24EI) * x*(L-x)*(L ²+x(L-x)) = q*x*(x ³-2Lx ²+L ³)/(24EI)
# 简化形式： y = (q*L ⁴)/(24EI) * (x/L)*(1-x/L)*((x/L)² - (x/L) + 1)
y_exact = (1/(24*EI)) * x * (L - x) * (L**2 - x*(L-x))

plt.figure(figsize=(10, 5))
plt.plot(x, -y*1000, 'b.-', label='数值解')
plt.plot(x, -y_exact*1000, 'r--', label='精确解')
plt.xlabel('位置 x')
plt.ylabel('挠度 (× 10 ⁻³)')
plt.title('简支梁在均布载荷下的弯曲')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

热传导稳态

一维稳态热传导：

边界条件可以是： - 温度给定（ Dirichlet） - 热流给定（ Neumann） - 对流边界（ Robin）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
def heat_conduction_steady(k, Q, L, T_left, T_right, N):
    """
    一维稳态热传导
    -k*T'' = Q(x)
    T(0) = T_left, T(L) = T_right
    """
    h = L / N
    x = np.linspace(0, L, N + 1)
    
    # 有限差分
    A = np.zeros((N + 1, N + 1))
    b = np.zeros(N + 1)
    
    # 边界条件
    A[0, 0] = 1
    b[0] = T_left
    A[-1, -1] = 1
    b[-1] = T_right
    
    # 内部节点
    for i in range(1, N):
        A[i, i-1] = -k / h**2
        A[i, i] = 2 * k / h**2
        A[i, i+1] = -k / h**2
        b[i] = Q(x[i])
    
    T = np.linalg.solve(A, b)
    return x, T

# 带内热源的杆
k = 1
L = 1
Q = lambda x: 10 * np.sin(np.pi * x)  # 正弦分布的热源
T_left, T_right = 0, 0

x, T = heat_conduction_steady(k, Q, L, T_left, T_right, N=50)

# 精确解： T'' = -10*sin(π x)/k
# T = 10/(k*π²) * sin(π x)
T_exact = 10 / (k * np.pi**2) * np.sin(np.pi * x)

plt.figure(figsize=(10, 5))
plt.plot(x, T, 'b.-', label='数值解')
plt.plot(x, T_exact, 'r--', label='精确解')
plt.xlabel('位置 x')
plt.ylabel('温度 T')
plt.title('带正弦热源的一维热传导')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

量子力学：薛定谔方程

一维定态薛定谔方程：

这是一个特征值问题！

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
def schrodinger_1d(V, x_min, x_max, N, num_states=5, hbar=1, m=1):
    """
    一维定态薛定谔方程
    返回能级和波函数
    """
    h = (x_max - x_min) / N
    x = np.linspace(x_min, x_max, N + 1)
    n = N - 1  # 内部点
    
    # Hamiltonian 矩阵
    H = np.zeros((n, n))
    
    kinetic_coeff = hbar**2 / (2 * m * h**2)
    
    for i in range(n):
        xi = x[i + 1]
        
        # 动能项（-ℏ²/(2m) d ²/dx ²）
        H[i, i] = 2 * kinetic_coeff + V(xi)
        if i > 0:
            H[i, i-1] = -kinetic_coeff
        if i < n - 1:
            H[i, i+1] = -kinetic_coeff
    
    # 求解特征值问题
    energies, wavefunctions = np.linalg.eigh(H)
    
    # 重构波函数（添加边界值 0）
    psi = []
    for j in range(num_states):
        wf = np.concatenate([[0], wavefunctions[:, j], [0]])
        # 归一化
        wf = wf / np.sqrt(np.sum(wf**2) * h)
        psi.append(wf)
    
    return x, energies[:num_states], psi

# 示例 1：无限深势阱
# V = 0 内部, V = ∞ 边界
# 精确能级： E_n = n ²π²ℏ²/(2mL ²)
V_well = lambda x: 0
L = 1
x, E_well, psi_well = schrodinger_1d(V_well, 0, L, N=200, num_states=4)

# 示例 2：谐振子势
# V = (1/2)kx ² = (1/2)m ω² x ²
omega = 1
V_harmonic = lambda x: 0.5 * omega**2 * x**2
x_harm, E_harm, psi_harm = schrodinger_1d(V_harmonic, -5, 5, N=200, num_states=4)
# 精确能级： E_n = (n + 1/2)ℏω

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 无限深势阱
axes[0, 0].set_title('无限深势阱的能级和波函数')
for i in range(4):
    # 波函数偏移显示
    axes[0, 0].plot(x, psi_well[i]**2 + E_well[i], label=f'n={i+1}, E={E_well[i]:.4f}')
    axes[0, 0].axhline(y=E_well[i], linestyle='--', alpha=0.5)
axes[0, 0].set_xlabel('x')
axes[0, 0].set_ylabel('|ψ|² + E')
axes[0, 0].legend()

# 无限深势阱的能级比较
axes[0, 1].set_title('无限深势阱：数值 vs 精确能级')
n_levels = np.arange(1, 5)
E_exact_well = n_levels**2 * np.pi**2 / (2 * L**2)
axes[0, 1].bar(n_levels - 0.2, E_well, 0.4, label='数值')
axes[0, 1].bar(n_levels + 0.2, E_exact_well, 0.4, label='精确')
axes[0, 1].set_xlabel('量子数 n')
axes[0, 1].set_ylabel('能量 E')
axes[0, 1].legend()

# 谐振子
axes[1, 0].set_title('量子谐振子的能级和波函数')
for i in range(4):
    axes[1, 0].fill_between(x_harm, psi_harm[i]**2 * 3 + E_harm[i], E_harm[i], alpha=0.5)
    axes[1, 0].axhline(y=E_harm[i], linestyle='--', alpha=0.5)
# 势能
axes[1, 0].plot(x_harm, V_harmonic(x_harm), 'k-', linewidth=2, label='V(x)')
axes[1, 0].set_xlabel('x')
axes[1, 0].set_ylabel('E')
axes[1, 0].set_ylim([0, 5])
axes[1, 0].legend()

# 谐振子能级比较
axes[1, 1].set_title('量子谐振子：数值 vs 精确能级')
E_exact_harm = (np.arange(4) + 0.5) * omega
axes[1, 1].bar(np.arange(4) - 0.2, E_harm, 0.4, label='数值')
axes[1, 1].bar(np.arange(4) + 0.2, E_exact_harm, 0.4, label='精确 (n+1/2)ℏω')
axes[1, 1].set_xlabel('量子数 n')
axes[1, 1].set_ylabel('能量 E')
axes[1, 1].legend()

plt.tight_layout()
plt.savefig('schrodinger.png', dpi=150, bbox_inches='tight')
plt.show()

print("\n 谐振子能级：")
print("数值\t\t 精确\t\t 误差")
for i in range(4):
    print(f"{E_harm[i]:.6f}\t{E_exact_harm[i]:.6f}\t{abs(E_harm[i]-E_exact_harm[i]):.2e}")

使用 SciPy 求解 BVP

SciPy 提供了solve_bvp函数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
from scipy.integrate import solve_bvp

# 示例： y'' = -e^y, y(0) = 0, y(1) = 0
def fun(x, y):
    """将二阶方程转化为一阶系统"""
    return np.vstack((y[1], -np.exp(y[0])))

def bc(ya, yb):
    """边界条件： ya 是 x=a 处的值， yb 是 x=b 处的值"""
    return np.array([ya[0] - 0,   # y(0) = 0
                     yb[0] - 0])  # y(1) = 0

# 初始网格和猜测
x = np.linspace(0, 1, 10)
y = np.zeros((2, x.size))  # 初始猜测

# 求解
sol = solve_bvp(fun, bc, x, y)

if sol.success:
    x_plot = np.linspace(0, 1, 100)
    y_plot = sol.sol(x_plot)[0]
    
    plt.figure(figsize=(10, 5))
    plt.plot(x_plot, y_plot, 'b-', linewidth=2)
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('solve_bvp 求解 Bratu 问题')
    plt.grid(True, alpha=0.3)
    plt.show()
else:
    print("求解失败:", sol.message)

# 更复杂的例子：带参数的问题
# 求解非线性边值问题，其中解依赖于参数
def fun_param(x, y, p):
    """p 是未知参数"""
    return np.vstack((y[1], -p[0] * np.exp(y[0])))

def bc_param(ya, yb, p):
    return np.array([ya[0], yb[0]])

# 这种情况下需要额外的条件来确定参数

练习题

基础题

练习 1：用叠加原理方法求解：

练习 2：用打靶法求解：

练习 3：用有限差分法（ N=20）求解：

并与精确解比较。

练习 4：考虑特征值问题：

用有限差分法求前 5 个特征值。

练习 5：验证对于,,，当 时无解。

中级题

练习 6：用多重打靶法（ 4 段）求解：

比较与单次打靶的稳定性。

练习 7：实现带 Neumann 边界条件的有限差分法：

练习 8：用配置法（ Chebyshev 节点）求解：

练习 9：求解 Sturm-Liouville 问题：

并绘制前 4 个特征函数。

练习 10：实现简支梁在集中载荷下的弯曲分析。

高级题

练习 11：考虑非线性特征值问题（ Bratu 问题）：

对不同的 值，研究解的存在性和多重性。绘制分支图。

练习 12：用有限元方法求解变系数问题：

练习 13：研究数值方法的超收敛现象：对于某些问题，在特定点上的误差比全局误差更小。

练习 14：实现后验误差估计：根据数值解估计真实误差，用于自适应网格细化。

练习 15：考虑奇异边值问题：

讨论在 处的奇异性处理。

总结

本章我们学习了边值问题的各种数值方法：

方法 优点 缺点 适用场景
叠加原理 简单，利用现有 IVP 求解器 仅限线性问题 线性 BVP
打靶法 灵活，可处理非线性 可能不稳定 一般 BVP
多重打靶 更稳定 实现复杂 敏感的 BVP
有限差分 直接，易于理解 需要解线性系统 规则区域
有限元 处理复杂几何 理论复杂 工程问题

实际建议：

  1. 先用scipy.integrate.solve_bvp尝试
  2. 如果失败，检查初始猜测
  3. 对于敏感问题，考虑多重打靶或有限差分
  4. 对于特征值问题，有限差分通常足够

下一章，我们将进入偏微分方程的世界——当未知函数依赖于多个变量时会发生什么？

  • 本文标题：常微分方程（十二）边值问题
  • 本文作者：Chen Kai
  • 创建时间：2019-06-03 10:45:00
  • 本文链接：https://www.chenk.top/%E5%B8%B8%E5%BE%AE%E5%88%86%E6%96%B9%E7%A8%8B%EF%BC%88%E5%8D%81%E4%BA%8C%EF%BC%89%E8%BE%B9%E5%80%BC%E9%97%AE%E9%A2%98/
  • 版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
 评论