Python案例 | 使用四阶龙格-库塔法计算Burgers方程

本文主要是介绍Python案例 | 使用四阶龙格-库塔法计算Burgers方程，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

引言

Burgers方程产生于应用数学的各个领域，包括流体力学、非线性声学、气体动力学和交通流。它是一个基本的偏微分方程，可以通过删除压力梯度项从速度场的Navier-Stokes方程导出。对于黏度系数较小的情况( $\nu =0.01/\pi$)，Burgers方程会导致经典数值方法难以解决的激波形成。在一个空间维度上，带Dirichlet边界条件的Burger方程为：
$$\begin{array}{rlr}& u_t+u{u}_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{array}$$

求解过程

1. 首先，定义一个函数：
   $$f(u,t,dx,\nu )=-u{u}_x+\nu u_{xx}$$

def f(u, t, dx, nu=0.01/np.pi):return -u*dudx(u, dx) + nu*d2udx2(u, dx)

2. 利用中心有限差分法，计算一阶导数$u_x$
   $$f^{\prime }(x_0)\approx \frac{f(x_0+△x)-f(x_0-△x)}{2△x}$$

def dudx(u, dx):"""Approximate the first derivative using the centered finite differenceformula."""first_deriv = np.zeros_like(u)# wrap to compute derivative at endpointsfirst_deriv[0] = (u[1] - u[-1]) / (2*dx)first_deriv[-1] = (u[0] - u[-2]) / (2*dx)# compute du/dx for all the other pointsfirst_deriv[1:-1] = (u[2:] - u[0:-2]) / (2*dx)return first_deriv

3. 利用中心有限差分法，计算二阶导数$u_{xx}$
   $$f^{\prime \prime }(x_0)\approx \frac{f(x_0+△x)-2f(x_0)+f(x_0-△x)}{△x^2}$$

def d2udx2(u, dx):"""Approximate the second derivative using the centered finite differenceformula."""second_deriv = np.zeros_like(u)  # 创建一个新数组second_deriv，其形状和类型与给定数组u相同，但是所有元素都被设置为 0。# wrap to compute second derivative at endpointssecond_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 pointssecond_deriv[1:-1] = (u[2:] - 2*u[1:-1] + u[0:-2]) / (dx**2)return second_deriv

4. 定义四阶龙格-库塔计算公式
   对一般微分方程有：
   $$\{\begin{array}{l}{y}^{\prime }=f(x,y)\\ y(x_0)=y_0\end{array}$$
   在x的取值范围内将其离散为$n$段，定义步长，令第$n$步对应的函数值为$y_n$。于是通过一系列的推导可以得到下一步的$y_{n+1}$值为
   $$y_{n+1}=y_n+\frac{h}{6}(K_1+2K_2+2K_3+K_4)$$
   其中
   $$\{\begin{array}{l}{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+h{K}_3)\end{array}$$

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)

5. Burgers方程计算
   位移初始边界条件：$x_0=-1$，$x_N=1$
   位移离散点个数：$N=512$
   时间初始边界条件：$t_0=0$ ，$t_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 methodu = 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 methodti = t0 + dt*idxU = u[idx, :]for step in range(1000):t = ti + h*stepU = rk4(f, U, t, dx, h)u[idx+1, :] = U

6. 计算结果可视化

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 pltdef 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 differenceformula."""first_deriv = np.zeros_like(u)# wrap to compute derivative at endpointsfirst_deriv[0] = (u[1] - u[-1]) / (2*dx)first_deriv[-1] = (u[0] - u[-2]) / (2*dx)# compute du/dx for all the other pointsfirst_deriv[1:-1] = (u[2:] - u[0:-2]) / (2*dx)return first_derivdef d2udx2(u, dx):"""Approximate the second derivative using the centered finite differenceformula."""second_deriv = np.zeros_like(u)  # 创建一个新数组second_deriv，其形状和类型与给定数组u相同，但是所有元素都被设置为 0。# wrap to compute second derivative at endpointssecond_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 pointssecond_deriv[1:-1] = (u[2:] - 2*u[1:-1] + u[0:-2]) / (dx**2)return second_derivdef 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 pointsdx = (xN - x0) / float(N - 1)  # space between each spatial pointdt = (tK - t0) / float(K)  # space between each temporal pointh = 2e-6  # time step for runge-kutta methodu = np.zeros(shape=(K, N))# u[0, :] = 1 + 0.5*np.exp(-(x**2))  # compute u at initial time stepu[0, :] = -np.sin(np.pi*x)for idx in range(K-1):  # for each temporal point perform runge-kutta methodti = t0 + dt*idxU = u[idx, :]for step in range(1000):t = ti + h*stepU = 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)

这篇关于Python案例 | 使用四阶龙格-库塔法计算Burgers方程的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！

---