计算实验报告样板.docx
- 文档编号:30159692
- 上传时间:2023-08-05
- 格式:DOCX
- 页数:14
- 大小:312.80KB
计算实验报告样板.docx
《计算实验报告样板.docx》由会员分享,可在线阅读,更多相关《计算实验报告样板.docx(14页珍藏版)》请在冰豆网上搜索。
计算实验报告样板
班号:
学号:
07
姓名:
成绩:
实验报告
微分方程综合实验
一、实验目的
1、了解求解微分方程的计算步骤;
2、掌握使用Matlab计算微分方程精确解及数值解步骤;
3、掌握使用PDETool及编程的方法计算偏微分方程数值解;
二、实验内容与结果
1、假设湖的容积为V立方米,若污染物在湖中分布均匀,即每单位湖水在t时刻所含污染物为统一量度,以x(t)记湖中污染物的总量,以r记湖水流出湖的速率,设湖的容积为常量。
问题:
若停止向湖中的流入污染物,多长时间后湖水的污染程度可减低至初始状态的5%?
分析:
因为湖的容积为常量,则可知,湖水的净流出量=湖水的净流入量
则湖中污染物的改变量=流入的污染物-流出的污染物
可知,在时间间隔[t,t+Δt]内,污染物的浓度为:
流出的污染物为:
而流入的污染物为0,在该时间间隔内,污染物的改变量dx为:
改写上述方程为:
即:
令x的初始值为x0,可以利用下列命令计算式
(1)的精确解:
symsxt
x=dsolve('Dx-u*x=0','x(0)=x0','t');
(a)写出你的计算结果:
x=
x0*exp(u*t)
(b)假设V=4×1010m3,每天流出量r=4×108m3/d,则污染物降低到初始状态5%,即x(t)=0.05x0,根据你计算结果,计算其所需时间。
写出你的计算命令和结果。
计算命令:
symsxt
x=dsolve('Dx+0.01*x=0','x(0)=100','t');
结果:
x=
x0*exp(-1/100*t)
(c)利用下述命令作出该函数的图形
ezplot(y)
图形:
2、对于二阶常微分方程,常遇到不存在精确解的情况,但可以计算其数值解。
计算常微分方程的特解:
先将其化为等价的一阶方程组,令z1=y,z2=y’,则上式转化为:
以myode.m为文件名保存上述等价式:
myode.m文件
functionz=myode(x,y)
z
(1)=y
(2);
z
(2)=x*(1-y
(1)^2)*y
(2)-x*y
(1);
z=z';
注意,最后需以列向量的方式给出导数结果。
将求解区间和初始值带入,在主程序中调用myode.m,计算原方程的数值解:
[x,y]=ode23('myode',1,15,[1;0])
注意,数值解保存在y的第一列中。
执行上述过程,并用命令plot作出函数数值解的曲线。
(a)作出函数曲线图:
3、利用Matlab工具箱计算偏微分方程的数值解
泊松方程的是最简单的椭圆型PDE问题,该问题公式为:
求解区域为单位圆盘,边界调节为在边界圆盘上U=0,该问题的精确解为:
采用Matlab中的pdetool计算该偏微分方程。
在命令窗口中输入pdetool命令,在新窗口的下拉框中选择“GenericScalar”模式。
新窗口如下图所示:
(1)单击Option菜单,选择“Grid”和“Snap”选项,单击工具栏上
按钮,画一个园,若该园不是标准园,则双击该园,打开对话框,指定其圆心的位置(0,0)和半径大小(均为1)。
(2)单击
按钮设置边界模式。
分割的几何边界显示出来,并且外边界指定为默认设置,即Dirichlet边界条件:
。
本例中采用默认设置,若边界条件不同,可以通过双击边界打开对话框,在其中输入相应的边界条件。
(3)单击
按钮,在其中定义PDE系数c,a和f。
本例中,它们均为常数:
c=1,f=1,a=0
(4)单击
按钮或选择Mesh菜单中的“InitializeMesh”选项,初始化显示三角形网格。
请将划分网格后的图附于下面。
网格图:
(5)单击
按钮或Mesh菜单中的“RefineMesh”选项,对现有网格进行改进加密。
(6)单击
进行问题求解,Matlab中可以用图形表示问题的解。
将计算结果的图形附于下面。
计算结果平面图:
(7)单击
,打开PlotSelection对话框,利用该对话框,可以选择不同类型的解图。
将计算结果的三维图附于下面。
三维计算结果:
(8)比较数值解与精确解间的误差。
在PlotSelection对话框中的Color项选择Property下拉框中的Userentry,然后在Userentry编辑区输入Matlab表达式u-(1-x.^2-y.^2)/4,可以获得解的绝对误差的图形表示。
将误差图附于下面。
误差结果图:
4、采用数值解法计算上述泊松方程。
(1)首先对问题进行定义:
g='circleg';%单位圆
b='circleb1';%边界上为零条件
c=1;a=0;f=1;%标准椭圆方程系数
(2)产生初始的三角形网格
[p,e,t]=initmesh(g);
(3)迭代直至得到误差允许范围内的解
error=[];err=1;
whileerr>0.01,
[p,e,t]=refinemesh(g,p,e,t);%继续剖分
u=assempde(b,p,e,t,c,a,f);%计算数值解
exact=(1-p(1,:
).^2-p(2,:
).^2)/4;%计算精确解
err=norm(u-exact',inf);%计算误差
error=[errorerr];
end
(4)以图形方式显示计算结果
subplot(2,2,1),pdemesh(p,e,t);
subplot(2,2,2),pdesurf(p,t,u)
subplot(2,2,3),pdesurf(p,t,u-exact')
注意,计算结果保存在变量u中,p,e,t中为网格数据。
以图形化的方式显示计算结果。
结果图形:
5、求解正方形区域
上的热传导方程
初始条件为
边界条件为Dirichlet条件
。
这里是抛物型方程,其中
。
编写程序如下:
(1)问题定义
g='squareg';%定义正方形区域
b='squareb1';%边界上为零条件
(2)产生初始的三角形网格
[p,e,t]=initmesh(g);
(3)定义初始条件
u0=zeros(size(p,2),1);
ix=find(sqrt(p(1,:
).^2+p(2,:
).^2)<0.4);
u0(ix)=1
(4)在时间段0到0.1的20上求解
nframe=20;
tlist=linspace(0,0.1,nframe);
u1=parabolic(u0,tlist,b,p,e,t,c,a,f,d);
(5)动画显示结果
forj=1:
nframe
pdesurf(p,t,u1(:
j));
mv(j)=getframe;
end
movie(mv,10)
将温度场传递方向以动画形式显示。
结果动画:
6、求解正方形区域
上的波方程
初始条件为
,
,边界条件为在
上满足Dirichlet条件
,在
上满足Neumann条件
。
这是双曲型方程,其中
。
(1)问题定义
g='squareg';%定义正方形区域
b='squareb3';%定义边界
c=1;a=0;f=0;d=1;
(2)产生初始三角形网格
[p,e,t]=initmesh(g);
(3)定义初始条件
x=p(1,:
)';y=p(2,:
)';
u0=atan(cos(pi*x));
ut0=3*sin(pi*x).*exp(cos(pi*y));
(4)在时间段0到5的31个点上求解
n=31;
tlist=linspace(0,5,n);
uu=hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d);
(5)动画显示结果
forj=1:
n
pdesurf(p,t,uu(:
j));
mv(j)=getframe;
end
movie(mv,10)
将计算结果以动画的形式表现出来。
计算结果:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算 实验 报告 样板