一维Burgers方程数值解法

一维Burgers方程

一维burgers方程为：

由于等式右边可以进行积分：

利用F = u**2，则方程为：

假设u初始为阶跃函数：

数值解法采用MacCormack格式：

但是这一解法，有失真的性质，后面具体介绍。

所以根据这一格式，可以直接数值求解，并利用matplotlib画出动态的数值图形，具体代码如下：

# -*- coding: utf-8 -*-
"""
Created on Tue Jan 20 14:32:23 2015
1D burges equation
@author: myjiayan
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation

def u_initial():
    first = np.ones(40)
    second = np.zeros(41)
    result = np.array(list(first)+list(second))
    return result

def computeF(u):
    F = 0.5*u**2
    return F


def maccormack(u,nt,dt,dx):
    un = np.zeros((nt,len(u)))
    un[0] = u.copy()
    ustar = u.copy()
    
    for t in xrange(1,nt):
        F = computeF(u)
        ustar[:-1] = u[:-1] - (F[1:]-F[:-1])*dt/dx
        ustar[-1] = 0.0
        Fstar = computeF(ustar)
        un[t,1:] = 0.5*(u[1:]+ustar[1:]-(Fstar[1:]-Fstar[:-1])*dt/dx)
        un[t,0] = 1.0
        u = un[t].copy()
    
    return un

if __name__ == '__main__':
    nx = 81
    nt = 70
    dx = 4.0/(nx-1)

    def animate(data):
        x = np.linspace(0,4,nx)
        y = data
        line.set_data(x,y)
        return line,

    u = u_initial()
    sigma = .5
    dt = sigma*dx

    un = maccormack(u,nt,dt,dx)
    fig = plt.figure();
    ax = plt.axes(xlim=(0,4),ylim=(-.5,2));
    line, = ax.plot([],[],lw=2);
    anim = animation.FuncAnimation(fig, animate, frames=un, interval=50)
    plt.show()

直接将代码保存为burgers.py文件，打开terminal：

$python burgers.py

数值结果就以动态的形式表现出来了。

很明显，数值结果失真了，数值结果中u的值竟然比1大，如何改进MacCormack格式呢？

我们在预测步上加一个人工项，新的格式为：

选取Ephsilon为0.5后改进格式得到的数值结果为：

比1大的值已经消失。