计算科学 - 尝试绘制一维波动方程以进行基准测试 - 吾爱随笔录

我正在尝试使用 python 绘制一维波动方程的参考解。

上面的链接说明如下：对于固定在右端并在左端自由并在左端受到锤击的杆 $t = 0$ s 其中锤击引起恒定速度 $V_0$ 从发生 $x<a<L$ ，位移如下：

A_{n} = \frac{2 V_{0} a \sin ((n - \frac{1}{2}) π a / L)}{v π (n - \frac{1}{2})^{2} π a / L}

$A_n = \frac{2V_0 a \sin((n-\frac 1 2)\pi a /L)}{v \pi (n-\frac 1 2)^2 \pi a/L}$

和

ψ (x, t) = \sum_{n = 1, \infty} \cos ((n - \frac{1}{2}) π x / L) \sin ((n - \frac{1}{2}) π t / τ)

$\psi(x,t) = \sum_{n=1,\infty}\cos((n-\frac 1 2)\pi x/L)\sin((n-\frac 1 2)\pi t/\tau)$

在哪里 $\tau = L/c$ , $\psi$ 是位移， $v = \sqrt{E/\rho}$ 和 $a$ 是初始速度的距离 $V_0$ 发生。

看来作者用过 $c$ 在 $\tau = L/c$ 而不是使用 $v$ 为波速？这让我很困惑，但我的直觉告诉我 $c = v = \sqrt{E/\rho}$

对于标准化杆位移 $\hat \psi(x,t) = (c / V_0 a) \psi(x,t)$ 使用前 100 个正常模式（即 $n = 1,100$ ) 并选择 $a/L = 0.1$ , 这就是我们得到的 $t/\tau = 1.28$

现在编码python：

import matplotlib.pyplot as plt
from numpy import *
plt.ion()

n = 100
L = 25
a = 2.5
E = 100
rho = 1
V0 = 3
c = (E/rho)**0.5
x = linspace(0, L, 100)
u = []
ti = 1.28*L/c
for xi in x:
    SUM = 0
    for i in range(1,n+1,1):
        Ai = V0*a*2*sin((i-0.5)*a*pi/L)/(c*pi*(i-0.5)**2*pi*a/L)
        SUM = SUM+Ai*cos((i-0.5)*pi*xi/L)*sin((n-0.5)*pi*ti/(L/c))
    u.append(SUM)


u = [x * (c/V0*a) for x in u]
plt.plot(x/L,u)
plt.show()

产生以下情节：

这实际上与上面的正确情节无关。

我已经检查并重新检查了方程式。我不确定我哪里出错了，希望这不是作者提供的错误？最重要的是最终结果，所以我希望我正确实施它？请帮忙。

编辑

根据下面的答案，这里是工作代码。

import matplotlib.pyplot as plt
from numpy import *
plt.ion()

n = 100
L = 25
a = 2.5
E = 100
rho = 1
V0 = 3
c = (E/rho)**0.5
x = linspace(0, L, 100)
u = []
ti = 1.28*L/c
for xi in x:
    SUM = 0
    for i in range(1,n+1,1):
        Ai = V0*a*2*sin((i-0.5)*a*pi/L)/(c*pi*(i-0.5)**2*pi*a/L)
        SUM = SUM+Ai*cos((i-0.5)*pi*xi/L)*sin((i-0.5)*pi*ti/(L/c))
    u.append(SUM)



u = [x * (c/V0/a) for x in u]
plt.plot(x/L,u)
plt.show()