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

半撅儿

已于 2024-09-05 13:33:50 修改

阅读量488

点赞数 3

分类专栏： Python学习笔记文章标签： python 龙格-库塔法 Burgers方程

于 2024-09-05 13:28:21 首次发布

本文链接：https://blog.csdn.net/ikhui7/article/details/141927874

版权

Python学习笔记专栏收录该内容

20 篇文章 4 订阅

订阅专栏

使用四阶龙格-库塔法计算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 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_{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 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

定义四阶龙格-库塔计算公式
对一般微分方程有：
$\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 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)

半撅儿

关注

3
点赞
踩
9

收藏

觉得还不错? 一键收藏
0
评论
Python案例 | 使用四阶龙格-库塔法计算Burgers方程

Burgers方程产生于应用数学的各个领域，包括流体力学、非线性声学、气体动力学和交通流。它是一个基本的偏微分方程，可以通过删除压力梯度项从速度场的Navier-Stokes方程导出。对于黏度系数较小的情况()，Burgers方程会导致经典数值方法难以解决的激波形成。
复制链接

扫一扫

专栏目录