数学建模第4章.docx
- 文档编号:2085696
- 上传时间:2022-10-26
- 格式:DOCX
- 页数:31
- 大小:339.80KB
数学建模第4章.docx
《数学建模第4章.docx》由会员分享,可在线阅读,更多相关《数学建模第4章.docx(31页珍藏版)》请在冰豆网上搜索。
数学建模第4章
欧拉方法和龙格-库塔方法求微分方程数值解,画出解的图形,对结果进行比较分析。
以下方程供选择:
d)幂级数解
解:
原方程化为
记,并用t代替x,则原方程化为:
;且;
于是可以用龙格-库塔法求解。
【模型求解】:
用Matlab作龙格-库塔法求解:
%chapter4—
%此函数是微分方程组
functionXdot=ch41dfun(t,x)
Xdot=[x
(2),-cos(t)*x
(1)]';
%
functionI=ch41d(a)
x0=[1,0]';
[t,x]=ode45('ch41dfun',[0,a],x0);
y=x(:
1);
plot(t,y,'r');holdon;gtext('y');
y1=1-1/jiecheng
(2)*t.^2+2/jiecheng(4)*t.^4-9/jiecheng(6)*t.^6+55/jiecheng(8)*t.^8;
plot(t,y1,'b');gtext('y1');
[t,y,y1]
holdoff;
运行程序可以得到
图1
>>ch41d(10)
图2
>>ch41d(15)
图3
>>ch41d(20)
图4
【结果分析】:
由图1得:
y(龙格-库塔方法)和y1(级数近似解)在0到大约的区间内是完全吻合的,从x=之后,两条曲线开始分离。
之后y的变化趋势可见图2-图4,呈幅度越来越大的上下震荡。
小型火箭初始重量为1400kg,其中包括1080kg燃料。
火箭竖直向上发射时燃料的燃烧率为18kg/s,由此产生32000N的推力,火箭引擎在燃料用尽时关闭。
设火箭上升时的空气阻力正比于速度的平方,比例系数为0.4kg/m。
求引擎关闭瞬间火箭的高度、速度、加速度,及火箭到达最高点时的高度和加速度,并画出高度、速度、加速度随时间变化的图形。
解:
解:
先解出速度的微分方程如下:
clear
symsamtv
m=1400-18*t;
a=**v^2)/m
a=
(18280+882/5*t-2/5*v^2)/(1400-18*t)
利用龙哥库塔算法画出从0时刻到60秒(燃料用尽)的图像
odefun=@(t,v)(-((882*t)/5-(2*v^2)/5+18280)/(18*t-1400));
[t,v]=ode45(odefun,[0:
:
60],[0]);
subplot(2,2,1)
plot(t,v);
gridon
a=diff(v)/;
t2=[0:
:
];
subplot(2,2,2)
plot(t2,a);
gridon
s=cumsum(v).*;
subplot(2,2,3)
plot(t,s);
gridon
所以v(60)=s
a(60)=0.9127m/s2
s(60)=12200m
某容器盛满水后,底端直径为d0的小孔开启(如图1),根据水力学知识,当水面高度为h时,谁从小孔中流出的速度为v=*(g*h)^(其中g为重力加速度,问哦小孔收缩系数)
1)若容器为倒圆锥形(如图1),现测得容器高和上底直径都为,小孔直径d为3cm,为水从小孔中流完需要多少时间;2min时水面高度是多少。
2)若容器为倒葫芦形(如图2),现测得容器高,小孔直径d为3cm,由底端(记x=0)向上每隔测出容器的直径D(m)如表1所示,问水从小孔中流完需要多少时间;2min时水面高度是多少。
图1
X/m
0
D/m
表1
【模型分析】
由题知,水从小孔中流出,不仅与容器有关,还与水流速度v=*(2*g*h)^有关。
第一小题容器是圆锥形,比较规则,但是由于水不断从小孔流出,容器中水的高度是不断变化的,水流速度没有一定的公式,所以要用到微积分解决。
由
(1)知,水面直径等于水深。
水深为h时,流量为(π/4)d^2*(gh)^,
*(g*h)^*π*(d0/2)^2*dt=π/4*h^2*dh
则水深下降dh所需时间:
dt=-[(π/4)h^2*dh]/[(π/4)d^2*(gh)^]=-[h^*dh]/[^2*(g)^]
水深由至0定积分得水从小孔流完的时间:
T(其中已知d=,g=*s(-2)
对于第二问:
设两分钟(120S)后水深为Xm,由
dt=-[(π/4)h^2*dh]/[*(π/4)*d^2*(gh)^]=-[h^*dh]/[^2*(g)^]
则=X^[*d^2*(g)^]
以d=,g=*s(-2代入上式得水深:
X
第二小题容器为倒葫芦形,比较不规则,比较复杂,不仅要考虑水不断从小孔流出,容器中水的高度是不断变化的,水流速度没有一定的公式,所以要用到微积分解决,还要注意表1的倒葫芦形的不断变化,水深的高度变化是不规则的但仍可以用微积分。
由
(2)知容器高,水深为h时,流量为(π/4)d^2*(gh)^,由于不同高度,倒葫芦形半径不同,用欧拉方程和龙格—库塔方法则水深下降dh所需时间:
dt=t(k+1)-t(k)=-[(π/4)h^2*dh]/[(π/4)d^2*(gh)^]=-[h^*dh]/[^2*(g)^]
然后利用循环fork=1:
length(L),t(k)=((h(k+1)-h(k))*(π/4)*d(k)^2)/*(π/4)*d^2*(g(k)))^,T=sum(t).可以求得水从小孔流完的总时间。
对于第二问:
设两分钟(120S)后水深为Xm,由
S=0,利用条件,当120-s<时s=s+t(k),x(k)
把d=,g=*s(-2)代入上式得水深:
X
【模型求解】
开始条件h=,d0=,g=*s(-2),d=
水深为h时,流量Q为(π/4)d^2*(gh)^,则水深下降dh所需
时间:
dt=-[(π/4)h^2*dh]/[(π/4)d^2*(gh)^]=-[h^*dh]/[^2*(g)^]
水深由至0定积分得水从小孔流完的时间:
T(其中已知d=,s)
设两分钟(120S)后水深为Xm,由
dt=-[(π/4)h^2*dh]/[*(π/4)*d^2*(gh)^]=-[h^*dh]/[^2*(g)^]
则=X^[*d^2*(g)^]
以d=,g=代入上式得水深:
X
用欧拉方程和龙格—库塔方法则水深下降dh所需时间:
dt=t(k+1)-t(k)=-[(π/4)h^2*dh]/[(π/4)d^2*(gh)^]=-[h^*dh]/[^2*(g)^]
k1=*sqrt(g*(x(n)))*d^2/*x(n)^8+*x(n)^*x(n)^6+*x(n)^*x(n)^4+*x(n)^*x(n)^2+*x(n)+.014892)^2;
k2=*sqrt(g*(x(n))-h*k1)*d^2/*(x(n)-h*k1)^8+*(x(n)-h*k1)^*(x(n)-h*k1)^6+*(x(n)-h*k1)^*(x(n)-h*k1)^4+*(x(n)-h*k1)^*(x(n)-h*k1)^2+*(x(n)-h*k1)+.014892)^2;
x(n+1)=x(n)-h*(k1+k2)/2;
【附录】
(1)
g=;
d=;
symsh
y=-h^/*d^2*sqrt(2*g));
T=int(y,h,,0);
eval(T)
x=((T-120)**d^2*sqrt(2*g)))^;
eval(x)
运行结果如下
>>sy1
ans=
ans=
(2)
clear;
g=;
d=;
k1=0;
k2=0;
h=;
x
(1)=;
forn=1:
1000
k1=*sqrt(g*(x(n)))*d^2/*x(n)^8+*x(n)^*x(n)^6+*x(n)^*x(n)^4+*x(n)^*x(n)^2+*x(n)+.014892)^2;
k2=*sqrt(g*(x(n))-h*k1)*d^2/*(x(n)-h*k1)^8+*(x(n)-h*k1)^*(x(n)-h*k1)^6+*(x(n)-h*k1)^*(x(n)-h*k1)^4+*(x(n)-h*k1)^*(x(n)-h*k1)^2+*(x(n)-h*k1)+.014892)^2;
x(n+1)=x(n)-h*(k1+k2)/2;
end
x(300)
t=0:
h:
1000*h;
plot(t,x);
axis([0,400,0,]);
运行结果如下
>>sy2
ans=
水从小孔流完需要263s,2min时水面高度是点头绪
放射性废物的处理:
将装满放射性废物的圆桶扔到水深300ft的海底,圆桶体积55gal,装满放射性废物时的圆桶重量为,在海中浮力为。
此外,下沉时圆桶还要受到海水的阻力,阻力与下沉速度成正比,比例系数为s/ft。
大量破坏性实验发现当圆桶速度超过40ft/s时,就会因与海底冲撞而发生破裂。
(1)建立解决上述问题的微分方程模型。
(2)用数值和解析两种方法求解微分方程,并回答谁赢得了官司
解:
【模型假设】
假设装满放射性废物的圆桶在水中仅受重力、浮力和阻力的作用。
设圆桶下沉距离为S,速度为v,由于阻力和下沉速度成正比,比例系数为k=**。
设圆桶质量为m,重力为G,浮力为F,阻力为R,则有R=kv=***v/。
根据牛顿第二定律可得微分方程组
对此微分方程组,利用龙格-库塔方法能给出数值解。
(一)用数值法求解微分方程
【算法设计】
有两个思路:
思路一是利用龙格-库塔方法直接解上面关于下沉距离S和下沉速度v的微分方程组。
思路二是利用龙格-库塔方法解关于速度v和时间t的一阶微分方程,即,再利用辛普森公式对速度求积分已得到下沉的距离S;分别求出S和v与时间的关系后,对比同一时间下的S和v,通过判断下沉距离S=300ft时速度是否大于40ft/s或判断速度达到40ft/s时下沉的距离是否大于300ft(实际上下沉距离不可能大于海深,这里是只是一种数学上的处理),则可以判断圆桶与海底冲撞时是否会破裂,从而判断官司的胜负。
思路一:
【程序】
functiondx=waste(t,x)%建立名为waste的函数M文件
m=*;
G=**;%圆桶所受重力
F=**;%圆桶所受浮力
R=***x
(2)/;%圆桶所受阻力
dx=[x
(2);(G-F-R)/m];%以向量形式表示方程
clearall
ts=0:
:
25;%终点试探
x0=[0,0];%初始值
opt=odeset;%默认误差限
[t,x]=ode45(@waste,ts,x0);%使用5级4阶龙格—库塔公式计算
[t,x]%输出t,S(t),v(t)
plot(t,x(:
2)),grid,title('v(t)-t曲线');
%作出v(t)图形,以便考察速度变化趋势
xlabel('t');
ylabel('v(t)');
figure;%在新窗口中作图
plot(t,x),axis([0250350]),grid,%按照数值输出作S(t),v(t)的图形
gtext('S(t)'),gtext('v(t)')
holdon,Smax=300*;plot([0,25],[Smax,Smax],'r-'),holdoff
%作出S=Smax的图形,以便比对
holdon,Vmax=40*;plot([0,
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数学 建模
![提示](https://static.bdocx.com/images/bang_tan.gif)