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

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

使用四阶龙格-库塔法计算Burgers方程

引言
求解过程
完整代码

引言

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

求解过程

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

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

利用中心有限差分法，计算一阶导数 $u_x$
$f'(x_0) \approx \frac{f(x_0+\bigtriangleup x) - f(x_0-\bigtriangleup x)}{2\bigtriangleup 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

利用中心有限差分法，计算二阶导数 $u_{xx}$
$f''(x_0) \approx \frac{f(x_0+\bigtriangleup x) - 2f(x_0) + f(x_0-\bigtriangleup x)}{\bigtriangleup 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

定义四阶龙格-库塔计算公式
对一般微分方程有：
$\begin{cases} y'=f(x,y)\\ y(x_0)=y_0 \end{cases}$
在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{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}$

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_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

计算结果可视化

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方程的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！