求解迎风格式的数值解的python代码

以下是一个简单的求解一维迎风格式的数值解的Python代码示例：

import numpy as np

# 定义迎风格式求解函数
def upwind_method(u, dx, dt, nt, c):
    un = np.zeros((nt, len(u)))
    un[0,:] = u
    
    for n in range(1, nt):  
        for i in range(1, len(u)): 
            un[n,i] = u[i] - c*dt/dx*(u[i]-u[i-1]) 
        
        u = un[n].copy()
        
    return un

# 定义初始条件和模拟参数
nx = 81
dx = 4.0/(nx-1)
nt = 25
dt = 0.02
c = 1

u = np.ones(nx)
u[int(0.5/dx):int(1/dx+1)] = 2 # 将一段区间设置为2

# 调用迎风格式函数求解数值解
un = upwind_method(u, dx, dt, nt, c)

# 绘制数值解
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(11,7), dpi=100)
ax = fig.gca(projection='3d')
X, Y = np.meshgrid(np.linspace(0,4,nx), np.linspace(0,nt*dt,nt))
ax.plot_surface(X, Y, un, cmap=plt.cm.viridis)
ax.set_xlabel('Distance')
ax.set_ylabel('Time')
ax.set_zlabel('u')
plt.show()

这个代码示例可以求解一个初始值为1，但在0.5到1之间的区间内为2的一维传输方程的数值解。在这个例子中，我们将空间网格数设置为81，时间步长为0.02，总共求解25个时间步长。通过迎风格式求解函数 `upwind_method()`，我们可以得到数值解 `un`。最后，我们使用 `matplotlib` 绘制数值解的三维图像。