偏微分方程数值实验报告八
实验题目：利用有限差分法求解
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);%以二元函数方式显示数值解
如果觉得《matlab的数值求解实验报告 偏微分方程数值及matlab实验报告》对你有帮助，请点赞、收藏，并留下你的观点哦！
