一维对流扩散方程的数值解法
对流-扩散方程是守恒定律控制方程的一种模型方程，它既是能量方程的表示形式，同时也可以认为是把压力梯度项隐含到了源项中去的动量方程的代表。因此，以对流-扩散方程为例，来研究数值求解偏微分方程的相容性、收敛性和稳定性具有代表性的意义。
1 数学模型
本作业从最简单的模型方程，即一维、稳态、无源项的对流扩散方程出发，方程如下： 22, 02f f f
U D x t x x
∂∂∂+=≤≤∂∂∂ (1)
初始条件 (),0sin(2)f x t A kx π==
(2)
解析解
()()()224,sin 2Dk t
f x t e
A k x Ut ππ-=-
(3)
式中，1,0.05,0.5,1U D A k ====
函数(3)描述的是一个衰减波的图像，如图1所示
t=0 t=0.5 t=1
图1 函数()()()224,sin 2Dk t
f x t e
k x Ut ππ-=- 的图像（U=1，D=0.05，k=1）
2 数值解法
2.1 数值误差分析
在网格点(),i n 上差分方程的数值解n
i f 偏离该点上相应的偏微分方程的精确解
(),f i n 的值，称为网格节点上的数值误差。
当取定网格节点数21N =时，观察差分方程的解与微分方程的解在不同时间步长下的趋近程度，其中时间步长分别取值0.05,0.025,0.0125,0.0005t ∆=。
（a ）21,0.05N t =∆= （b ）21,0.025N t =∆=
（c ）21,0.0125N t =∆= （d ）201,0.0005N t =∆=
图2 数值误差随步长的变化情况
从图2的(a)~(d)可以定性的看出，数值误差与步长的大小有关。在满足稳定性条件的前提下，数值误差随着时间步长的减小而减小，同时，图（d ）表示增大网格的分辨率也有助于减小网格误差。
为了对数值误差有一个定量的认识，接下来取定时间步长为0.0005t ∆=，分别算出
11,21,41,61,81,101,121,161N =时，指标E =1所示。
表1 不同网格节点数下指标E 的值
由表1可以看出，指标E 的值随着网格节点数的增加呈先大幅减小后平缓变化的趋势，该趋势示于图3中。
（a ） （b ）
图3 不同网格节点数下指标E 的值
在图3的（a ）和（b ）中，随着网格节点数的增加，曲线最后都趋于水平，即此时再增加网格节点数，对于提高数值求解的精确性作用不大。 2.2 截断误差分析
对于方程
()5f x x =
(4)
其一阶导数和二阶导数的精确解分别为4
5f x x
∂=∂和23220f x x ∂=∂。
一阶导数和二阶导数的中心差分表示式分别为
211()2i i f f f
o h x h
+--∂=+∂和221122
2()i i i f f f f
o h x h
+--+∂=+∂，2()o h 为截断误差。 当空间步长分别取值为0.2,0.1,0.05,0.025,0.0125h =时，节点()1i x =取处一阶导数和二阶导数的精确解与离散解的差值见表2。
表2 截断误差随空间步长的变化
图示如下
图4 截断误差随空间步长的变化
由图4可以看出，截断误差随着步长的减小而减小，当步长趋近于0时，截断误差趋近于0。
3 结论
（1）在满足稳定性条件的前提下，数值误差随着时间步长的减小而减小。而对于空间步长的减小，数值误差呈先大幅减小后变平缓的趋势，即存在一个网格无关解的网格节点数。
（2）对于导数在空间上的离散，二阶精度的中心差分格式能较好的满足数值求解精确性的要求，当选取合适的步长时，其精确性也将显著提高。
%% =============================================================== clear
clc
format long
%% ===============================================================
U=1;D=0.05;k=1;A=0.5;
T=0.5;dt=0.05;m=T/dt;
L=2;N=21;h=L/(N-1);
%% ===============================================================
E0=0; E=0; %误差估计
%% ===============================================================
x=linspace(0,L,N);y=[1:N];
for i=1:N
f0(i)=A*sin(2*pi*k*x(i)); %初始条件
end
%% ===============================================================
for is=1:m,is
%% =============================================================== for i=2:N-1
f(i)=f0(i)-(U*dt/(2*h))*(f0(i+1)-f0(i-1))+(D*dt/(h^2))*(f0(i+1)-2 *f0(i)+f0(i-1));
end
%f(1)=f0(1)-(U*dt/(2*h))*(f0(2)-f0(N-1))+(D*dt/(h^2))*(f0(2)-2*f0 (1)+f0(N-1)); %边界条件
%f(N)=f0(N)-(U*dt/(2*h))*(f0(2)-f0(N-1))+(D*dt/(h^2))*(f0(2)-2*f0 (N)+f0(N-1)); %边界条件
f(1)=exp(-4*D*k^2*pi^2*dt*is)*A*sin(2*pi*k*(0-U*dt*is)); %边界条件 f(N)=exp(-4*D*k^2*pi^2*dt*is)*A*sin(2*pi*k*(L-U*dt*is)); %±边界条件 f0=f;
end
%% ===============================================================
for i=1:N
fexact(i)=exp(-4*D*k^2*pi^2*T)*A*sin(2*pi*k*(x(i)-U*T));
E0=(f(i)-fexact(i))^2+E0;
end
E=h*sqrt(E0)
%% ===============================================================
f1=exp(-4*D*k^2*pi^2*T)*A*sin(2*pi*k*(x-U*T));
plot(y,f1,'r',y,f,'b','linewidt',2),legend('Exact','Numerical'),axis( [1 N -2 2]),grid on
%plot(y,f1,'r','linewidt',4),legend('Exact'),axis([1 N -2 2]),grid on %hold on
%plot(y,f,'b','linewidt',2),legend('Numerical')
