山东大学数学学院数学实验作业题

数学实验

成员签名：

曹云20070901005

李宪锋20070901061

李晓翻20070901062

施尚20070901110

实验二教堂顶部曲面面积的计算方法

实验题目：

教堂顶部曲面面积的计算方法

实验目的：

本试验主要涉及微积分，通过试验将复习曲面面积的计算、重积分和Taylor展

开等知识；另外将介绍重积分的数值计算法和取得函数近似解析表达式的摄动方

法。

实验内容：

思考下面这个实际问题并借助数学软件完成后面4个题的解答：

某个阿拉伯国家有一座著名的伊斯兰教堂，它以中央大厅的金色巨大拱形圆

顶名震遐迩。因年久失修，国王下令将教堂顶部重新贴金箔装饰。据档案记载，

大厅的顶部形状为球面，其半径为30m。考虑到可能的损耗和其他技术因素，实

际用量将会比教堂顶部面积多1.5%.据此，国王的财政大臣拨出了可制造

5750m有规定厚度金箔的黄金。建筑商人哈桑略通数学，他计算了一下，觉得

黄金会有盈余。于是，他以较低的承包价得到了这项装饰

工程，但在施工前的测量中，工程师发现教堂顶部实际上并非是一个精确的半球

面而是半椭圆球面，其半立轴恰是30m,而半长轴和半短轴分别是30.6m

和29.6m。这一来哈桑犯了愁，他担心黄金是否还有盈余？甚至可能短缺。最后

的结果究竟如何呢？

1.用近似格式(2.10)计算教堂顶部面积，与用格式(2.8)计算的结果

相比较；

2.试用数学软件直接计算面积(2.3)；

3.在俄国沙皇的宫廷宝藏中，有许多复活节蛋，它们大都以金银制作，装

饰着或者内藏着各种钻石。其中有一中较大的金“蛋"，"蛋''壳的外层表面是一

个椭球面，其半长轴、半短轴和半立轴分别为8cm、5.2m和5cm。“蛋”壳的

厚度为0.24cm,重量是1680go

用所学的知识解决这只复活节蛋的壳是否用纯金制作的。(金的密度是19.2g/cm)

4.建筑商人哈桑在对另一座伊斯兰建筑物顶部表面进行装饰时，他碰到的

是一个类似半球面、然而又具有一些其他变化规律的曲面，哈桑这次仍要对

该建筑物的顶部贴以金箔，我们可以确切地用球坐标表示该曲面方程，为

