偏微分方程数值及matlab实验报告8.docx
- 文档编号:27514686
- 上传时间:2023-07-02
- 格式:DOCX
- 页数:10
- 大小:790.93KB
偏微分方程数值及matlab实验报告8.docx
《偏微分方程数值及matlab实验报告8.docx》由会员分享,可在线阅读,更多相关《偏微分方程数值及matlab实验报告8.docx(10页珍藏版)》请在冰豆网上搜索。
偏微分方程数值及matlab实验报告8
实验题目:
用Lax-Wendroff格式求解方程:
(1)
(精确解
)
数值边值条件分别为:
(a)
(b)
(c)
请将计算结果与精确解进行比较。
实现算法:
1.网格剖分:
对求解区域
作均匀网格剖分.
节点:
其中空间和时间步长:
2.算法实现
将
在节点
处作泰勒级数展开
(2)
考虑在节点
处
(1)的微分方程,有:
将上述两式代入
(2)式,得
对
的一阶、二阶导数用中心差商代替
代入整理后得到
略去误差项,以
代替
,得到如下差分格式
(3)
(3)式就是Lax-Wendroff格式,其截断误差为
,节点如图
令
,就得到
(1)式的Lax-Wendroff格式的公式
(4)
(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]=pde.space_grid(NS);
[T,tau]=pde.time_grid(NT);
N=length(X);M=length(T);
U=zeros(N,M);
%初值条件
U(:
1)=pde.u_initial(X);
a=pde.a;
r=a*tau/h;
%边值条件
ifa>=0%左边值条件
U(1,:
)=pde.u_left(T)
else
U(end,:
)=pde.u_right(T)%右边值条件
end
fori=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,Idonotknowyour',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_initial,...
'u_left',@u_left,'u_right',@u_right,'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)
u=1+sin(4*pi*t);
end
end
functionshowsolution(X,T,U)
%%SHOWSOLUTION以二元函数方式显示数值解
%输入参数
%X长度为NS+1的列向量,空间网格剖分N
%T长度为NT+1的行向量,时间网格剖分M
%UN*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
%UN*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(0.01);
title(['T=',num2str(T(i)),'时刻的温度分布'])
End
%一维双曲线有限差分方法主测试脚本
pde=model_data()
[X,T,U]=advection_fd1d(100,200,pde,'a');
UE=pde.u_exact(X,T);
showvarysolution(X,T,U,UE);%以随时间变化方式显示数值解
showsolution(X,T,U);%以二元函数方式显示数值解
[X,T,U]=advection_fd1d(100,200,pde,'b');
UE=pde.u_exact(X,T);
showvarysolution(X,T,U,UE);%以随时间变化方式显示数值解
showsolution(X,T,U);%以二元函数方式显示数值解
[X,T,U]=advection_fd1d(100,200,pde,'c');
UE=pde.u_exact(X,T);
showvarysolution(X,T,U,UE);%以随时间变化方式显示数值解
showsolution(X,T,U);%以二元函数方式显示数值解
THANKS!
!
!
致力为企业和个人提供合同协议,策划案计划书,学习课件等等
打造全网一站式需求
欢迎您的下载,资料仅供参考
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分方程 数值 matlab 实验 报告
![提示](https://static.bdocx.com/images/bang_tan.gif)