偏微分方程数值实验报告八

实验题目：利用有限差分法求解

u

( x)

u(x)

f (x),

u(

1)

0, u(1)

0.

真解为

u( x)

e x2

(1

x2 )

实现算法：对于两点边值问题

d2u

f , x

l ,

dx 2

(1)

u(a),u(b)

,

其中 l ( a, b) (a

b), f 为 l

[ a,b] 上的连续函数，

，

为给定常数 .

其相应的有限差分法的算法如下：

1.对求解区域做网格剖分，得到计算网格.在这里我们对区间l 均匀剖分 n 段，每个剖分单元

b

a

的剖分步长记为 h

.

n

2.对微分方程中的各阶导数进行差分离散，得到差分方程

.运用的离散方法有：

方法一 :用待定系数和泰勒展开进行离散

d 2u( xi )

i 1u( xi 1)

i u( xi )

i 1u( xi 1)

d( xi )2

方法二：利用差商逼近导数

d2u( xi )u( xi 1 )2u( xi )u( xi 1 )

d( xi )2h2

将(2) 带入 (1)可以得到

u(xi 1) 2u(xi ) u(xi 1 )

)

Ri (u) ,

h2

f ( xi

其中 Ri (u) 为无穷小量，这时我们丢弃 Ri

(u) ，则有在 xi 处满足的计算公式：

u(xi 1) 2u( xi ) u( xi 1 )

1,..., n 1

h2

f ( xi )， i

3.根据边界条件，进行边界处理

.由 (1)可得

u0

, un

(2)

(3)

(4)

称(3)(4)为逼近 (1) 的差分方程，并称相应的数值解向量

U n 1 为差分解， ui 为 u( xi ) 的近似值 .

4.最后求解线性代数方程组，得到数值解向量U n 1 .

实验题目： 用 Lax-Wendroff格式求解方程：

ut

2ux

0, x (0,1), t

0,

u( x,0)

1

sin 2

x, x

[ 0,1],

u(1,t )

1

sin 4

t.

(精确解 u

1 sin 2 ( x

2t). )

数值边值条件分别为：

(a) u0n 1

u0n 2 (u1n

u0n ).

h

(b ) u0n

u1n.

(c) u0n 1 - 2u1n 1

u2n 1

0.

请将计算结果与精确解进行比较。

实现算法：

1. 网格剖分：

对求解区域 G

(0,1)

(0,T ] 作均匀网格剖分 .

节点：

xj

jh ,

j

0,1,..., N

t k

jh ,

k

0,1,..., M

其中空间和时间步长：

h

1 ,

T .

N

M

算法实现

u(xi , tk 1 ) 在节点 ( xi ,t k ) 处作泰勒级数展开

u(xi , tk 1)

u(xi , tk )

[

u

k

2

2u k

O (

3

)

]i

[

t

2 ]i

t

2!

考虑在节点 ( x j , tk ) 处( 1)的微分方程，有：

u

a

u

0.

t

x

2u

(

a

u

a

(

u

2

2u

t

2

)

) a

x

2 .

t

x

x

t

将上述两式代入( 2)式，得

u( xi ,t k 1 ) u( xi ,t k ) a [ u ] ik

a2 2 [

2 u2 ]ik

O ( 3 )

x

2

x

对 x 的一阶、二阶导数用中心差商代替

1)

2)

u k

1

2

[

x ]i

2h

[( u( xi

1, tk )

u( xi 1, tk ))]

O(h

)

[

2u k

1

2u( xi ,tk )

u( xi 1, tk ))]

2

)

x

2 ]i

h

