使用有限差分法求解二阶微分方程时代码输出异常求助

问题重现 #

尝试用有限差分法求解二阶微分方程，构建矩阵后用np.linalg.solve求解，结果除边界点外全部为0，代码如下：

import numpy as np
from scipy.sparse import csc_matrix

# Define the dimensions of the grid
Lx = 2   # Length of the grid in the x-direction
Nx = 11    # Number of grid points in the x-direction

fa = 2
fb = 0.2

I = np.eye(Nx)

def compute_derivatives(Nx, dx):
    # Compute grid spacing
    x = np.linspace(0, Lx, Nx)
    dx = Lx/(Nx-1)

    # Initialize arrays for first and second order derivatives
    DX = np.zeros(Nx)
    DX2 = np.zeros(Nx)

    # Compute first and second order derivatives using central differences
    for i in range(1, Nx-1):
        DX[i] = (x[i+1] - x[i-1]) / (2*dx)
        DX2[i] = (x[i+1] - 2*x[i] + x[i-1]) / (dx**2)

    return DX, DX2

DX,DX2 = compute_derivatives(Nx,Lx)

A = DX2 + 5*DX + 6*I
b = csc_matrix((Nx, 1)).toarray()
A[[0, Nx-1], :] = 0
A[0,0]=1
A[Nx-1,Nx-1] = 1
b[0] = fa
b[-1] = fb

f = np.linalg.solve(A, np.asarray(np.array(b), dtype=float))
print(f)

核心错误分析 #

1. 差分算子类型错误：混淆了「函数在网格点的导数值」和「差分算子矩阵」。有限差分法中，一阶/二阶导数是通过矩阵-向量乘法作用在未知函数向量f上的，必须构建Nx×Nx的差分矩阵，而非一维数组。
2. 错误计算网格点x的导数：compute_derivatives函数计算的是网格坐标x的一阶/二阶导数，而非未知函数f的差分算子。对于均匀网格，x的一阶导数恒为1，二阶导数恒为0，导致：
   • DX中间元素全为1，DX2中间元素全为0
   • 构建的矩阵A中间行仅对角元素为0 + 5*1 +6*1=11，其余元素为0
   • 结合全0的中间项b，解方程后自然得到f[i]=0

修正后的代码 #

重新构建正确的差分算子矩阵，使用中心差分格式实现一阶/二阶导数：

import numpy as np

# 网格参数
Lx = 2
Nx = 11
dx = Lx / (Nx - 1)
x = np.linspace(0, Lx, Nx)

# 边界条件
fa = 2
fb = 0.2

def build_diff_matrices(Nx, dx):
    # 一阶中心差分矩阵（DX @ f 为f的一阶导数）
    DX = np.zeros((Nx, Nx))
    for i in range(1, Nx-1):
        DX[i, i-1] = -1 / (2*dx)
        DX[i, i+1] = 1 / (2*dx)
    
    # 二阶中心差分矩阵（DX2 @ f 为f的二阶导数）
    DX2 = np.zeros((Nx, Nx))
    for i in range(1, Nx-1):
        DX2[i, i-1] = 1 / dx**2
        DX2[i, i] = -2 / dx**2
        DX2[i, i+1] = 1 / dx**2
    
    return DX, DX2

# 构建差分矩阵
DX, DX2 = build_diff_matrices(Nx, dx)
I = np.eye(Nx)

# 组装方程矩阵 A*f = b
A = DX2 + 5*DX + 6*I
b = np.zeros((Nx, 1))

# 施加边界条件
A[0, :] = 0
A[0, 0] = 1
A[-1, :] = 0
A[-1, -1] = 1
b[0] = fa
b[-1] = fb

# 求解
f = np.linalg.solve(A, b)
print("求解结果：")
print(f)

关键修正点说明 #

• 将DX和DX2改为Nx×Nx的矩阵，对应中心差分的系数：
  • 一阶差分：第i行仅i-1和i+1列有非零值，系数为-1/(2dx)和1/(2dx)
  • 二阶差分：第i行的i-1、i、i+1列有非零值，系数为1/dx²、-2/dx²、1/dx²
• 移除不必要的稀疏矩阵操作，直接用稠密矩阵处理小规模问题
• 提前计算dx，避免函数内重复计算

内容的提问来源于stack exchange，提问作者Chat0924