差分格式稳定性及数值效应比较实验.doc

差分格式稳定性及数值效应比较实验

5090719044 张 赟 F0907102

一 实验目的：
1.以一阶线性双曲线方程为例，使用Matlab工具分析4种差分格式的误差。
2.了解4种差分格式的稳定性

二 实验问题：
对于一阶线性双曲型方程：

取a=1,2,4, h=0.1, τ=0.08, 对不同的差分格式（迎风格式，Lax-Friedrichs格式，Lax-Wendroff格式，修正迎风格式）及不同的a值进行迭代计算。通过将计算结果与精确解来进行比较，来讨论分析差分格式的稳定性。

三 实验原理：
1. 迎风格式：
这种格式的基本思想是简单的，就是在双曲型方程中关于空间偏导数用在特征线方向一侧的单边差商来代替，格式如下：

运算格式：

2. Lax-Friedrichs格式：

运算格式：

3. Lax-Wendroff格式：
这种格式构造是采用 Taylor 级数展开和微分方程本身得到，运算格式：

4. 修正迎风格式（目标点范围跟踪格式）：

其中是取整数部分，=。根据之后的理论分析可以得到这是一个无条件稳定结构。

四 四种格式理论分析：
通过求差分格式的增长因子G(τ , k)，来判定差分格式是否稳定。

1. 迎风格式：
记，则，
得,
即。
所以。
则在，满足von Neumann条件，格式稳定。
以下格式用相同方法求解稳定性条件。

2. Lax-Friedrichs格式：
，在时稳定。

3. Lax-Wendroff格式：
，在时稳定。

4. 修正迎风格式（目标点范围跟踪格式）：
，
其中，的成立条件为。
而恒成立，故格式无条件稳定。

五 实验结果：
a=1（）

迎风格式 Lax-Friedrichs格式

Lax-Wendroff格式 修正迎风格式

a=2（）

迎风格式 Lax-Friedrichs格式

Lax-Wendroff格式 修正迎风格式
a=4（）

迎风格式 Lax-Friedrichs格式

Lax-Wendroff格式 修正迎风格式

六 总结：
本次实验，通过4种差分格式求解T=4时的解并与解析解画图比较，可以看出：
(1) a=1（aλ=0.8＜1）时，迎风格式，Lax-Friedrichs格式，修正迎风格式的计算结果与解析解近似情况较好，而Lax-Wendroff格式则在间断点处出现了波前波，形成双波现象，这符合Lax-Wendroff格式为二阶迭代格式的性质。
(2) a=2（aλ=1.6＞1）时，迎风格式，Lax-Friedrichs格式，Lax-Wendroff格式都出现了比较强烈的震荡，震荡区域主要在(2,3)的区间内。这三种震荡中，Lax-Friedrichs格式震荡较小，为级别，迎风格式与Lax-Wendroff格式的震荡则较大，为级别。与之相对应的是修正迎风格式，保持着稳定的性质，不过间断点从4到了13。
(3) a=4（aλ=3.2＞1）时，迎风格式，Lax-Friedrichs格式，Lax-Wendroff格式的震荡更加强烈，分别到达，，的级别。震荡区间也变为(0,3)。修正迎风格式则仍然保持着原有的稳定性不变，只是间断区间变为(60,90)。
由上得出，稳定性对差分格式求解偏微分方程有重大意义。一个差分格式是否好，是否可用，首先要判定它是否稳定并找到稳定性条件。修正迎风格式强大的稳定性在解决一阶线性双曲线方程中有着很强的实用价值。

七 程序：
迎风格式：
function yingfeng(a,h,t,minx,maxx)
m=(maxx-minx)/h;
T=4;
p=t/h;
n=T/t;
u1=ones(m+n+1,1);
u1(n+1:m+n+1)=0;
u2=u1;
for i=1:1:n
for j=i+1:1:m+n+1
u2(j)=a*p*u1(j-1)+(1-a*p)*u2(j);
end
u1=u2;
end
y1=u2(n+1:m+n+1);
x1=minx:h:maxx;
plot(x1,y1,'--')

Lax-Friedrichs格式：
function Friedrichs(a,h,t,minx,maxx)
m=(maxx-minx)/h;
T=4;
p=t/h;
n=T/t;
u1=ones(m+2*n+1,1);
u1(n+1:m+2*n+1)=0;
u2=u1;
for i=1:1:n
for j=i+1:1:m+2*n+1-i
u2(j)=0.5*(1+a*p)*u1(j-1)+0.5*(1-a*p)*u2(j+1);
end
u1=u2;
end
y1=u2(n+1:m+n+1);
x1=minx:h:maxx;
plot(x1,y1,'--')

Lax-Wendroff格式：
function Wendroff(a,h,t,minx,maxx)
m=(maxx-minx)/h;
T=4;
p=t/h;
n=T/t;
u1=ones(m+2*n+1,1);
u1(n+1:m+2*n+1)=0;
u2=u1;
for i=1:1:n
for j=i+1:1:m+2*n+1-i
u2(j)=0.5*a*p*(1+a*p)*u1(j-1)+(1-a*p)*(1+a*p)*u2(j)+0.5*a*p*(a*p-1)*u1(j+1);
end
u1=u2;
end
y1=u2(n+1:m+n+1);
x1=minx:h:maxx;
plot(x1,y1,'--')

修正迎风格式：
function gaijinyingfeng(a,h,t,minx,maxx)
m=(maxx-minx)/h;
T=4;
p=t/h;
n=T/t;
u1=ones(m+n+1,1);
u1(n+1:m+n+1)=0;
u2=u1;
for i=1:1:n
for j=i+floor(a*p)+1:1:m+n+1
u2(j)=(a*p-floor(a*p))*u1(j-floor(a*p)-1)+(1-a*p+floor(a*p))*u2(j-floor(a*p));
end
u1=u2;
end
y1=u2(n+1:m+n+1);
x1=minx:h:maxx;
plot(x1,y1,'--')