2 [( u(xi 1 ,t k )

O (h

代入整理后得到

u( xi ,tk 1 ) u( xi ,t k )

a [ u( xi 1 ,t k ) u( xi 1 ,t k )]

2h

a

2

2

O( h2 )

2 h2 )

O( 3 )

2

[ u( xi 1 ,tk )

2u( xi , tk )

u( xi

1 , tk )]

O (

2h

略去误差项，以

uik 代替 u( xi ,t k ) ，得到如下差分格式

a (uik 1

2

2

uik 1

uik

uik 1 )

a

2 (uik 1

2uik

uik 1)

( 3)

2h

2h

( 3)式就是 Lax-Wendroff 格式，其截断误差为

O(

2

h2 ) ，节点如图

| a |

，就得到(

1)式的 Lax-Wendroff

格式的公式

令 r

h

uik 1

uik

r (uik 1

uik 1 )

r 2

(uik 1

2uik

uik

1)

( 4)

2

2

( 4)式是二阶精度的差分格式.

程序代码：

function[X,T,U] = advection_fd1d (NS ,NT ,pde,bd)

WAVE_EQUATION_FD1D利用有限差分方法计算一维双曲线方程

输入参数：

NS整型，空间剖分段数

NT整型，时间剖分段数

pde结构体，带求解的微分方程模型的已知数据，

%如边界、初始、系数和右端项等条件.

%bd数值边值条件

输出参数 :

X长度 NS+1 的列向量，空间网格剖分

T长度 NT+1 的行向量，时间网格剖分

%U (NS+1)*(NT+1)矩阵， U(:,i)表示第 i个时间层网格剖分上的数值解

[X,h] = (NS);

[T,tau] = (NT);

N = length(X);M = length(T);

U = zeros(N,M);

初值条件

U(:,1) = (X);

a = ;

r = a*tau/h;

边值条件

ifa >= 0% 左边值条件

U(1,:) = (T)

else

U(end,:) = (T)%右边值条件

end

for i = 2:M

U(2:end -1,i) =U(2:end-1,i-1)-r*(U(3:end,i-1)-U(1:end-2,i-1))/2+...

r^2*(U(3:end,i-1)-2*U(2:end-1,i-1)+U(1:end-2,i-1))/2;

switch(bd)

case { 'a0'}

a0();

case { 'b'}

b();

case { 'c'}

c();

otherwise

disp(['Sorry, I do not know your ', bd]);

end

end

functiona0()

U(1,i)=U(1,i-1)-r*(U(2,i-1)-U(1,i-1));

end

functionb()

U(1,i)=U(2,i-1);

end

functionc()

U(1,i)=2*U(2,i)-U(3,i);

end

end

functionpde = model_data ()

%MODEL_DATA数据模型

TI = 0;

TF = 1;

SI = 0;

SF = 1;

pde = struct('u_exact',@u_exact, 'u_initial'

'u_left',@u_left,'u_right',@u_right,

,@u_initial,

'time_grid', ...

...

@time_grid,

'space_grid'

,@space_grid,'advection_fd1d_error'

,@advection_fd1d_error,

'a

' ,-2);

function[T,tau] = time_grid(NT)

T = linspace(TI,TF,NT+1);

tau = (TF-TI)/NT;

end

function[X,h] = space_grid(NS)

X = linspace(SI,SF,NS+1)'

h = (SF-SI)/NS;

end

functionU = u_exact(X,T)

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

U = 1+sin(2*pi*(x+2*t));

end

functionu = u_initial (x)

u = 1+sin(2*pi*x);

end

functionu = u_right(t)

=1+sin(4*pi*t);

end

end

functionshowsolution (X,T,U)

%% SHOWSOLUTION以二元函数方式显示数值解

% 输入参数

% X 长度为 NS +1的列向量，空间网格剖分N

% T 长度为 NT +1的行向量，时间网格剖分M

U N*M 矩阵， U(:,i) 表示第 i 个时间层网格部分上的数值解

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

mesh (x,t,U');

xlabel ('X' );

ylabel ('T' );

zlabel ('U(X,T)');

end

functionshowvarysolution (X,T,U,UE)

%% SHOWVARYSOLUTION显示数值解随着时间的变化

输入参数

%X长度为 NS +1的列向量，空间网格剖分N

%T长度为 NT +1的行向量，时间网格剖分M

U N*M矩阵， U(:,i) 表示第 i 个时间层网格部分上的数值解

M = size (U ,2) ;

figure

xlabel ('X' );

ylabel ('U' );

s = [X(1) ,X(end),min(min(U)),max(max(U))];

axis (s);

fori = 1:M

plot (X,U(:,i));

axis (s);

pause ;

title (['T=',num2str(T(i)),' 时刻的温度分布' ])

End

一维双曲线有限差分方法主测试脚本

pde=model_data()

[X,T,U]=advection_fd1d(100,200,pde,'a');

UE=(X,T);

showvarysolution(X,T,U,UE);%以随时间变化方式显示数值解

showsolution(X,T,U);

%以二元函数方式显示数值解

[X,T,U]=advection_fd1d(100,200,pde,

'b');

UE=(X,T);

showvarysolution(X,T,U,UE);

%以随时间变化方式显示数值解

showsolution(X,T,U);%以二元函数方式显示数值解

[X,T,U]=advection_fd1d(100,200,pde, UE=(X,T);

'c');

showvarysolution(X,T,U,UE);%以随时间变化方式显示数值解

showsolution(X,T,U);%以二元函数方式显示数值解