计算科学 - 为什么我对一阶波动方程的模拟不稳定？ - 吾爱随笔录

为什么我对一阶波动方程的模拟不稳定？

计算科学有限差分稳定波传播

2021-12-13 10:09:34

根据等式

\frac{\partial y}{\partial t} = - a \frac{\partial y}{\partial x}

$\frac{\partial y}{\partial t} = -a\frac{\partial y}{\partial x}$

我在python中模拟了这个。我使用了中心微分，并根据线性方程的 Von-Neumann 稳定性标准确定了步长。 $\Delta t = \frac{\Delta x^2}{2a}$

这就是我得到的：https ://gyazo.com/86024db42f71a6b5cb34eca7eb0d115f

import os
import numpy as np

import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.animation as animation

alpha = 1
seconds = 1000
dx = 1
dt = .1/alpha*dx**2
iter = int(np.round(seconds/dt))
print(iter, dt)


x = np.linspace(0,1000,1000)
# z = np.linspace(0,990, 10, dtype= int)
y = np.zeros([iter,1000])
y[0,:] = 50*np.sin(x[:]*np.pi/1000)
dy_t = np.zeros(1000)


for t in range(0,iter-1):

    for x in range(1, 998):

        dy_t[x] = -alpha*(y[t, x+1]-y[t,x-1])/(2*dx)


    for x in range(1,998):
        y[t+1,x] = y[t,x]+dy_t[x]*dt


# animation
fig, ax = plt.subplots()
line, = ax.plot(y[0,:])
def init():
    line.set_ydata([np.nan] * 1000)
    return line,

def animate(i):
    line.set_ydata(y[i,:])
    return line,

ani = animation.FuncAnimation(
fig, animate, range(0, iter), init_func=init, interval=5, blit=True, save_count=50)
plt.show()

2个回答

空间中心差和时间前向欧拉平流方程的数值解是无条件不稳定的。因此，您看到的行为是预期的。

这是 Gil Strang ( MIT 18.086 ) 的精彩讲座，他在其中讨论了这种方法的不稳定性。他还展示了如何修改简单的中心差分法以产生条件稳定的 Lax-Friedrichs 和 Lax-Wendroff 方法。

假设您真的想在空间中使用中心差分方案并在时间上使用前差分，我将回答这个问题。

考虑问题

u_{t} + a u_{x} = 0, x \in (x_{L}, x_{R}), t \in [0, T]

$u_t + a u_x = 0, \quad x \in (x_L,x_R), \quad t \in [0,T]$ 和一个数字方案

u_{j}^{n + 1} = H_{j} (u^{n}; a, Δ t, h)

$u_j^{n+1} = H_j(u^n; a, \Delta t, h)$ 如果存在，则称该方案是弱稳定的

C = C (T)

$C = C(T)$ 这样

‖ u^{n} ‖ \leq C ‖ u^{0} ‖, n Δ t \leq T

$\| u^n \| \le C \| u^0 \|, \qquad n \Delta t \le T$ 一个等价的定义是如果，那么我们说它是稳定的（强稳定）。FTCS 方案从来都不是强稳定的。它在条件下是弱稳定的这不是一个有用的方案： (1) 时间步太小。(2) 长时间积分无用。

‖ u^{n + 1} ‖ \leq (1 + C Δ t) ‖ u^{n} ‖

$\| u^{n+1} \| \le ( 1 + C \Delta t) \| u^n \|$

C = 1

$C=1$

\frac{u_{j}^{n + 1} - u_{j}^{n}}{Δ t} + \frac{u_{j + 1}^{n} - u_{j - 1}^{n}}{2 h} = 0

$\frac{u_j^{n+1} - u_j^n}{\Delta t} + \frac{u_{j+1}^n - u_{j-1}^n}{2h} = 0$

Δ t \leq {(\frac{h}{a})}^{2}

$\Delta t \le \left( \frac{h}{a} \right)^2$

您的代码中似乎有一些错误，您必须自己修复。我在下面给出了我自己的实现，其中时间步长被选择为用 Python2 运行它；对于 Python3，您必须进行一些更改。

Δ t = C F L {(\frac{h}{a})}^{2}

$\Delta t = CFL \left( \frac{h}{a} \right)^2$

import numpy as np
import matplotlib.pyplot as plt

def smooth(x):
    return np.sin(2*np.pi*x)

def update_ftcs(nu, u):
    unew = np.empty_like(u)
    unew[0] = u[0] + 0.5*nu*(u[-2] - u[1])
    unew[1:-1] = u[1:-1] + 0.5*nu*(u[0:-2] - u[2:])
    unew[-1] = unew[0]
    return unew

def solve():
    scheme     = 'FTCS'
    N          = 100
    cfl        = 0.95
    Tf         = 1.0
    uinit      = smooth
    xmin, xmax = 0.0, 1.0
    a          = 1.0

    h = (xmax - xmin)/N
    dt= cfl * (h / np.abs(a))**2
    nu= a * dt / h

    x = np.linspace(xmin, xmax, N+1)
    u = uinit(x)

    fig = plt.figure()
    ax = fig.add_subplot(111)
    line1, = ax.plot(x, u, 'r')
    line2, = ax.plot(x, u, 'b')
    ax.set_xlabel('x'); ax.set_ylabel('u')
    plt.legend(('Numerical','Exact'))
    plt.title('N='+str(N)+', CFL='+str(cfl)+', Scheme='+scheme)
    plt.draw(); plt.pause(0.1)
    wait = raw_input("Press enter to continue ")

    t, it = 0.0, 0
    while t < Tf:
        u = update_ftcs(nu, u)
        t += dt; it += 1
        if it%100 == 0:
            line1.set_ydata(u)
            line2.set_ydata(uinit(x-a*t))
            plt.draw(); plt.pause(0.1)
    plt.show()

solve()

其它你可能感兴趣的问题

上一篇在 C++ 中从文件存储/读取数据的简单方法下一篇奇异矩阵的不良SVD重建