实验4常微分方程数值解.docx
《实验4常微分方程数值解.docx》由会员分享，可在线阅读，更多相关《实验4常微分方程数值解.docx（20页珍藏版）》请在冰点文库上搜索。
实验4常微分方程数值解
实验4常微分方程数值解
化学工程系化32
坂井优（日本留学生）2013080091
【实验目的】
1.掌握用MATLAB软件求微分方程初值问题数值解的方法；
2.通过实例学习用微分方程模型解决简化的实际问题；
3.了解欧拉方法和龙格-库塔方法的基本思想和计算公式，稳定性等概念。
【实验内容】
1.题目5：
（题目略）
模型及其求解：
设圆桶的质量为m，圆桶受到的重力为G，在海水中受到的浮力为F1，受到的阻力为F2，圆桶下沉的速度为v，加速度为a，距海面的距离为s，设重力加速度为g。
因为圆桶所受海水阻力与下沉速度成正比，设比例系数为k.
dv/dt=（1-m1/m-kv/m）g
ds/dt=v
将圆桶与海底接触时的速度求解出来，与极限速度进行比较，若大于极限速度，则圆桶与海底碰撞后会破裂；若小于极限速度，则圆桶与海底碰撞后不会破裂。
源程序如下：
functiondv=Cylinder（t,v）
m1=213.3403;m=239.245;k=0.1191;g=9.80;
p=1-m1/m-k*v/m;
dv=p*g;
end
clc;clearall;
ts=0:
0.1:
1200;
n=length（ts）;
v0=0;
opt=odeset（'RelTol',1e-6,'AbsTol',1e-9）;
[t,v]=ode45（@Cylinder,ts,v0,opt）;
[t,v];
plot（t,v,'r-','LineWidth',2）,grid;
xlabel（'t','FontSize',16）;ylabel（'v（t）','FontSize',16）;
title（'速度随时间变化的关系','FontSize',16）;
gtext（'v（t）','FontSize',16）;pause;
s=cumtrapz（t,v）;
plot（t,s,'LineWidth',2）,grid;
xlabel（'t','FontSize',16）;ylabel（'s（t）','FontSize',16）;
title（'距离随时间变化的关系','FontSize',16）
gtext（'s（t）','FontSize',16）;
pause
H=91.44;v0=12.192;
plot（v,s,'LineWidth',2）;axis（[020,0150]）;
line（[0,20],[H,H],'Color','y','LineStyle','--'）;
line（[v0,v0],[0,200],'Color','g','LineStyle','--'）;
xlabel（'v','FontSize',16）;ylabel（'s','FontSize',16）;
title（'距离与速度的关系','FontSize',16）
gtext（'s~v','FontSize',16）;
图：
计算结果：
由图像可知，圆桶在下沉到海底时速度已经超过圆桶破裂时的极限速度，故其与海底碰撞后会发生破裂。
即这场官司的赢家是工程师。
2.题目6：
一只小船度过宽为d的河流，目标是起点A正对着的B点，已知河水流速v1与船在静
止的水中的速度v2之比为k。
（1）建立描述小船航线的微分方程模型，求其解析解；
（2）设d=100m,v1=1m/s,v2=2m/s，用数值方法求渡河所需时间、任意时刻小船位
置及航行曲线，作图，并与解析解比较；
（3）若流速v1=0,0.5,1.5,2（m/s），结果将如何。
模型及其求解：
假设驾驶船者并不知道水流速度，如果他知道，则只要有v1故船头的方向始终对着目标B点，如书上图。
以B点为原点可得速度v的x，y方向的分量的如下关系式：
由此建立一个常微分方程组，初始条件为。
（x,y）=（d,0）。
（1）-求其解析解：
由本题的微分方程表达式，1、2两式相除得：
令
=p,从而dx=pdy+ydp,于是上式变为
=
于是ln（
）=kln（cy）
解得：
故由
=p,x=
又初始值为A点的坐标：
（d,0）；
即：
x=
这就是小船航线的解析解表达式。
（2）代入d=100m,v1=1m/s,v2=2m/s，
用matlab进行求解，源程序如下：
建立boat.m
functiondx=boat（t,x）
m=sqrt（x
（1）^2+x
（2）^2）
a=1;b=2
dx=[a-b*x
（1）./m;-b*x
（2）./m];
end
ts=0:
80;
x0=[0;-100];
[t,x]=ode23s（@boat,ts,x0）;
[t,x]
plot（x（:
1）,x（:
2））;
grid;
title（'y-x曲线'）;
xlabel（'x/m'）;ylabel（'y/m'）;
plot（t,x（:
1））;
title（'x-t曲线'）;xlabel（'t/s'）;ylabel（'x/m'）;
grid;
plot（t,x（:
2））;
title（'y-t曲线'）;xlabel（'t/s'）;ylabel（'y/m'）;
grid,;
t
x
y
0
0
-100
1
0.990036802
-98.00000216
2
1.959741409
-96.00021829
3
2.908828047
-94.00077163
4
3.836821656
-92.00197073
5
4.742902446
-90.00418957
…
…
…
64
2.680322705
-0.283146382
65
1.687817487
-0.111541602
66
0.690109125
-0.018499931
67
-2.28E-08
0
68
2.15E-08
0
69
-6.19E-10
0
…
…
…
计算结果
部分表格如下：
图：
作其解析解：
x=
图像k=0.5,d=100；
源程序如下：
在MATLAB中编写如下M文件，命名为guohe.m：
fuctionx=guohe（y）
x=-0.5.*（-100）.^0.5*y.^0.5+0.5.*（-100）.^（-0.5）.*y.^1.5
编写如下程序作图：
y=[0:
-0.1:
-100];%以-0.01为步长
x=guohe（y）;%算出每个节点上相应的x的值
plot（x,y）,grid,%作y-x的图形
gtext（'x'）;
gtext（'y'）;
计算结果：
y-x的解析解曲线
（3）v1=0时
将原程序中v1=1改为v1=0即可
源程序如下：
建立boat.m
functiondx=boat（t,x）
m=sqrt（x
（1）^2+x
（2）^2）
a=0;b=2
dx=[a-b*x
（1）./m;-b*x
（2）./m];
end
ts=0:
50;
x0=[0;-100];
[t,x]=ode23s（@boat,ts,x0）;
[t,x]
plot（x（:
1）,x（:
2））;
grid;
title（'y-x曲线'）;
xlabel（'x/m'）;ylabel（'y/m'）;
plot（t,x（:
1））;
title（'x-t曲线'）;xlabel（'t/s'）;ylabel（'x/m'）;
grid;
plot（t,x（:
2））;
title（'y-t曲线'）;xlabel（'t/s'）;ylabel（'y/m'）;
grid;
计算结果
从图中可以看出，船是沿直线匀速从A驶向B的，这是水速为0的必然结果，从而也证明了模型本身的正确性。
v1=0.5时
将原程序中v1=1改为v1=0.5即可,由于ode23计算相对太慢，因此采用ode15计算。
源程序如下：
建立boat.m
functiondx=boat（t,x）
m=sqrt（x
（1）^2+x
（2）^2）
a=0.5;b=2
dx=[a-b*x
（1）./m;-b*x
（2）./m];
end
ts=0:
70;
x0=[0;-100];
[t,x]=ode15s（@boat,ts,x0）;
[t,x]
plot（x（:
1）,x（:
2））;
grid;
title（'y-x曲线'）;
xlabel（'x/m'）;ylabel（'y/m'）;
plot（t,x（:
1））;
title（'x-t曲线'）;xlabel（'t/s'）;ylabel（'x/m'）;
grid;
plot（t,x（:
2））;
title（'y-t曲线'）;xlabel（'t/s'）;ylabel（'y/m'）;
grid;
计算结果：
与v1=1m/s时相比，小船在x方向上行使的最大值减小了近一半。
这是因为河水流速减小后，船在x方向上的分速度减小的缘故。
与v1=0时比较，小船过河花的时间并没有增加很多，但路程长了很多，这说明小船的平均实际速度大于静水船速v2。
这是因为实际速度为河水与船静水速度的矢量叠加造成的。
v1=1.5时
将原程序中v1=1改为v1=1.5即可,由于ode23计算相对太慢，因此采用ode15计算。
源程序如下：
建立boat.m
functiondx=boat（t,x）
m=sqrt（x
（1）^2+x
（2）^2）
a=1.5;b=2
dx=[a-b*x
（1）./m;-b*x
（2）./m];
end
ts=0:
120;
x0=[0;-100];
[t,x]=ode15s（@boat,ts,x0）;
plot（x（:
1）,x（:
2））;
grid;
title（'y-x曲线'）;
xlabel（'x/m'）;ylabel（'y/m'）;
plot（t,x（:
1））;
title（'x-t曲线'）;xlabel（'t/s'）;ylabel（'x/m'）;
grid;
plot（t,x（:
2））;
title（'y-t曲线'）;xlabel（'t/s'）;ylabel（'y/m'）;...
grid;
计算结果：
与v1=0.5，1m/s时比较，从图中可以看出，小船几乎都是在[30s，35s]时间范围内到达x方向的最大值，然后逆流驶向B点。
所以v1越小，小船在逆水回到B的时间就越短。
v1=1.5时
将原程序中v1=1改为v1=2即可,由于ode23计算相对太慢，因此采用ode15计算。
源程序如下：
建立boat.m
functiondx=boat（t,x）
m=sqrt（x
（1）^2+x
（2）^2）
a=1.5;b=2
dx=[a-b*x
（1）./m;-b*x
（2）./m];
end
ts=0:
200;
x0=[0;-100];
[t,x]=ode15s（@boat,ts,x0）;
plot（x（:
1）,x（:
2））;
grid;
title（'y-x曲线'）;
xlabel（'x/m'）;ylabel（'y/m'）;
plot（t,x（:
1））;
title（'x-t曲线'）;xlabel（'t/s'）;ylabel（'x/m'）;
grid;
plot（t,x（:
2））;
title（'y-t曲线'）;xlabel（'t/s'）;ylabel（'y/m'）;
grid;
计算结果：
发现曲线的变化趋势明显不同于前几个，且最终小船无法到达B点，而到达了距离B点50米的位置。
分析：
原因在于v1=v2，则k=1，代入原式，有
，这是一个抛物线，将d=100,y=0代入解得x=50，与实验结果相符。
当到达（0,50）点后，速度矢量合成的结果便是和速度为0。
所以船会在该点处停止。
升级会员 