12步解N-S方程之第5步（2）

在5（1）中我们知道了Burger方程的离散形式及其初始条件应该如何给定。在数学理解上这没有什么难度，但是如何通过程序来给定初始条件却有点麻烦，因为初始条件中含有偏导数项。在这里我们要借助python的sympy符号库来帮助我们完成这些工作。

# -*- coding: utf-8 -*-
#author wubin
#2014/09/26
import numpy as np
import sympy
#tell python we want all of its output to be rendered using latex
#以LATEX格式打印公式
from sympy import init_printing
init_printing(use_latex=True)
#
x,nu,t=sympy.symbols('x nu t')
phi=sympy.exp(-(x-4*t)**2/(t+1)/nu/4)+sympy.exp(-(x-4*t-2*np.pi)**2/(t+1)/nu/4)
phi
#定义偏导数
phiprime=phi.diff(x)
phiprime
#
#用lambdify函数转为可python调用的函数
#
from sympy.utilities.lambdify import lambdify
u=-2*nu*(phiprime/phi)+4
print u
ufunc=lambdify((t,x,nu),u)  #python使用lambdify将u转为t,x,nu三个变量的函数，可以被直接调用
print ufunc(1,4,3)          #测试一下，打印输出
#
import matplotlib.pyplot as plt
###
nx=101
nt=100
dx=2*np.pi/(nx-1)
nu=0.07
dt=nu*dx
#
x=np.linspace(0,2*np.pi,nx)
un=np.empty(nx)
t=0
#
u=np.asarray([ufunc(t,x0,nu) for x0 in x])
# 打印解析解的曲线图
plt.figure(figsize=(11,7),dpi=150)
plt.plot(x,u,marker='o',lw=2)
plt.xlim([0,2*np.pi])
plt.ylim([0,12])
plt.savefig("F:\PYTHON\python_cfd\step5.png",dpi=150)
plt.show()
####数值解
for n in range(nt):
un=u.copy()
for i in range(nx-1):
u[i]=un[i]-un[i]*dt/dx*(un[i]-un[i-1])+nu*dt/dx**2*\
(un[i+1]-2*un[i]+un[i-1])
u[-1]=un[-1]-un[-1]*dt/dx*(un[-1]-un[-2])+nu*dt/dx**2*\
(un[0]-2*un[-1]+un[-2])  #
u_analytical=np.asarray([ufunc(nt*dt,xi,nu)for xi in x])
#
plt.figure(figsize=(11,7), dpi=150)
plt.plot(x,u, marker='o', lw=2, label='Computational')
plt.plot(x, u_analytical, label='Analytical')
plt.xlim([0,2*np.pi])
plt.ylim([0,10])
plt.legend()
plt.savefig("F:\PYTHON\python_cfd\step5_1.png",dpi=150)
plt.show()

图1 解析解曲线（锯齿状曲线）

图2 解析解与数值解对比

至此，一维的例子就到此结束了。在后面的学习中，我们将开始处理2维问题。