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 方法和正确的空间离散化相差太多，不能一概而论。

应用 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,

等等