Python数值方法建模悬链线：差分格式近似结果偏差问题

问题回顾 #

参考弦振动数值模拟示例，采用数值方法对悬链线建模，但差分格式近似效果差，尤其是初始时间步结果与真实解不符（预期初始步应一致），以下是针对性的问题排查与修正方案：

核心问题与修正措施 #

1. 初始时间步的致命错误 #

• 原代码从t=0开始循环，此时调用t-1会索引到数组末尾的初始0值，直接导致第一步计算完全偏离真实解。必须先初始化t=1的时间步，再从t=1开始循环。
• 初始条件仅设置了t=0的位移，未考虑初始速度（真实解初始速度为0），需用泰勒展开近似t=1的位移：

  # 初始化t=1的位移（初始速度为0）
  for i in range(1, N-1):
      # 从控制方程计算初始加速度
      acc = g * ( (y_ctr[i+1,0]-y_ctr[i-1,0])/(2*dx) + x[i]*(y_ctr[i+1,0]-2*y_ctr[i,0]+y_ctr[i-1,0])/dx**2 )
      y_ctr[i,1] = y_ctr[i,0] + 0.5*dt**2 * acc
  

2. 差分格式推导错误 #

悬链线控制方程展开为：$\frac{\partial^2 y}{\partial t^2} = g\left(\frac{\partial y}{\partial x} + x \frac{\partial^2 y}{\partial x^2}\right)$，原代码的差分格式系数组合错误，正确的中心差分格式应为：

for i in range(1, N-1):
    term1 = (y_ctr[i+1,t] - y_ctr[i-1,t])/(2*dx)  # 一阶空间导数中心差分
    term2 = x[i] * (y_ctr[i+1,t] - 2*y_ctr[i,t] + y_ctr[i-1,t])/dx**2  # 二阶空间导数中心差分
    y_ctr[i,t+1] = 2*y_ctr[i,t] - y_ctr[i,t-1] + dt**2 * g*(term1 + term2)

3. 边界条件缺失 #

悬链线顶端（x=0）固定，位移为0；底端（x=L）自由，空间导数为0，需明确处理：

# 顶端固定
y_ctr[0, :] = 0
# 底端自由，用二阶差分离散
y_ctr[-1, t+1] = 2*y_ctr[-2, t+1] - y_ctr[-3, t+1]

4. 时间尺度错误 #

原代码真实解的时间项用t直接计算，但t是时间步索引，真实时间应为t*dt，否则频率完全偏离：

hc = A * j0(w*bessel_scale) * np.cos(w*t*dt)

5. 数组维度不匹配 #

原代码nT=200但tend=1000，会导致索引越界，需保证tend <= nT-1，或调整nT = tend + 1。

修正后完整代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import j0, jn_zeros

# 可调参数
mode = 1
tend = 100  # 调整为小于nT-1
dt = 0.1

# 常量定义
g = 9.81
L = 10
A = 0.05
x_end = 200
x = np.linspace(0, L, x_end)
bessel_scale = 2 * np.sqrt(x/g)
w = jn_zeros(0, mode)[mode-1] * np.sqrt(g/L) / 2
dx = L / x_end

# 初始化数组
nT = tend + 2  # 确保时间维度足够
N = np.size(x)
y_ctr = np.zeros([N, nT])

# 初始位移（覆盖顶端固定点）
for i in range(N):
    y_ctr[i,0] = A * j0(w*bessel_scale[i])
y_ctr[0,0] = 0  # 顶端固定

# 初始化t=1的位移
for i in range(1, N-1):
    acc = g * ( (y_ctr[i+1,0]-y_ctr[i-1,0])/(2*dx) + x[i]*(y_ctr[i+1,0]-2*y_ctr[i,0]+y_ctr[i-1,0])/dx**2 )
    y_ctr[i,1] = y_ctr[i,0] + 0.5*dt**2 * acc
y_ctr[-1,1] = 2*y_ctr[-2,1] - y_ctr[-3,1]  # 底端自由

# 时间循环
for t in range(1, tend):
    hc = A * j0(w*bessel_scale) * np.cos(w*t*dt)
    # 更新内部点
    for i in range(1, N-1):
        term1 = (y_ctr[i+1,t] - y_ctr[i-1,t])/(2*dx)
        term2 = x[i] * (y_ctr[i+1,t] - 2*y_ctr[i,t] + y_ctr[i-1,t])/dx**2
        y_ctr[i,t+1] = 2*y_ctr[i,t] - y_ctr[i,t-1] + dt**2 * g*(term1 + term2)
    # 更新边界
    y_ctr[-1, t+1] = 2*y_ctr[-2, t+1] - y_ctr[-3, t+1]
    
    # 绘图
    plt.clf()
    plt.xlabel('水平位移')
    plt.xlim((-2*A, 2*A))
    plt.ylabel('竖直位置')
    plt.ylim((0, L))
    plt.title(f'悬链线振动: 模态={mode}')
    plt.plot(y_ctr[:,t], x, label='中心差分法')
    plt.plot(hc, x, 'k', label='真实解')
    plt.legend()
    plt.pause(0.1)

plt.show()

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