数学建模MATLAB培训2

年份 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990
人口/.百万 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4
人口指数增长模型: x 两边取对数,得 ln x(t)
利用 a polyfit(x,
常 用 x=fzero（‘f’, [a,b]） [x，fv，ef，out]=fzero（＠f，x0，opt，p1，p2，…）
⑵、fsolve：主要用于求解非线性方程组 常 用 x=fsolve（‘f’, x0） f 须为函数名
[x，fv，ef，out，jac]=fsolve（＠f，x0，opt，p1，p2，…）
已知20世纪美国人口统计数据如下，试计算表1 中这些年 份的人口增长率。
表1 20世纪美国人口统计数据
年份 1790 1800 1810 1820 1830 1840 1850 1860 1870 1880 1890
人口/百万 3.90 5.30 7.20 9.60 12.90 17.10 23.20 31.40 38.60 50.20 62.90
二分法: 是基于连续函数的性质。设f（x）是连续函数，若对于a<b，有
f(a)·f(b)<0,则在(a,b)内f(x)至少有一个零点。（参考书王沫然P258程序erfen,收
敛速度较慢.）
牛顿切线法、割线法：xk 1
xk
f (xk ) f ( xk )
1.2 Matlab中的库函数
⑴、fzero：用于求单变量方程的根，所采用的算法主要是 二分法、割线法。
建立数学模型应考虑分有燃料和燃料已用尽两个阶段。设燃料 燃烧率为常量，燃料燃烧的时间可以确定第一阶段的时间。第二 阶段的初值由第一阶段的终值确定。引擎关闭的时刻即第一阶段 终止第二阶段开始的时刻。火箭先加速，后减速，达到最高点的 时刻火箭的运动有什么特征？
建立数学模型，明确所要求出的量
M0=1400kg，Q0=1080kg，r=18kg/s，k=0.4kg/m，N=32000N 第一阶段，燃料燃烧推动阶段，为变加速运动，初值问题为
⑶、roots：用于单变量多项式求根。 r=roots（c） c为多项式系数（按降幂排列）
输出 r 为f（x）=0的全部根（包括复根） c=poly（r） 与roots的功能相反。roots与poly互为逆函数。
例1、小张夫妇欲贷款50万元买房，他咨询了两家银行，第 一家银行开出的条件是每月还4500元，15年还清；第二家 银行开出的条件是每年还45000元，20年还清，从利率方 面看哪家银行较优惠（简单地假设年利率=月利率×12）
[tdxx] b cos , dy b sin plodtt(t,x);grid dt
即 gfitgddeuxxtrte(';x(t)')(;cgtexxtb()'(y2c(t)')(;xpa)atusey)2
pldoty(x(:,1),x(:,2)b);(garitd y)
名径存中P盘。(x,，y)保存在搜索R路(c, y)
(t) r•t
y, m)
x0 e rt ln( x0
求得 r,
x)0阻x滞(t)增长1模型( x:xm0xm
1)e
rt
x, y为实验数据 , m为拟合多项式次数 ,方法最小二乘法
解得 : r 0.2743 /10年, x0 4.1884
阻滞增长模型:
x(t)
xm
1 ( xm 1)ert
龙格-库塔方法的MATLAB实现
•
x(t
)
f (t, x)
x(t0 ) x0
其中
x (x1, x2,, xn )', f ( f1, f2,, fn )', x0 (x01, x02,, x0n )'
先对所解方程或方程组（高阶方程）编写ODE文件（M
函数文件，Ordinary Differential Equation）保存在MATLAB
result1 = 1.0e+004 * 0.0060 1.2190 result2 = 1.0e+004 * 0.0071 1.3115
0.0267 -0.0005
0.0002
14000
300
20
12000
10
250
0 10000
200 -10
8000
150
-20
6000
-30 100 4000
0.3500 12.8170 4.5552 15.0000 7.0000
0.4000 13.9806 6.1773 15.0000 8.0000
0.4500 14.7451 8.0273 15.0000 9.0000
0.5000 15.0046 9.9979 15.0000 10.0000
3、数据拟合、polyfit、人口预报.
常微分方程数值解：计算（精确）解y(x)在一系列离散点 x0
< x1 < x2 < … < xn <…上的近似值 ，通常选取相等的计算步
长 h。 例2、求解下列微分方程：
(1)
y y
' y 2 (0) 1
x '
(2)
x2 y(
y'
2
' xy'(x2
) 2, y'( )
2
0.25) y 2
0
1、非线性方程求解 ,fzero ,贷款买房利率. f (x) 0
包含三角函数、指数函数及对数函数等的超越方程与n （≥2）次代数方程一起统称为非线性方程。
如
x6 2x4 6x3 13x2 8x 12 0
x2 2x 4ln x 0
5x ex 0
1.1 求解非线性方程的方法
图形法：绘图确定根的近似位置;
y0=[h1,v1];
%第二阶段初值
[t2,y]=ode45('odehj2',[T:T+ 20],y0);
%求第二阶段解
s=max(find(y(:,2)>=0));
%速度为零时的时间分量
hm=y(s,1);
%最大高度
t=[t1;t2(2:end)];h=[x(:,1);y(2:end,1)];
%高度与时间
解、设年利率为r，贷款额A，每年还B，贷款期限n年，则
B 1 r
B (1 r)2
B (1 r)n
A
1 (1 1 ) A r (1 r)n B
(1
1 (1 r)n
)
A B
r
0
2、常微分方程数值解、ode45、火箭问题
设有常微分方程初值问题：
y' f (x, y) y( x0 ) y0
常微分方程结论: f 适当光滑，对y 满足Lipschitz条件,则方程 存在唯一解
subplot(1,3,1),plot(t,h,'b-'),axis([0 80 0 14000]) %高度函数
subplot(1,3,2),plot(t,v,'r-'),axis([0 80 -10 300])
%速度函数
subplot(1,3,3),plot(t,ac,'g-'),axis([0 80 -30 40]) %加速度函数
-40
2000
50 -50
0
0
-60
0
50
100
0
50
100
0
50
100
设 a=20,b=40,c=15;
例有先f1一u、n建艘c海t立i走o上n一私缉d船x个私=正j函is以海i(数t,一防xM)定某文速部件度缉向私正艇北上方的向雷行达驶发，现缉正私东艇方立向即C海以里最处大
速a度=2前0;往b=拦40截;c=。15用; 雷达进行跟踪时，可保持缉私艇的速度方向始 终模sd指型=xs=向，q[r(tc走确(-(xc私定(-1x船缉)()1/)。私s);^(建艇a2*+t立追(-ax*任上(2t-)意走x)(/s2时私])*)刻船^b2; 缉的)%; 私位以艇置列的，向位求量置出形和追式缉上表私的示艇时方航间程线。的数学
0.1000 3.9854 0.2924 15.0000 2.0000
0.1500 5.9445 0.6906 15.0000 3.0000
0.2000 7.8515 1.2899 15.0000 4.0000
0.2500 9.6705 2.1178 15.0000 5.0000
0.3000 11.3496 3.2005 15.0000 6.0000
模再型编设写在一t=个0时脚刻本缉M私文艇发件现走 Y
私c船lea，r 走私船以速度a平行于Y轴
正ts向=0行:0驶.05，:0缉.5私; 艇按速度b以指
向x走0=私[0船0]的'; 方向行驶(b>a)，在任
Q(c, at)
意[时t x刻]=otd缉e45私(@艇ji位si,于ts,点x0P);(x,y),则 M文件以ex4t1作为文件
O
cX
gtdextt('x');g(tecxt('yx') 2 (at y)2
设走私船在t时刻的位置为x1(t),y1(t),则
x1(t) c 15, y1(t) at 20t
时刻t 缉私艇的位置 走私船的位置
0
0
0 15.0000
0
0.0500 1.9984 0.0698 15.0000 1.0000
x0
dx rx dt x(0) x0
r r sx
dx rx(1 x )
dt
xm
x(0) x0
dx dt
Biblioteka Baidu
x
x
r(1
) xm
x(0) x0
例、汽车刹车距离. 汽车司机在行驶过程中发现前方出 现突发事件,会紧急刹车.人们把从司机决定刹车到车完全停 止这段时间内汽车行驶的距离,称为刹车距离. 现做了一组实 验：用固定牌子的汽车,由同一司机驾驶,在不变的道路、气 候条件下，对不同的车速测量其刹车距离，得到的数据如下 表所示，试建立刹车距离与车速之间的数学模型.
v=[x(:,2);y(2:end,2)];
%速度函数
L=length(h);dt=(T+12)/(L-1);
%用三点公式求加速度
ac(1)=(-3*v(1)+4*v(2)-v(3))/2/dt;
%起始点加速度
for k=2:L-1
ac(k)=(v(k+1)-v(k-1))/2/dt;
%中间点加速度
end
ac(L)=(v(L-2)-4*v(L-1)+3*v(L))/2/dt;
%最后点加速度
a1=ac(length(t1)-1);
%关闭引擎时加速度
am=ac(length(t1)+s-1);
%最高点加速度
result1=[T h1 v1 a1]
%引擎关闭时高度,速度和加速度
result2=[T+t2(s) hm am] %达到最高点时间高度和加速度
求当 x(t) 0 时的高度hm和加速度a。
huojian.m
clear
global M Q r N k
%声明全局变量
M=1400;Q=1080;N=32000;r=18;k=0.4;T=Q/r;x0=[0 0];
[t1,x]=ode45('odehj1',[0:T],x0);
%求第一阶段解
h1=x(end,1);v1=x(end,2); %引擎关闭时火箭高度和速度
设置的搜索路径，然后在解方程的脚本文件中调用：
[t, x]=ode23(@f,ts,x0,options) 或
ts=[a, b]表示自变量的 初值到终值的区间
[t, x]=ode45('f',tspan,x0,options) 或ts=[t0,t1,t2,…tf]
options用于设置误差限（缺省时相对误差为10-3，绝对误差 为10-6），用语句实现为：
车速km/h 20 刹车距离m 6.5
40
60
80 100 120 140
17.8 33.6 57.1 83.4 118.0 153.5
数据分析: 作图 plot（v，d，‘ro’） grid on 刹车距离由反应距离和制动距离两部分构成；反应距 离由反应时间和车速决定；制动距离与制动器的作用 力、车重、车速及道路、气候等因素有关
options=odeset('reltol',rt,'abstol',at)
火箭发射
小型火箭初始重量为1400kg，其中包括1080kg燃料。火箭竖直 发射时，燃料燃烧率为18kg/s，由此产生推力32000N，火箭引擎在 燃料用尽时关闭。设火箭上升时，空气阻力与速度的平方成正比， 比例系数为0.4kg/m。求引擎关闭瞬间火箭的高度、速度、加速度及 火箭到达最高点时的高度和加速度，并画出高度、速度和加速度随 时间变化的图形。
(M 0 rt)x N (M 0 rt)g kx2 x(0) x(0) 0
0 t T (T Q0 ) r
求x=T 时，hT =x(T ), vT x(T )和aT x(T )
第二阶段，燃料烧尽，火箭作减速运动，初值问题为
(M0 Q0 )x kx2 (M0 Q0 )g x(T ) hT , x(T ) vT