引言
Burgers方程产生于应用数学的各个领域,包括流体力学、非线性声学、气体动力学和交通流。它是一个基本的偏微分方程,可以通过删除压力梯度项从速度场的Navier-Stokes方程导出。对于黏度系数较小的情况(
ν
=
0.01
/
π
\nu = 0.01/ \pi
ν=0.01/π),Burgers方程会导致经典数值方法难以解决的激波形成。在一个空间维度上,带Dirichlet边界条件的Burger方程为:
u
t
+
u
u
x
−
ν
u
x
x
=
0
,
x
∈
[
−
1
,
1
]
,
t
∈
[
0
,
1
]
u
(
0
,
x
)
=
−
s
i
n
(
π
x
)
u
(
t
,
−
1
)
=
u
(
t
,
1
)
=
0
\begin{align*} & u_t + uu_x - \nu u_{xx} = 0 , x \in [-1,1], t \in [0,1] & \\ & u(0,x) = -sin(\pi x) & \\ & u(t,-1) = u(t,1) = 0 & \end{align*}
ut+uux−νuxx=0,x∈[−1,1],t∈[0,1]u(0,x)=−sin(πx)u(t,−1)=u(t,1)=0
求解过程
- 首先,定义一个函数:
f ( u , t , d x , ν ) = − u u x + ν u x x f(u,t,dx,\nu)= -uu_x + \nu u_{xx} f(u,t,dx,ν)=−uux+νuxx
def f(u, t, dx, nu=0.01/np.pi):
return -u*dudx(u, dx) + nu*d2udx2(u, dx)
- 利用中心有限差分法,计算一阶导数
u
x
u_x
ux
f ′ ( x 0 ) ≈ f ( x 0 + △ x ) − f ( x 0 − △ x ) 2 △ x f'(x_0) \approx \frac{f(x_0+\bigtriangleup x) - f(x_0-\bigtriangleup x)}{2\bigtriangleup x} f′(x0)≈2△xf(x0+△x)−f(x0−△x)
def dudx(u, dx):
"""
Approximate the first derivative using the centered finite difference
formula.
"""
first_deriv = np.zeros_like(u)
# wrap to compute derivative at endpoints
first_deriv[0] = (u[1] - u[-1]) / (2*dx)
first_deriv[-1] = (u[0] - u[-2]) / (2*dx)
# compute du/dx for all the other points
first_deriv[1:-1] = (u[2:] - u[0:-2]) / (2*dx)
return first_deriv
- 利用中心有限差分法,计算二阶导数
u
x
x
u_{xx}
uxx
f ′ ′ ( x 0 ) ≈ f ( x 0 + △ x ) − 2 f ( x 0 ) + f ( x 0 − △ x ) △ x 2 f''(x_0) \approx \frac{f(x_0+\bigtriangleup x) - 2f(x_0) + f(x_0-\bigtriangleup x)}{\bigtriangleup x^2} f′′(x0)≈△x2f(x0+△x)−2f(x0)+f(x0−△x)
def d2udx2(u, dx):
"""
Approximate the second derivative using the centered finite difference
formula.
"""
second_deriv = np.zeros_like(u) # 创建一个新数组second_deriv,其形状和类型与给定数组u相同,但是所有元素都被设置为 0。
# wrap to compute second derivative at endpoints
second_deriv[0] = (u[1] - 2*u[0] + u[-1]) / (dx**2)
second_deriv[-1] = (u[0] - 2*u[-1] + u[-2]) / (dx**2)
# compute d2u/dx2 for all the other points
second_deriv[1:-1] = (u[2:] - 2*u[1:-1] + u[0:-2]) / (dx**2)
return second_deriv
- 定义四阶龙格-库塔计算公式
对一般微分方程有:
{ y ′ = f ( x , y ) y ( x 0 ) = y 0 \begin{cases} y'=f(x,y)\\ y(x_0)=y_0 \end{cases} {y′=f(x,y)y(x0)=y0
在x的取值范围内将其离散为 n n n段,定义步长,令第 n n n步对应的函数值为 y n y_n yn。于是通过一系列的推导可以得到下一步的 y n + 1 y_{n+1} yn+1值为
y n + 1 = y n + h 6 ( K 1 + 2 K 2 + 2 K 3 + K 4 ) y_{n+1}=y_n+\frac{h}{6} (K_1+2K_2+2K_3+K_4) yn+1=yn+6h(K1+2K2+2K3+K4)
其中
{ K 1 = f ( x n , y n ) K 2 = f ( x n + h 2 , y n + h 2 K 1 ) K 3 = f ( x n + h 2 , y n + h 2 K 2 ) K 4 = f ( x n + h , y n + h K 3 ) \begin{cases} K_1=f(x_n, y_n) \\ K_2=f(x_n+\frac{h}{2}, y_n+\frac{h}{2}K_1) \\ K_3=f(x_n+\frac{h}{2},y_n+\frac{h}{2}K_2) \\ K_4=f(x_n+h,y_n+hK_3) \end{cases} ⎩ ⎨ ⎧K1=f(xn,yn)K2=f(xn+2h,yn+2hK1)K3=f(xn+2h,yn+2hK2)K4=f(xn+h,yn+hK3)
def rk4(f, u, t, dx, h):
"""
Fourth-order Runge-Kutta method for computing u at the next time step.
"""
k1 = f(u, t, dx)
k2 = f(u + 0.5*h*k1, t + 0.5*h, dx)
k3 = f(u + 0.5*h*k2, t + 0.5*h, dx)
k4 = f(u + h*k3, t + h, dx)
return u + (h/6)*(k1 + 2*k2 + 2*k3 + k4)
- Burgers方程计算
位移初始边界条件: x 0 = − 1 x_0=-1 x0=−1, x N = 1 x_N=1 xN=1
位移离散点个数: N = 512 N=512 N=512
时间初始边界条件: t 0 = 0 t_0=0 t0=0 , t K = 500 t_K=500 tK=500
时间离散点个数: K = 500 K=500 K=500
x = np.linspace(x0, xN, N) # evenly spaced spatial points
dx = (xN - x0) / float(N - 1) # space between each spatial point
dt = (tK - t0) / float(K) # space between each temporal point
h = 2e-6 # time step for runge-kutta method
u = np.zeros(shape=(K, N))
# u[0, :] = 1 + 0.5*np.exp(-(x**2)) # compute u at initial time step
u[0, :] = -np.sin(np.pi*x)
for idx in range(K-1): # for each temporal point perform runge-kutta method
ti = t0 + dt*idx
U = u[idx, :]
for step in range(1000):
t = ti + h*step
U = rk4(f, U, t, dx, h)
u[idx+1, :] = U
- 计算结果可视化
plt.imshow(u.T, interpolation='nearest', cmap='rainbow',
extent=[t0, tK, x0, xN], origin='lower', aspect='auto')
plt.xlabel('t')
plt.ylabel('x')
plt.colorbar()
plt.show()
完整代码
""" Solving the Burgers' Equation using a 4th order Runge-Kutta method """
import numpy as np
import matplotlib.pyplot as plt
def rk4(f, u, t, dx, h):
"""
Fourth-order Runge-Kutta method for computing u at the next time step.
"""
k1 = f(u, t, dx)
k2 = f(u + 0.5*h*k1, t + 0.5*h, dx)
k3 = f(u + 0.5*h*k2, t + 0.5*h, dx)
k4 = f(u + h*k3, t + h, dx)
return u + (h/6)*(k1 + 2*k2 + 2*k3 + k4)
def dudx(u, dx):
"""
Approximate the first derivative using the centered finite difference
formula.
"""
first_deriv = np.zeros_like(u)
# wrap to compute derivative at endpoints
first_deriv[0] = (u[1] - u[-1]) / (2*dx)
first_deriv[-1] = (u[0] - u[-2]) / (2*dx)
# compute du/dx for all the other points
first_deriv[1:-1] = (u[2:] - u[0:-2]) / (2*dx)
return first_deriv
def d2udx2(u, dx):
"""
Approximate the second derivative using the centered finite difference
formula.
"""
second_deriv = np.zeros_like(u) # 创建一个新数组second_deriv,其形状和类型与给定数组u相同,但是所有元素都被设置为 0。
# wrap to compute second derivative at endpoints
second_deriv[0] = (u[1] - 2*u[0] + u[-1]) / (dx**2)
second_deriv[-1] = (u[0] - 2*u[-1] + u[-2]) / (dx**2)
# compute d2u/dx2 for all the other points
second_deriv[1:-1] = (u[2:] - 2*u[1:-1] + u[0:-2]) / (dx**2)
return second_deriv
def f(u, t, dx, nu=0.01/np.pi):
return -u*dudx(u, dx) + nu*d2udx2(u, dx)
def make_square_axis(ax):
ax.set_aspect(1 / ax.get_data_ratio())
def burgers(x0, xN, N, t0, tK, K):
x = np.linspace(x0, xN, N) # evenly spaced spatial points
dx = (xN - x0) / float(N - 1) # space between each spatial point
dt = (tK - t0) / float(K) # space between each temporal point
h = 2e-6 # time step for runge-kutta method
u = np.zeros(shape=(K, N))
# u[0, :] = 1 + 0.5*np.exp(-(x**2)) # compute u at initial time step
u[0, :] = -np.sin(np.pi*x)
for idx in range(K-1): # for each temporal point perform runge-kutta method
ti = t0 + dt*idx
U = u[idx, :]
for step in range(1000):
t = ti + h*step
U = rk4(f, U, t, dx, h)
u[idx+1, :] = U
# plt.imshow(u, extent=[x0, xN, t0, tK])
plt.imshow(u.T, interpolation='nearest', cmap='rainbow',
extent=[t0, tK, x0, xN], origin='lower', aspect='auto')
plt.xlabel('t')
plt.ylabel('x')
plt.colorbar()
plt.show()
if __name__ == '__main__':
# burgers(-10, 10, 1024, 0, 50, 500)
burgers(-1, 1, 512, 0, 1, 500)