1第八章 其它模型

第八章 其它模型

【教学目的】：理解一些综合模型建立的思路和方法，培养解决综合问题的能力。 【教学重点难点】：

教学重点：综合应用各类数学建模方法解决较复杂的实际问题。

教学难点：针对特定的数学建模问题，如何合理使用数学方法和工具。 【课时安排】：6学时

【教学方法】：采用多媒体教学手段，配合实例教学法，通过对典型例题的讲解启发学生思维，并给与学生适当的课后思考讨论的时间，加深知识掌握的程度。安排一定课时的上机操作。 【教学内容】：

8.1 导弹跟踪问题

设位于坐标原点的甲舰向位于x轴上点A（1，0）处的乙舰发射导弹，导弹始终对准乙舰。如果乙舰以最大的速度v0沿平行于y轴的直线行驶，导弹的速度是5v0，求导弹运行的曲线。当乙舰行驶多远时，导弹将它击中？

理论知识介绍：

对于微分方程，如果能够得到解析形式的解固然助问题的分析和应用，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以得到解析解（又称通解），而绝大多数变系数方程、非线性方程都是所谓“解不出来”的，至少现行阶段求不出通解。即使看来非常简单的一个方程，例如

dy2＝y2＋x，我们是求不到其解析解的。于dx是，对于微分方程解决实际问题来说，研究其稳定性和数值解法就是一个十分重要的手段。

1）．微分方程数值解法

考虑一阶常微分方程组初值问题

y?(t)＝f(t,y(t))， t0＜t＜tf （1） y(t0)＝y0

其中y＝(y1，y2，…，ym)?，f＝(f1，f2，…，fm)?，y0 ＝(y10，y20，…，

ym0)?。所谓数值解法，就是寻求解y(t)在一系列离散节点t0＜t1＜…＜tn≤tf上的近似

值yk（k＝0，1，…，n）。称hk＝tk?1－tk为步长，通常取为常量h。最简单的数值解法是Euler法。

Euler法的思路极其简单：在节点处用差商近似代替导数

y?(tk)≈

y(tk?1)?y(tk) （2）

h这样导出计算公式（称为Euler格式）

yk?1＝yk＋hf（tk，yk）， k＝0，1，…，n （3） 它能求解各种形式的微分方程，Euler也称法折线法。之所以称为欧拉法是因为，就几何角度而言，所求得的近似解是初值问题精确解的折线逼近，而且此折线的起点是初值条件所对应的点。

由于Euler法对［tk，tk?1］上y(t)的平均斜率k＝

y(tk?1)?y(tk)仅用左端点tk的斜

h率作近似，精度只有一阶，所以较低，改进方法有二阶Runge－Kutta法、四阶Runge－Kutta法、五阶Runge－Kutta－Felhberg法和线性多步法等，这些方法可用于求解高阶常微分方程（组）初值问题。

Runge－Kutta法利用［tk，tk?1］上多个点斜率的加权平均作k的近似。从而大大提高了计算精度。四阶Runge－Kutta法格式为 yn?1＝yn＋

h（k1＋2k2＋2k3＋k4） 6k1＝f（tn，yn）

hhk2＝f（tn＋，yn＋k1） （4）

22hhk3＝f（tn＋，yn＋k2）

22k4＝f（tn＋h，yn＋hk3）

它具有四阶收敛精度。

2）．数值算法稳定性

上述提到收敛精度是单步意义上的，计算中还涉及到误差的传播问题。不适当的步长可能导致计算误差恶性发展而使计算失败，这种现象称为数值不稳定。考虑一阶方程

y?＝?y， y(0)＝1

