python - 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 四阶离散化。空间离散化是二阶有限差分。当我运行代码时，我得到的直线转变为正弦波。但我给出了正弦波作为初始条件。为什么它一开始是直线？我想让正弦波继续前进。你知道如何让正弦波向前发展吗？我感谢您的帮助。提前致谢!

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 作为时间步长中的主要变量进行计算，然后使用 dt 步长从 tmin 构造算术序列。其他方法也是可能的，但必须使 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,

等等