x-Rsin(9(l+0.1|sin6同)cos<p,

'y-Rsin^(1+0.1|sin6^>|)sinip,0<^><2^

TT

z=Reos。Q<0<—

I2

其中R=30(m),(请考虑一下，这是怎样地一个曲面？)如果由技术和损耗的

因素将使用料比实际面积多1.6%,那么装饰这个顶部至少需要多少金箔？

试用数值方法和摄动方法分别求解这个问题，并将两种方法的结果比较。

（注意：这里给出的曲面方程是参数形式的，因此首先需要弄清这种情况下曲面

的计算式有什么变化。）

采用方法：1.取椭圆中心为坐标原点建立直角坐标系，则教堂顶部半椭

圆球面的方程可写为：

z闻-，啜（2.1）

其中R=30,a=30.6,b=29.6,而其表面积为

s=-^-zjdxdy（2.2）

D

这里积分区域D为1/

LL

通过简单的计算容易得到一立+”

Ia2b2

引进变量代换x~arcos。，y-/?rsin0

2

„2->(cos0sin'

RT+

则有5=0诏f[1+------g~---------Lhrdr(2.3)

这个积分相当复杂，不过关于变量r还是可以积出初等函数的表达式，有兴趣的

读者可以试一试，若记

Jcos20sin28、/_

"=匕「+"（24）

那么（2.3）中关于r的积分

3上(1+)1-

(〃<1)

22yli-〃

（2.5）

1〃

—H----.arcsin(〃>1)

、227/7^1

这里M=1的情况要对表达式求极限。注意到M的表达式（2.4）,若将式（2.5）

带入式(2.3)得到的是一个极为复杂的积分式。事实上，这是一个无法以初等

函数形式来表达的积分，因此我们必须使用近似方法来处理它。考虑到这一积

分形式相当复杂，我们宁可直接对式(2.3)来进行处理。

2.数值积分方法：

对于二重积分，可以如同一元函数定积分那样，将区域划分为小块，然后在每个

小区域上对被积函数作近似简化求积，再把所得的值求和即可。

3.摄动方法：

简单地说，摄动方法就是对解析式中的小参数进行展开，从而求得近似解析解的

方法，应用于积分计算，常常是采取将被积函数(或其部分)展开的方法。

使用的主要程序：

程序1：

m=18;

a=8.0-0.14;

b=5.2-0.14;

R=5.0-0.14;

h=l/(2*m);

k=2*pi/(2*m);

e=0:k:2*pi;

t=(0:h:1),;

%算式(2.13)

f=sqrt(t.A2*ones(size(t))'+RA2*(1-t.A2)*((cos(e)/a).A2+(sin(e)/b).A2));

clearIij;

forj=2:2:2*m

fori=2:2:2*m

%算式(2.10)

Iij(i,j)=k*h/9*(f(i-1,j-1)+f(i+1,j-1)+f(i-1,j+1)+f(i+1,j+1)...

+4*(f(i,j-l)+f(i-1,j)+f(i+Lj)+f(i,j+1))-

+16*f(i,j));

end

end

1=sum(sum(Iij));

S=2*a*b*I;

L=0.24*S;

sprintfC不是。纯金蛋应重％7.2f克，该蛋壳密度为%5.2f(g/cm3)。

\n',19.2*L,1680/L)

程序2：

m=15;

R=30;

k=pi/6/(2*m);

h=pi/2/(2*m);

u=0:k:pi/6;

v=(0:h:pi/2),;

forj=l:2*m+l

fori=l:2*m+l

f(i,j)=sqrt(l/100*RA2*sin(v(j))A2*(101+20*sin(6*u(i))+35*cos(6*u(i))A2)...

*(l/100*RA2*cos(v(j))A2+l/5*RA2*cos(v(j))A2*sin(6*u(i))...

-l/100*RA2*cos(v(j))A2*cos(6*u(i))A2+RA2)...

-9/2500*sin(v(j))A2*RA4*cos(v(j))A2*cos(6*u(i))A2*(10+sin(6*u(i)))A2);

end

end

clearIij;

forj=2:2:2*m

fori=2:2:2*m

%算式(2.10)

Iij(i,j)=k*h/9*(f(i-l,j-l)+f(i+Lj-l)+f(i-Lj+l)+f(i+Lj+l)...

+4*(f(i,j-l)+f(i-l,j)+f(i+l,j)+f(i,j+l)).・・

+16*f(ij));

end

end

1=sum(sum(Iij));

S=12*1;

ans=sprintf(1表面积％7.2f(m2),需金箔％7.2f(m2)\n',S,(l+0.016)*S)

text(-50,-70,l,ans)

实验结果:

1.近似格式2.8计算的结果:

mSms

25621.42165679.83

45679.78245679.82

65679.89445679.81

105679.891005679.81

近似格式2.10计算的结果:

msms

25700.54165679.81

45679.88245679.81

65679.81445679.81

105679.811005679.81

2.用数学软件直接计算面积2.3得：S=5679.82

3.由算式2.10(见程序1)计算得：

不是。纯金蛋应重2004.25克，该蛋壳密度为16.09(g/cm3)。

4.由程序2计算得：

表面积为6454.59(m2),需要金箔6557.87(m2)

实验三导弹跟踪问题

实验目的:

本实验主要涉及常微分方程的建模和求解；介绍两种微分方程的数值方法：Euler

法和改进的Euler法；还介绍了仿真方法。

实验内容：

1.应用数学软件或编制计算程序对问题（3.12）-（3.14）进行数值计算，先运用

Euler法，与表3.2以及表3.3的数据比较，并以更小的步长计算结果；再用改进

的Euler法计算（步长与Euler法相同）。

2.在本实验介绍的计算过程中，我们是计算到"<",”+1之”即停止，然后取

L，小，这样做法可能会有不小的误差。有时甚至会出现整体步长改小而结果却

未必能改进的情况。由于Euler法或改进的Euler法的计算格式中每一步值的取

得仅仅依赖上一步的值，因此在计算过程中改变步长是可行的，即当计算到

为<”,"+]2”而y远大于H时，可缩小步长（例如为原来的十分之一）以xy

作为新起点继续进行迭代。试用这种变步长方法来改进在任务1中得到的结果。

3.如果当基地发射导弹的同时，敌艇立即由仪器发觉。假定敌艇为一高速快艇，

它即刻以135km/h的速度与导弹方向垂直的方向逃逸，问导弹何时何地击中快

艇？试建立数学模型并求解。

采用方法：

主要公式：

数学模型：

解析方法：

导弹轨迹方程：

l(H-y)A+,,AH

X=-[r------------------------------JH---------

24"(X+1)1-21-A2

设导弹击中敌舰于点（L,H）：

数值方法：

Euler格式:

则八为所求

匕

改进的Euler方法：

主要程序：

1.

Euler法

h=0.0005;

H=120;Vw=450;Ve=90;

clearxy;

tk=0;

k=1;

%(3.23)

x(1)=0;y(1)=0;

whiley(k)<H

%(3.21)

x(k+1)=x(k)+

Vw*h*(Ve*tk-x(k))/sqrt((Ve*tk-x(k))人2+(H-y(k))人2);

%(3.22)

y(k+1)=y(k)+Vw*h/sqrt(1+((Ve*tk-x(k))/(H-y(k)))A2);

k=k+1;tk=tk+h;

end

x

y

11

sprintf(k=%dztk=%7.4f\n,k-1,tk)

1

sprintf(L=%8.4f,T=%8.4f\n*,x(k)zx(k)/Ve)

改进的Euler法

h=0.0005;

H=120;Vw=450;Ve=90;

clearxy;

tk=0;

k=1;

%(3.28)

x(1)=0;y(1)=0;

whiley(k)<H

tkl=tk+h;

%(3.26)

xkl=x(k)+

Vw*h*(Ve*tk-x(k))/sqrt((Ve*tk-x(k))A2+(H-y(k))A2);

%(3.27)

ykl=y(k)+Vw*h/sqrt(1+((Ve*tk-x(k))/(H-y(k)))A2);

%(3.24)

x(k+1)=0.5*(xkl+x(k)+

Vw*h/sqrt(1+((H-ykl)/(Ve*tkl-xkl))人2));

%(3.25)

y(k+1)=0.5*(ykl+y(k)+

Vw*h/sqrt(1+((Ve*tkl-xkl)/(H-ykl))A2));

tk=tk+h;k=k+1;

end

x

y

sprintf(*k=%d,tk=%7.4f\n*,k-1,tk)

sprintf(1L=%8.4f,T=%8.4f\n*,x(k),x(k)/Ve)

2.

h=0.01;

H=120;Vw=450;Ve=90;

clearxy;

tk=0;

k=1;

%(3.28)

x(1)=0;y(1)=0;

whileh>0.00001

ify(k)>H

tk=tk-h;tkl=tk-h;k=k-l;h=h/10;

end

tkl=tk+h;

%(3.26)

xkl=x(k)+

Vw*h*(Ve*tk-x(k))/sqrt((Ve*tk-x(k))A2+(H-y(k))A2);

%(3.27)

ykl=y(k)+Vw*h/sqrt(1+((Ve*tk-x(k))/(H-y(k)))A2);

%(3.24)

x(k+l)=0.5*(xkl+x(k)+

Vw*h/sqrt(1+((H-ykl)/(Ve*tkl-xkl))人2));

%(3.25)

y(k+1)=0.5*(ykl+y(k)+

Vw*h/sqrt(1+((Ve*tkl-xkl)/(H-ykl))人2));

tk=tk+h;k=k+l;

end

11

sprintf(k=%dztk=%7.4f\n,k-lztk)

1

sprintf(L=%8.4f,T=%8.4f\n*,x(k)zx(k)/Ve)

3.

h=0.0005;

H=120;Vw=450;Ve=135;

elf

axis([-535-10130])

holdon

title(1pD^Oep^p-z^iiOOY1)

11

plot(0,Hzbo)

plot(0,0,1r.1)

pause

clearXwYwXeYe;

tk=0;s=0;

k=1;

Xw(1)=0;Yw(1)=0;

Xe(1)=0;Ye(1)=H;

while(Xw(k)-Xe(k))A2+(Yw(k)-Ye(k))A2>0.4

Xw(k+1)=Xw(k)+

A

Vw*h*(Xe(k)-Xw(k))/sqrt((Xe(k)-Xw(k))八2+(Ye(k)-Yw(k))2);

Yw(k+1)=Yw(k)+

Vw*h/sqrt(1+((Xe(k)-Xw(k))/(Ye(k)-Yw(k)))人2);

Xe(k+l)=Xe(k)+

Ve*h/sqrt(1+((Xe(k)-Xw(k))/(Ye(k)-Yw(k)))八2);

Ye(k+1)=Ye(k)-

Ve*h*(Xe(k)-Xw(k))/sqrt((Xe(k)-Xw(k))人2+(Ye(k)-Yw(k))^2);

s=s+sqrt((Xe(k+1)-Xe(k))A2+(Ye(k+1)-Ye(k))A2);

Wx(1)=Xw(k);

Wx(2)=Xw(k+1)

Wy(1)=Yw(k);

Wy(2)=Yw(k+1)

Ex(l)=Xe(k);

Ex(2)=Xe(k+1)

Ey(l)=Ye(k);

Ey(2)=Ye(k+1)

1

plot(Xe(k)zYe(k),*wo)

1

plot(Xw(k)zYw(k)zw.*)

plot(Xe(k+1),Ye(k+1),*bof)

f1

plot(Xw(k+1),Yw(k+1)zr.)

plot(Ex,Ey,1)

forrp=0:10

plot(Wx,Wy,1y1)

plot(Wx,Wy,1w')

plot(Wx,Wy,1r*)

end

k=k+1;tk=tk+h;

end

plot(Xe(k),Ye(k),*ro1)

plot(Xe(k),Ye(k),'y*')

text(Xe(k)-1,Ye(k)-8,!!1)

sprintf(1k=%d,tk=%7.4f\n*,k-1,tk)

ans=sprintf(1X=%8.4f,Y=%8.4f,

T=%8.4f\n,,Xe(k),Ye(k),s/Ve)

text(10,10,ans)

holdoff

pause

close

clearall

实验结果：

1.用Euler法：当令u=0.0005时,L=25.0763,T=0.2786。

用改进的Euler法：当令T=0.0005时，L=25.0608,T=0.2785。

用更小的步长，所得结果更接近解析方法的结果。

2.所得结果为：L=24.9563,T=0.2773

3.敌舰被击中的位置为(33.0906,110.2464)

实验六：个人住房抵押贷款和其他金融问题

实验题目：

个人住房抵押贷款和其他金融问题

实验目的：

本实验涉及微积分和线性代数，通过实验复习数列，函数方程求根和与线性代数方

程组有关的某些知识:主要是介绍与经济生活中某些常见重要问题有关的离散形

式数学模型-差分方程。

实验内容：

1、实际问题：

1998年12月，中国人民银行公布了新的存，贷款利率水平，其中贷款利率如表

1,表2和表3分别列出了上海商业银行报章公布的个人住房商业抵押贷款年利

率和上海商业银行提供的个人住房商业性抵押贷款（万元）还款额的部分数据。

表1

贷款期限半年1年3年5年5年以上

利率%6.126.396.667.207.56

表2

贷款期限1年2年3年4年5年

利率6.1206.2556.3906.5256.660

表3

贷款期年12345

月1224364860

月还款期到期一次还444.356305.9896237.2649196.4118

本付息

本息总额10612.0010664.5411015.6311388.7111784.71

2、数学模型:

以商业贷款1000元为例，一年期贷款年利率为6.12%,到期一次还本付息总

计10612.00元，二年期贷款的年利率为6.255%,月还款数444.3560元恰为本息

总额10664.54元的1/24,这是怎么产生的呢？

设贷款后第k个月是欠款余数位阿Ak元，月还款为m元，则由人变化到

Ak+1,有还款数和利息因素，月利率设为r,从而得到

Ak+i=(l+k)Ak-m,k=0,1,2.......(1)

开始的贷款数A。=10000(2)

成为数学模型。

月利率采用将年利率R=0.06255平均，即r=0.06255/12=0.0052125(3)

(1)称为差分方程。

3、问题的解法与讨论：

a.月还款额

二年期贷款满足Aw=0(4)

令Bk=A1Ai(5)

由(1)得及产(l+r)Bk于是得Bk=B«+r)i,k=l,2……(6)

Ai(—A()=BI+B2+...Bk

=B」l+(l+r)+...+(l+r)k-1]

k,

=(A,-A0)[[(l+r)-]/r]

=[(l+r)A°-m-Aj[[(l+r)i]/r]

从而得到差分方程的解：Ak=A°(l+r)-m/r[(l+r)i],k=0,1,2,……(7)

将Ak,A0,r,k=24的值代入得,m=444.3560(元)。

b.还款周期

如果按月还款的话，显然要比按年付款的钱少。考虑到人们的收入一般均以月

薪方式获得，因此逐月归还法对于贷款这是合适的。

若改为逐年归还方法，情况如何呢？以二年为例，则（7）中的r代为

R=0.06255,k=2,Ao=10000,得到m~=5473.8673（元）。

本息总额为2m~=10947.73（元）。

c.平衡点

若令Ag=Ak=A,可解得A=m/r称之为差分方程（1）的平衡点。

当A0=m/r时，衡有Ak=m/r,k=0,1,2,……，则每月还款额恰抵上利息，则所

欠额保持不变。

当A。稍大于或小于m/r时，Ak随着k的增大越来越远离m/r,这种平衡点称为

不稳定的。

对一般的差分方程Ak*尸f（A）k=0,1,2,……（9）

当初始值稍大于或小于差分方程的平衡点A时，若Ak-A,（当k-8）则

称A为稳定的，否则，称A为不稳定的。

4.其他问题

a.养老保险

数学模型：养老保险是与人们生活密切相关的一种保险类型。通常保险公司

会提供多种方式的养老金计划让投保人选择，在计划中详细列出保险费和养老金

的数额。例如某保险公司的一份材料指出：在每月交费200元至60岁开始领养

老金的约定下，男子若25岁起投保，届时月养老金2282元，若35岁起投保，

月养老金1056元；若45岁起投保，月养老金420元。我们来考虑这三种情况所

交保险费获得的利率。

Fk+i=Fk(l+r)+p,k=0,1,2...,N(11)

Fk+i=Fk(l+r)-q,k=N+l,N+2,....M(12)

其中P，q分别是60岁前所交的月保险费和60岁起所领月养老金数目，r是

交保险金获得利率，N,M分别是自投保起至停交保险费和至停交养老金的时间。

射手名平均值75岁。以25岁起投保为例，则有p=20均q=2282；N=420,M=600o

不难得到

kk

Fk=Fo(l+r)+p/r[(l+r)-l],k=0,1,2……N

(13)

k-NkN

Fk=FN(1+r)-q/r[(1+r)--1],k=N+l,N+2,....M

(14)

(13)中k=N,(14)中k=M

且R=0。则得关于r方程

(1+r)M-(1+q/p)(1+r)"'+q/p=O

记x=l+r,代入数据解出根x=1.00485,r=0.00485o对于35岁和45岁起投保

情况，得月利率分别为0.00461和0.00413o

b.金融公司的支付基金的流动

某金融机构为保证现金充分支付，设立一笔总额$540万元的基金，分开放置

在位于A城和B城的两公司。发现每过一周，A城公司有10%支付基金流动到B

公司，B公司则有12%支付基金流动到A城公司。此时，A城公司基金额为$260

万，B城公司基金额$280万。按此规律，两公司支付基金数额变化趋势如何？

设此后第k周结算时，A和B基金数分别是ak和bk(万元)，则有

ak+I=0.9ak+0.12bk,k=0,1,2...

bk+i=0.lak+0.88bk,k=0,1,2...(17)

3O=260,bo=28O

(18)

迭代可求得各周末时小和bk的数值。表(4)给出1至12周末两公司基金数

(单位:万美元)：

表(4)

kkk

akbk3kbkakbk

1267.6000272.40005284.5716255.42849290.8536249.1464

2273.5280266.47206284.7658253.234210291.6658248.3342

3278.1518261.84827288.4773251.522711292.2993247.7007

4281.7584258.24168289v8123250.187712292.7935247.2065

由表知，A城公司支付基金数在逐步增加，但增幅逐步减小;

B城公司支付基金数则正好相反，但ak是否有上界，bk是否有下界？(附程序)

5.问题与讨论(实验任务)

(1)确定表2中二、三、四期贷款利率是如何产生的，推导出相应的一至五年万

元贷款的还款额表。

(2)制定一张完整的个人住房商业贷款利率和还款表，贷款期从一年至二十年。

表中包含贷款期(年，月)、年利率、月利率月还款额和本息总额。

(3)小李夫妇曾准备申请商业贷款10万元用于购置住房，每月还款880.66元,

25年还清。此时，房产商介绍的一家金融机构提出：贷款10万元，每月还款440.33

元，22年还清。贷款时预付4000元。试分析两种情况哪个更合适？

(4)从还款周期比较看出逐月还款比逐年还款付出较少本息总额则逐月还款情

况如何？考虑是否有必要采取尽可能短的周期还款？

6.实验中所用的程序和方法

实验任务1

对比表6.1和6.2可知：表一中的半年利率、一年利率、三年利率分别是表二中

的一年、三年、五年利率。而表二中的利率呈线性关系。不难发现表一中的利率

亦成线性关系。表二中的二三四年利率就是根据此线性关系产生。依据前面的同

样的解差分方程的方法将前面2年期对应的数据分别改成1年、3年、4年和5

年期，分别计算ml、m3、m4、m5,将得到的数据与表3进行对比，发现数据是相

同的。这验证了表3。

实验任务2

程序如下：

Ao=10000;

sdr=[6.126.396.667.207.567.56];

sdk=[0.5135620];

fdk=1:9;

fdr=0.135*(fdk-l)+6.12;

fork=10:20

fdk(k)=k;

fdr(k)=7.20;

end

m(l)=0;

s(l)=10612.00;

r=fdr/12/100;

fork=2:20

m(k)=(A0*r(k)*(1+r(k))'(k*12))/(((1+r(k))"(k*12))-1);

s(k)=k*12*m(k);

end

tb=sprintfC贷款为d元还款表\n\n',AO);

tb=[tbsprintfC年月年利率月利率月还款额本息总额\n'n')];

fork=l:20

tb=[tbsprintf('%2d%3d%6.4f%6.4f%8.4f%8.2f\n',...

k,12*k,fdr(k),fdr(k)/12,m(k),s(k))];

end

tb

plot(sdk,sdr,fdk,fdr)

pause

close

clearall

程序结果如下表:

年月年利率月利率月还款额本息总额

1126.12000.51000.000010612.00

2246.25000.5212444.356010664.54

3366.39000.5325305.989611015.63

4486.52500.5438237.264911388.71

5606.66000.5550196.411811784.71

6726.79500.5663169.507412204.53

7846.93000.5775150.584812649.13

8967.06500.5888136.660913119.44

91087.20000.6000126.078213616.44

101207.20000.6000117.141914057.02

111327.20000.6000109.892414505.79

121147.20000.6000103.907314962.65

131567.20000.600098.894315427.52

141687.20000.600094.644515900.28

151807.20000.600091.004716380.84

161927.20000.600087.859816869.08

172047.20000.600085.122017364.88

182167.20000.600082.722817868.12

192287.20000.600080.608218378.66

202407.20000.600078.734918896.38

实验任务3函数文件：

%函数返回值rate为实际享受到的贷款年利率，％参数A。贷款总额，m每次还

款额，y贷款年数，c每年还款次数。

程序1：

A0=100000;

m=880.66;

y=25;

c=12;

rate=ep6_f2(AO,m,y,c)

k=y*c;

P(D=AO;

p(2)=-m-AO;

p(k+l)=m;

x=roots(p);

rl=(x-l)*c;

nl=1;

forn=1:k

ifrl(n)>0

r(nl)=rl(n);

nl=nl+1;

end

rate=r;

sprintf('还款年利率为%f',r);

end

程序二

A0=96000;

m=440.33;

y=22;

c=24;

rate=ep6_f2(AO,m,y,c)

k=y*c;

p(l)=AO;

p(2)=-m-AO;

p(k+l)=m;

x=roots(p);

rl=(x-l)*c;

nl=1;

forn=1:k

ifrl(n)>0

r(nl)=rl(n);

nl=nl+1;

end

rate=r;

sprintfC还款年利率为%f',r);

end

运行后1利率为0.0959,2利率为0。0969,所以并不优惠。

实验任务4

程序如下：

文件1

functionout=ep6_f1(q,k)

Ao=10000;

fdk=1:9;

fdr=0.135*(fdk-l)+6.12;

forii=10:20

fdk(ii)=ii;

fdr(ii)=7.20;

end

r=q*fdr(q)/k/100;

m=(A0*r*(l+r)"k)/((1+r)-1);

[sprintf%d年期贷款,分%d次还款:'，q,k),・・.

sprintf('每次还款%9.4f元，本息总额为%9.2f元。'，m,m*k)]

out=m*k;

文件2

s(1)=ep6_f1⑵2*1);

s(2)=ep6_f1(2,2*2);

s(3)=ep6_f1(2,2*4);

s(4)=ep6_f1(2,2*12);

s(5)=ep6_fl(2,2*12*4);

s(6)=ep6_f1(2,2*12*8);

s(7)=ep6_f1(2,2*365);

plot(s)

pause

close

程序结果为:

2年期贷款,分2次还款:每次还款5473.8673元,本息总额为10947.73元。

2年期贷款,分4次还款:每次还款2698.4778元,本息总额为10793.91元。

2年期贷款,分8次还款:每次还款1339.5528元,本息总额为10716.42元。

2年期贷款，分24次还款:每次还款444.3560元,本息总额为10664.54元。

2年期贷款，分96次还款:每次还款110.8859元,本息总额为10645.04元。

2年期贷款，分192次还款:每次还款55.4260元,本息总额为10641.79元。

2年期贷款，分730次还款:每次还款14.5745元，本息总额为10639.39元。

分析运行结果可知，综合考虑还款因素与人力因素，按月还款最为符合实际情况。

(

四、实验结果：

3、对实验任务3来说，看似第二种方法花的总钱数少，但实际上第一种情况对

小李夫妇更优惠。

4、对还款周期来说，从实验结果来看，缩短还款周期确定所付的本息总额在逐

渐减少，但减少的幅度会越来越少，最后时，几乎上下只差几元钱，这显然麻烦,

所以有时不必采取更短周期，只需按月就行。

五、实验总结：

《个人住房抵押贷款和其他金融问题》这个数学实验比较贴近现实，和我们的

日常生活息息相关，所以在做这个实验的时候我们感到比较熟悉。从这个实验可

以看出，现在的保险一般都没有在银行存款的利息高。最后，我们凭着一些生活

中的经验进行了一些简单的分析，希望老师能指出我们实验中的错误。通过互相

的配合，我们在实验中积累了不少经验，也会给我们进行以后的数学实验和数学

建模打下更好的基础。

实验九合金工厂的生产规划问题

实验题目：

合金工厂的生产规划问题

实验内容：

本实验涉及线性代数，通过实验复习矩阵的运算、初等变换和解线性方程组等知

识，另外介绍了运筹学中一类重要的问题：线性规划问题和求解的基本方法。

实际问题

某合金工厂生产甲、乙两种合金，生产每吨甲种合金需用A元素20kg、B元素

40kg和C元素90kg,而生产每吨乙种合金需用A元素100kg、B元素80kg和

C元素60kg。由于A、B、C三种元素都是原料市场上十分紧缺的货品，工厂每

月所能得到的这些元素的供应量分别为200kg、200kg、和360kg。工厂生产每

吨甲种合金的利润为30万元，生产每吨乙种合金的利润为40万元。试问：该工

厂应如何安排生产，才能获得最大利润？

实验任务

①药房有两种复合维生素制剂，甲种每粒含维生素A、B各lg,D、E各4g和

C5g,乙种每粒含维生素A3g,B2g,D1g,E3g和C2g,一顾客每天需摄入维

生素A不超过18g、B不超过13g、D不超过24g和E至少12g,问

（1）每天应服两种维生素各多少才能满足需要而且尽可能摄入较多的维生素

C?

（2）甲种复合维生素每粒1.5元，乙种复合维生素每粒1元，选择怎样的服法

才能花最少的钱而又满足每天的需要，此时顾客摄入的维生素C是多少？试用

图解法求解。

②某工厂准备制造一批设备，每件设备需要7根2m长和2根7m长的钢梁，钢

梁由150根长15m的钢材截成，问如何截取钢材能使产生的废料最少？

③试对上述两个任务用单纯形表求解。

④农场有A、B和C三块地，分别是200km2、400km2和600km2,计划种植

水稻、大豆和玉米，要求三种作物的最低收获量分别为375t、120t和750t。估计

各块地种植三种作物的单产（单位：t/km2）如表9.9所示：

表9.9

ABC

水稻11.2509.7509.000

大豆6.0006.7505.250

玉米15.00013.50012.750

应如何制订种植计划能使总产量最高？又若作物的售价为水稻1100元/t,大

豆1950元/3玉米950元/t,那么应如何制订种植计划能使总收益最高？

采用方法：

数学模型

设工厂每月生产甲种合禽兀乙种卷金利润为U万元，那么

u=30%,+40尤2（9.1）

显然，问题是要求何时有

maxu-max(30X1+40x2)(9.2)

而为,%2满足约束条件

20%+100x2<200

40X,+80X2<200（9.3）

90%+60/m360

%1>0,x2>0

其中最后一行的不等式是反映产量是非负的。

这就是实际问题的数学模型。它是求一个线性函数在非负自变量受到线性不

等式（或等式）的约束时极值问题，称之为线性规划；所求极值问题的解称为线

性规划的最优解。

处理方法

图解法

我耕在平面上考察问题（9.的%卜+鱼翔三瓦一次方程

代表了坐标平面上的一条直线，C标匚宛2秘冰噂式则代表了以

此直线为界的半平面，因此约束条件（9.3）意味着五个半平面的交集。图9.1中的

阴影部分给出了这个集合，它是一个包含边界的凸多边形OPQRS.显然，代表上

述线性规划的最优解的点必定位于这个集合中，故称此集合（确切地说是此点集

对应的坐标向量集）为线性规划的容许集。

图9.1

再来看由式（9.1）给出的目标函数M=3O玉+40々如果将U视作参数，那

么30办+4（比2=〃代表了斜率为的直线，随着U的增加或减少，直线向右上

方或左下方平移（见图9.2）容易理解：若直线经过容许集的某个顶点而此时U

再增加将使直线离开容许集，则此临界状态的直线所对应的参数值就是所求问题

的最大值，所过顶点的坐标就是问题的最优解。从图9.2可看出最优解应为R点。

图9.2

由于R是直线*+80^=200与90%+60x,=360的交点，可得最优

解为

%=3.5,羽=0.75,此时U有最大值为U=135.这说明应安排每月生产甲、乙合

金分别为3.5t、0.75t,才能获得最大利润135万元。

图解法告诉我们，线性规划的最优解往往在容许集的顶点取得。因此在确定

了容许集是有界的条件下，也可以不画图而求出容许集所有顶点，再将目标函数

在这些顶点上的值加以比较来求出最有解。这种做法一般在三个变量时采用，因

为三维空间中的多面体一般是不容易画的。

然而上述方法虽然简单，但在实际中却很少应用，因为在源于实践的真正的

线性规划问题中，变量通常不是两三个而往往是几十甚至几百个，所以既无法画

图也很难求出所有的容许集顶点，因此我们介绍了一种更普遍有效的方法。

单纯形法

单纯形法的基本思路是：线性规划(通常是求最小值的形式)若有最有解，

其必定在容许集(相应几何空间中的凸多面体)的顶点达到，故从某一个定点出

发，沿着凸多面体的棱向另一顶点迭代，使得目标函数的值下降，经过有限次迭

代，将达到最优解点。

I.标准型

首先我们通过令丫=-U,把问题变海於0石-40%,在约束下的极小

值，再引进新变量，将线性规划(9.2),(9.3)化为本：

minv=min(-30%]—40%+Ox,+0x4+Ox5)(9.4)

满足约束条件

20西+100无2+七=200

40%+80%+8=200

(9.5)

90%+60X2+X5=360

于稣工霁驻负之外，其他的约束条件都化为等式，且约束等式

的个数聚式右端项非负。我们将符合这样条件的线性规划形式

称为标准行型。

H.单纯形法

，20、'100、11)'200

L己40，=80=0,々4=1，&5=0,b=200

<90；<60,O、360

那么约束条件(9.5)中的等式组成的方程组可写为

4%+a2x2+/为+a4x4+a5x5=b(9.6)

在向量al、a2、a3、a4、a5中选取三个成为线性无关组，通常总可以(或

采取一些方法)取到单位向量，例如这里可以选取a3、a4、a5,称之为基，相应

的变量x3,x4,x5称为基变量，而其余的变量xl,x2称为非基变量。记Nl=(al,

a2),那么式(9.6)即

x.

=b-N1(9.7)

}1*2,

令

\X2j

可得⑹(200、

x=b、=b200(9.9)

^54J[360,

注意，bl的分量非负，故(9.8),(9.9)不仅构成了(9.6)的一个解，而且

也是(9.5)的一个解，称之为初始容许解。这就是单纯形法的出发点。

现考虑是否和如何转向下一个点，为说明一般的原则，把目标函数V写成

v=q由+c2x2+c3x3+c4x4+c5x5(9.10)

对应基变量J的系数为(C3,1,G)，而非基变量系数

记为

或=(。，。2)将(9.7)代入上式，则、、

V=(cl,c2)+(。3，。4，。5)八N、

7(工2JJ

/再、

=c箭一4N-4)

记(9.11)

那么丫=%一工巴与(912)

;=1

注意，对应基变量的总是为零。以上显然是将目标函数写成仅含非基变量的

形式，不过在我们的(9.4).(9.5)中目标函数已有了这形式。从而容易看出：

(1)若非基变量的系数均不大于零，那么初始容许解就是最优解。

(2)可以证明：在所有大于零的所对应的非基变量aj中，只要有一个向量的

所有分量小于或等于零，则这个线性规划无解。

在我们的问题中环=30,5=40,它们对屏矽均无非正分量，因此

不能得出以上两种结论。此时，就将进行所谓城麓

(3)换基的选择是先从(9.11)给出的所有大于零的系数中比较出最大者，则

其所对应的非基变量将成为基变量弓在速余同巡J

故心将进入基变量。那么谁将离开基变量？记非基向量勺和相应的bl分别

为

aj="4/，/=1,2;b=b，

其中的第一个下标i是与基变量的下标一致的。考做相应的比

值

b

W，i=3,4,5(9.13)

ai2

中分母为正的那些量，其中最小者的下标所对应的变量不再作为基变量。这

岬200.200.360

3100480560

最小者为“，购离开基变量，湖乂,不将成为一组新携变量。

所理应的'将被称为主元。

(4)新基变量在约束条件(9.5)中的系数向量不是单位向量，但只要利用初等

变换可以化到这样的结果。注意，要将主元所在方程中的新基变量的系数变为lo

读者容易验证：约束(9.5)可转化为

11c

-XH---尤？+X,=2

5x1100-

4…

24M—Xo+x—40

v15-d4(9.14)

3

78xj—+网=