?t)恒正且渐近稳定。现取步长h用Euler法（3）式进其中?＜0，从而它的解y(t)＝exp(行求解

yk?1＝（1＋?h）yk

当h＜－1/?时，数值解正确反映了解函数性质；当－1/?＜h＜2/?时，数值解仍稳定但发生振荡；当h＞2/?时，数值解不稳定且振荡趋于无穷。一般地，对一阶方程组 y?＝Ay＋b

Euler法数值稳定性条件为｜1＋?h｜＜1，其中?为矩阵A的任意任意特征值。四阶Runge－Kutta法稳定性条件稍宽但基本上处于同一数量级。

3）．常微分方程的解的稳定性

考虑一阶常微分方程初值问题（1），称（1）式的特解?(t)是稳定的，若对??＞0存在?＞0，使当｜y0－?(t0)｜＜?时，以y0为初值的解满足｜y(t)－?(t)｜＜?（t＞；若还有｜y(t)－?(t)｜→0（t→∞），称?(t)是渐近稳定的。稳定性表明初值y0的小t0）

扰动不会对解产生明显影响。微分方程描述的运动的轭线密切依赖于初值，而初值的计算或测定实际上不可避免地出现误差和干扰，如果描述运动的微分方程的特解是不稳定的，初值的微小误差或干扰，将招致“差之毫厘，记得谬以千里”的严重后果，因此不稳定的特解不宜作为设计的依据；反之，稳定的特解都才是我们最感兴趣的。

4）．求微分方程（组）数值解的MATLAB命令

ode23 2/3阶Runge－Kutta法； ode15s 刚性方程组的解法；

ode45 4/5阶Runge－Kutta法； ode23s 刚性方程组的解法； [T,Y] = solver( 'odefun ' , tspan, y0, options )；

solver取上述四个函数之一，不同的函数代表不同的内部算法，ode45是最常用的求解微分方程命令，它采用变步长四阶Runge－Kutta法和五阶Runge－Kutta－Felhberg法，ode23与ode45类似，只是精度低一些。对于刚性方程，宜采用ode15s或ode23s；

字符串odefun是用以表示待解方程写成的M－文件名；

tspan = [ t0 tf ]表示自变量初值t0和终值tf，y0表示函数初始向量；

options用于设定求解的误差限（该参数可以缺省，缺省时相对误差为103，绝对误差－

为106）；

输出列向量T表示节点（t0, t1, ... , tn）T，输出矩阵Y表示微分方程（组）数值解； 可以在MATLAB帮助文件中查阅有关命令的更多信息。 问题求解： 1．问题的分析

设任意时刻t，乙舰的坐标为（X(t)，Y(t)），导弹的坐标为（x(t)，y(t)）。

导弹速度恒为v，从原点射出，且速度是横向和纵向距离的一阶导数的矢量和，则有

－

?dy??dx?2 ??＋??＝?v? （5）

?dt??dt? 由于导弹头始终对准乙舰，帮弹头的速度向量平行于乙舰与导弹头位臵的差向量，即

22dx＝?（X－x） ?＞0 （6） dtdy ＝?（Y－y）

dt 故由（5）、（6）可得

(?（X－x）)2＋(?（Y－y）)2＝?v? （7）

2由（7）式可解得

?＝

v22 （8）

(X?x)?(Y?y)由题目已知乙舰以最大速度v0沿直线x＝1运动，为便于分析，我们可暂且令v0＝1，则导弹速度v＝5v0＝5，乙舰的坐标可重新表示为

X(t)＝1 （9） Y(t)＝v0t＝t

t＝0时，甲舰位于原点，则导弹坐标的初始值可确定

x(0)＝0 （10） y(0)＝0

因此，由（6）、（8）、（9）、（10）式可知导弹运动轭迹的参数方程为

dx5＝(1?x)

22dt(1?x)?(t?y)dy5 ＝(t?y)

22dt(1?x)?(t?y) 2．MATLAB计算机求解

建立M文件eq2.m如下

function dy=eq2(t,y)

dy=zeros(2,1);% 初始化dy的两个分量 dy(1)=5*(1-y(1))/sqrt((1-y(1))^2+(t-y(2))^2); dy(2)=5*(t-y(2))/sqrt((1-y(1))^2+(t-y(2))^2); 在MATLAB命令框中输入命令

>> t0=0;tf=2; % 确定自变量的初始值和终值 >> [t,y]=ode45('eq2',[t0,tf],[0,0]);

>> plot(y(:,1),y(:,2),'-')% 画出导弹运动轭迹图

>> x=0:0.01:1.5; >> y=0.2:0.2;

>> hold on% 在同一个图形座标上继续作图 >> plot(x,y,'-')% 画一条经过y值为0.2的横线

由上述的代码可得到导弹的运动轨迹图1，由图1我们可以看出，导弹大致在（1，0.21）处击中乙舰。

21.81.61.41.210.80.60.40.2000.511.5

图1 导弹运动轭迹图

3．结果分析：

对于引例问题，我们用常微分方程的数值解法得到了一个近似解，实际上该问题可以参照实验四交通管理问题中的初等积分方法求得导弹运动轭迹的解析解为

555 y＝ －(1?x)5＋ (1?x)5＋

24128其分析过程与我们上面讨论的类似，求解也不难，读者不妨自己动手做一下。

4655，即当乙舰航行到（1，）处时被甲舰发射的导弹击中，被击2424y5中的时间为t＝＝。若取v0＝1，则在t＝0.21时被击中。这与我们通过数值方法

v024v0当x＝1时，y＝求解的结果非常近似。

8.2 个人住房抵押贷款

随着经济的发展,金融问题正越来越多地进入普通市民的生活，贷款、保险、养老金和信用卡等都涉及金融问题，个人住房抵押贷款是其中最重要的一项。1998年12月,中国人民银行公布了新的存、贷款利率水平,其中贷款利率如下表所列:

表1 中国人民银行贷款利率表

贷款期限 半年 利率﹪ 6.12

一年 三年 6.39

6.66

五年 7.20

五年以上 7.56

当贷款期处于表中所列相邻年限之间时利率为对应相邻两数中较大者。其后,上海商业银行对个人住房商业性贷款利率作出相应调整。表2和表3分别列出了上海市个人住房商业抵押贷款年利率和商业抵押贷款(万元)还款额的部分数据(仅列出了五年)。

表2 上海市商业银行住房抵押贷款利率表

贷款期限 一年 利率﹪ 6.12

表3 上海市商业银行住房抵押贷款分期付款表

贷款期限 一年

二年

三年 305.99 11015.63

四年 237.26 11388.71

五年 196.41 11784.71

月还款（元） 到期一次还清 444.36

10664.54 本息总和（元） 10612.00

一个自然的问题是,表2和表3是如何依据中央人民银行公布的存、贷款利率水平制定

的？

我们以商业贷款10000元为例,一年期贷款的年利率为6.12﹪,到期一次还本付息总计10612.00元,这很容易理解。然而二年期贷款的年利率为6.225﹪,月还款数444.36元为本息和的二十四分之一,这后两个数字究竟是怎样产生的？是根据本息总额算出月还款额，还是恰好相反？让我们稍微仔细一些来进行分析。由于贷款是逐月归还的,就有必要考察每个月欠款余额的情况。

设贷款后第k个月时欠款余额为Ak元，月还款m元,则由Ak变化到Ak+1,除了还款额外,还有什么因素呢？无疑就是利息。但时间仅过了一个月,当然应该是月利率,设为r,从而得到

Ak?1?Ak=r Ak-m 或者

Ak+1=(1+r)Ak-m , k=1，2… (8.1)

初时条件

A0=10000 (8.2)

这就是问题的数学模型。其中月利率采用将年利率R=0.06255平均。即

r=0.06255/12=0.00512125 (8.3)

若m是已知的,则由(8.1)式可以求出Ak中的每一项,我们称(8.1)式为一阶差分方程。 模型解法与讨论 （1）月还款额

二年期的贷款在24个月时还清,即

A24=0 (8.4)

为求m的值,设

Bk?Ak?Ak?1，k=1,2,… (8.5) 易见

Bk?1?(1?r)Bk

二年 6.225

三年 6.390

四年 6.525

五年 6.660

于是导出Bk的表达式

Bk=B1(1+r)k-1, k=1,2,… (8.6) 由(8.5)式与(8.6)式得

Ak-A0=B1+B2+…+Bk=B1[1+(1+r)+…+(1+r)k-1]

(1?r)k?1(1?r)k?1 =(A1-A0)[]=[(1+r)A0-m-A0][ []

rr

从而得到差分方程(8.1)的解为

Ak=A0(1+r)k-m[(1+r)k-1]/r, k=1,2,… (8.7) 将A24,A0,r的值和k=24代入得到m=444.36(元),与表8.3中的数据完全一致,这样我们就了解了还款额的确定方法。

当然还款额表的制定依赖于年利率表,而后者又是怎样制定的呢？尽管我们无法获知银行方面的各种考虑,但还是可以通过比较分析得出一些结论。首先注意表8.2商业性贷款利率中有两个数据与中央银行公布的表8.1中的数据相同,不过相应的贷款年限则放宽了一档:6.12﹪是一年期,而在表8.1中是上一档半年期,6.66﹪是五年期而在表8.1中是上一档三年期。其次再考察表8.2商业性贷款二、三、四年期的利率,我们把这三个数字是如何得到的问题留给读者。

依据这两个结论,请读者自己制定出住房商业性贷款直至二十年的利率表和还款额表。 (2)还款周期

我们看到个人住房贷款是采用逐月归还的方法,虽然依据的最初利率是年利率。那么如 果采用逐年归还的方法,情况又如何呢？仍然以二年期贷款为例,显然,只要对(8.7)式中的利率r代之以年利率R=0.06255,那么由k=2,A2=0,A0=10000,则可以求出年还款额应为

~ =5473.87(元) m这样本息和总额为

~=10947.73(元) 2m远远超出逐月还款的本息总额。考虑到人们的收入一般都以月薪方式获得,因此逐月归还对于贷款者来说是比较合适的。读者还可以讨论缩短贷款周期对于贷款本息总额的影响。 （3）平衡点

回到差分方程(8.1),若令Ak+1=AK=A,可解出

A=m/r (8.8) 称之为差分方程的平衡点或称之为不动点。显然,当初值A0=m/r时,将恒有Ak=m/r,k=0,1,…。

在住房贷款的例子里,平衡点意味着如果贷款月利率r和月还款额m是固定的,则当初贷款额稍大于或小于m/r时,从方程(8.1)的解的表达式(8.7)中容易看出，欠款额Ak随着k的增加越来越远离m/r,这种情况下的平衡点称为不稳定的,对一般的差分方程

Ak+1=f(Ak), k=0,1,2，… (8.9)

当初始值稍大于或小于平衡点的值A时,若

limAk?A, (8.10)

k???则称A为稳定的,否则称A为不稳定的。判别平衡点A是否稳定的一个方法是考察导数f'(A): 当f'(A)<1时,A是稳定的; 当f'(A)>1时,A是不稳定的。

8.3 曲柄滑块机构的运动规律

曲柄滑块机构是一种常用的机械结构,它将曲柄的转动转化为滑块在直线上的往复运动,是气压机、冲床、活塞式水泵等机械的主机构。图1 为其示意图。

记曲柄OQ的长为r,连杆QP的长为l, 当曲柄绕固定点O以角速度w旋转时, 由连杆带动滑块P在水平槽内做往复直线运动。假设初始时刻曲柄的端点Q位于水平线段OP上, 曲柄从初始位臵起转动的角度为?,而连杆QP与OP的锐夹角为?(称为摆角) 。在机械设计中要研究滑块的运动规律和摆角的变化规律, 确切的说,要研究滑块的位移,速度和加速度关于?的函数关系,摆角?及其角速度和角加速度关于?的函数关系, 进而 (1) 求出滑块的行程s(即滑块往复运动时左、右极限位臵间的距离);

(2) 求出滑块的最大和最小加速度(绝对值), 以了解滑块在水平方向上的作用力; (3) 求出?的最大和最小角加速度(绝对值), 以了解连杆的转动惯量对滑块的影响; 在求解上述问题时,我们假定:

r?100(mm),l?3r?300(mm),??240(转/min)

（求解作为练习。）

8.4 教堂顶部曲面面积计算

某个阿拉伯国家有一座著名的伊斯兰教堂，它以中央大厅的金色巨大拱形圆顶名震遐迩。因年久失修，国王下令将教堂顶部重新贴金箔装饰。据档案记载，大厅的顶部形状为半球面，其半径为30m。考虑到可能的损耗和其他技术因素，实际用量将会比教堂顶部面积多1.5％.据此， 国王的财政大臣拨出了可制造 5750m 有规定厚度金箔的黄金。 建筑商人哈桑略通数学，他计算了一下，觉得黄金会有盈余。于是，他以较低的承包价得到了这项装饰 工程，但在施工前的测量中，工程师发现教堂顶部实际上并非是一个精确的半球面而是半椭圆球面， 其半立轴恰是 30 m ， 而半长轴和半短轴分别是30.6m和29.6m。

取椭圆中心为坐标原点建立直角坐标系，则教堂顶部半椭圆球面的方程可写为 x2y2z?R1?2?2?2.1? ab

其中R=30，a=30.6 ，b=29.6，而其表面积为 22s?1?z?zdxdyxy

D

x2y2??1这里积分区域D为

a2b2

通过简单的计算容易得到

R2x2R2y2 ?4422 1?zx?zy?1?a2b2xy 1?2?2ab

引进变量代换

x?arcos?,y?brsin?

则教堂顶部曲面面积为

2 sin2??22?cos??Rr??22? 2x1ab???abrdrS?d?1? 001?r2

????

（1） 利用数值积分方法，用梯形法和simpson法两种近似格式计算教堂顶部曲面面积； （2） 利用摄动的方法近似计算教堂顶部曲面面积； （3） 试用数学软件直接计算教堂顶部曲面面积。

理论

（1）数值积分方法

对于二重积分，可以如同一元函数定积分那样，将区域划分为小块，然后在每个小区域上对被积函数作近似简化求积 ，再把所得的值求和即可。

考虑矩形区域D上的二重积分，将 D 划分 mn 个相等的小矩形 ，s和t分别是s和t方向的分点，那么小矩形上的积分可写为

tjsi dtf?s,t?dstj?1si?1记 si?i?t??f?s,t?ds si?1则 tjI??i?t?dtij tj?1若对这两个单积分都用梯形法 ，就有

k

?i?t???f?si,t??f?si?1,t??2

而

kh

Iij?fsi?1,tj?1?fsi,tj?1?fsi?1,tj?fsi,tj4

mn这样便可求得在D上的积分 I 的近似值 I?Iiji?1j?1

cd? , 而记 当将分点增加一倍使得 k ? , h 2m2n

si?ik,tj?jh(i?1,2,,2m;j?1,2,,2n)

那么对 tj?1si?1I?dtf(s,t)ds(i?1,3,L,2m?1;j?1,3,L,2n?1)ij tj?1si?1的两次积分都用 Simpson 法，就得到

si?1k

f(s,t)ds?[f(si?1,t)?4f(si,t)?f(si?1,t)]si?13

khIij?{f(si?1,tj?1)?f(si?1,tj?1)?f(si?1,tj?1)从而

9

?f(si?1,tj?1)?4[f(si,tj?1)?f(si?1,tj)?

f(si?1,tj)?f(si,tj?1)]?16f(si,tj)}

于是有积分 I 的近似值为

2m?12n?1

I?Iij,i,j取奇数 i?1j?1（2）摄动方法

简单地说，摄动方法就是对解析式中的小参数进行展开，从而求得近似解析解的方法，应用于积分计算，常常是采取将被积函数（或其部分）展开的方法，通过一个简单例子来说明这方法。

? 对于 ? (0,1), 计算

? d?2?????????????????????I(?)??01??sin2?

利用Taylor公式，关于参数 展开，有

1 ?12?(1??sin?)2 21??sin? 11?31?3?5?1?(??sin2?)?2(??sin2?)2?3(??sin2?)3

222!23!

1?3?5?724n(2n?1)!!?(??sin?)??(?1)(??sin2?)n 4n24!2n! n?o(?)

sin??1余项的写法考虑到了

1? (2n?1)!!n2n 我们称级数 (? 1) n ( ? ? sin ? ) 为函数 1 ? sin 2 ? 的渐近级数，通常 ?2n!n?0

1应用渐进级数的有限项来近似函数。例如，在这里把前n+1 项替换 代入

1??sin2?

???(2n?1)!!d?2sin2n?d??) ?积分式 I ( ? 2 (2.14) 那么由于 002(2n)!!1??sin2?

可得

??19225312254[(2n?1)!!]2n? I(?)??1??????????n??2?464256163842n!(2n)!!?

当 ? ? 0.6 时，用上式前5项计算 的误差不超过0.01.

???I(?)

【思考与讨论】

按揭贷款的还款时间与还款期数对总还款额的影响。

【作业布臵】：

调查当地的房产价格和银行利率，并算出面积为100平方米的商品房在首付3成的情况下办理银行贷款进行15年和20年按揭所需还款总额和月还款金额。

本文来源：https://www.bwwdw.com/article/57u2.html