对流扩散方程

徐州工程学院

课程设计报告

课程名称偏微分方程数值解

课题名称对流扩散方程

的迎风格式的推导和求解专业信息与计算科学

班级10信计3

姓名学号

指导教师杨扬

2013年 5 月23 日

一、实验目的：

进一步巩固理论学习的结果，学习双曲型对流扩散方程的迎风格式的构造

方法，以及稳定的条件。从而进一步了解差分求解偏微分方程的一些基本概念，掌握数值求解偏微分方程的基本过程。在此基础上考虑如何使用Matlab 的软件进行上机实现，并针对具体的题目给出相应的数值计算结果。

二、实验题目：

⎪⎩

⎪

⎨⎧-=-==<<<<+=+);2/1exp(),1();exp(),0();2/exp()0,(10,10,11t t u t t u x x u t x f u b u a u xx x t 其中a1=1,b1=2,

)2/exp(),(t x t x f --=。

用迎风格式求解双曲型对流扩散方程,观差分解对真解的敛散性()2/exp(t x u -=

三、实验原理：

1、用迎风格式求解双曲型对流扩散方程，迎风格式为：

)

01(21

1

)01(2112

1

1112

1

11

1<++-=-+->++-=-+--+++-+-+a f h

u u u b h

u u a u u a f h

u u u b h

u u a u u n

j n

j n j n j n

j

n j n

j

n j n j n

j n j n j n

j n j n j

n j τ

τ

若令,/*1,/*12h b h a r

τμτ==

则迎风格式可整理为：

>

<<++-+-+=><>++++--=-+++-+2)01()()21(1)01()()21(111111a f u u r u r u a f u u r u r u n j

n j n j n j n j n

j n j n j n j n j τμμμτμμμ2、稳定条件：

）

（）

（01),*11*2/(01),*11*2/(2

2<-≤>+≤a h a b h a h a b h ττ（*） 四、数值实验的过程、相关程序及结果：

本次的实验题目所给出的边界条件是第一边界条件，直接利用所给的边界条件，我们可以给出界点处以及第0层的函数值，根据a1的正负性，使用相应的<1>或者<2>式，求出其他层的函数值。误差转化成图的形式，并输出最大值。 针对三种不同的输入对应输出结果 ：

A: a1=1;b1=2;a=1;b=1;h=0.1;k=0.001;

结果一：

1.误差最大值：

e =

7.9402e-004

2.误差图如下图所示：

B: a1=-1;b1=2;a=1;b=1;h=0.1;k=0.001; 结果二：

1.误差最大值：

e =

0.0682

2.误差图：

C: a1=-1;b1=-0.1;a=1;b=1;h=0.1;k=0.001; 结果三：

1.误差最大值：

e =

6.2221e+005

2.误差图：

五、实验结论：

通过上机实现，进一步直观了解流扩散方程的稳定具有很强的条件性，只要在a1,b1，h和 满足（*）式时才是稳定的，如结果一、二，否则会出现结果三的情形，误差相当大。本次实验，熟悉并掌握了差分格式的一般构造方法，理清了具体的步骤，提高了利用计算机解决问题的能力。

附：Matlab源代码：

1. function z=ft(x)%求下边界

z=exp(x/2);

2.function z=fx1(t)%求左边界

z=exp(-t);

3.function z=fx2(t)%求右边界

z=exp(1/2-t);

4.function z=f(x,t)%求右端函数

z=-exp(x/2-t);

5 .function z=fu(x,t)%求真解

z=exp(x/2-t);

6. function [X,T,z]=upwindL(a1,b1,a,b,h,k)%用迎风格式求解upwindL(1,2,1,1,0.1,0.1)

x=0:h:a;t=0:k:b;

[T,X]=meshgrid(t,x);

m=length(x);n=length(t);

r1=a1*k/h;r2=b1*k/h^2;

uu=zeros(m,n);%储存数值解

z=uu;%储存误差

for i=1:m%求下边界

uu(i,1)=ft(x(i));

end

for j=2:n%求左右边界

uu(1,j)=fx1(t(j));

uu(m,j)=fx2(t(j));

end

%迎风格式求内点，从下往上

if(a1>0)

for j=2:n

for i=2:m-1%从左往右

uu(i,j)=(1-r1-2*r2)*uu(i,j-1)+(r1+r2)*uu(i-1,j-1)+r2*uu(i+1,j-1)+k*f(x(i),t(j-1));%求数值解z(i,j)=abs(uu(i,j)-fu(x(i),t(j)));%求误差

end

end

else

for j=2:n

for i=2:m-1%从左往右

uu(i,j)=(1+r1-2*r2)*uu(i,j-1)+(r2-r1)*uu(i+1,j-1)+r2*uu(i-1,j-1)+k*f(x(i),t(j-1));%求数值解z(i,j)=abs(uu(i,j)-fu(x(i),t(j)));%求误差

end

end

end

%主函数，用于输出

7. [X,T,z]=upwindL(a1,b1,a,b,h,k);

mesh(T,X,z)

e=max(max(z))

title('误差图')

xlabel('x轴')

ylabel('t轴')

zlabel('z轴')