python 中平流方程的四阶龙格-库塔规划

Question

%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from math import pi

# wave speed
c = 1
# spatial domain
xmin = 0
xmax = 1
#time domain
m=500; # num of time steps 
tmin=0
T = tmin + np.arange(m+1);
tmax=500

n = 50 # num of grid points

# x grid of n points
X, dx = np.linspace(xmin, xmax, n+1, retstep=True);
X = X[:-1] # remove last point, as u(x=1,t)=u(x=0,t)

# for CFL of 0.1
CFL = 0.3
dt = CFL*dx/c


# initial conditions
def initial_u(x):
    return np.sin(2*pi*x)

# each value of the U array contains the solution for all x values at each timestep
U = np.zeros((m+1,n),dtype=float)
U[0] = u = initial_u(X);




def derivatives(t,u,c,dx):
    uvals = [] # u values for this time step
    for j in range(len(X)):
        if j == 0: # left boundary
            uvals.append((-c/(2*dx))*(u[j+1]-u[n-1]))
        elif j == n-1: # right boundary
            uvals.append((-c/(2*dx))*(u[0]-u[j-1]))
        else:
            uvals.append((-c/(2*dx))*(u[j+1]-u[j-1]))
    return np.asarray(uvals)


# solve for 500 time steps
for k in range(m):
    t = T[k];
    k1 = derivatives(t,u,c,dx)*dt;
    k2 = derivatives(t+0.5*dt,u+0.5*k1,c,dx)*dt;
    k3 = derivatives(t+0.5*dt,u+0.5*k2,c,dx)*dt;
    k4 = derivatives(t+dt,u+k3,c,dx)*dt;
    U[k+1] = u = u + (k1+2*k2+2*k3+k4)/6;

# plot solution
plt.style.use('dark_background')
fig = plt.figure()
ax1 = fig.add_subplot(1,1,1)

line, = ax1.plot(X,U[0],color='cyan')
ax1.grid(True)
ax1.set_ylim([-2,2])
ax1.set_xlim([0,1])
def animate(i):
    line.set_ydata(U[i])
    return line,

我想在 Python 中编写一个平流方程，它是 (∂u/∂t) +c (∂u/∂x) = 0。时间应该用 Runge-kutta 4 阶离散化。空间离散化是二阶有限差分。当我运行我的代码时，我得到了转换为正弦波的直线。但我给出了初始条件正弦波。为什么一开始是直线？我想让正弦波向前移动。你知道如何让正弦波向前移动吗？我感谢您的帮助。提前致谢！

Answer 1

虽然从表面上看你的计算步骤与 RK4 方法有关，但它们与 RK4 方法和正确的 space 离散化相差太多，无法一一列举。

应用 ODE 积分方法的传统方法是使用一个函数 derivatives(t, state, params)，然后应用它来计算 Euler 步长或 RK4 步长。在你的情况下它将是

def derivatives(t,u,c,dx):
    du = np.zeros(len(u));
    p = c/(2*dx);
    du[0] = p*(u[1]-u[-1]);
    du[1:-1] = p*(u[2:]-u[:-2]);
    du[-1] = p*(u[0]-u[-2]);
    return du;

那你可以做

X, dx = np.linspace(xmin, xmax, n+1, retstep=True);
X = X[:-1] # remove last point, as u(x=1,t)=u(x=0,t)

m=500; # number of time steps
T = tmin + np.arange(m+1);

U = np.zeros((m+1,n),dtype=float)
U[0] = u = initial_u(X);
for k in range(m):
    t = T[k];
    k1 = derivatives(t,u,c,dx)*dt;
    k2 = derivatives(t+0.5*dt,u+0.5*k1,c,dx)*dt;
    k3 = derivatives(t+0.5*dt,u+0.5*k2,c,dx)*dt;
    k4 = derivatives(t+dt,u+k3,c,dx)*dt;
    U[k+1] = u = u + (k1+2*k2+2*k3+k4)/6;

这使用dt作为时间步长中计算的主要变量，然后从tmin构造等差数列，步长为dt。其他方法也是可能的，但必须使 tmax 和时间步数兼容。

到目前为止的计算现在应该是成功的，可以在动画中使用了。以我的理解，你不会在每一帧中产生一个新的情节，你只绘制一次图形，然后改变线数据

# animate the time data
line, = ax1.plot(X,U[0],color='cyan')
ax1.grid(True)
ax1.set_ylim([-2,2])
ax1.set_xlim([0,1])

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

等等

python 中平流方程的四阶龙格-库塔规划

Programming of 4th order Runge-Kutta in advection equation in python

python

numerical-methods

runge-kutta