东南大学《数学实验》报告
学号 姓名 成绩
实验内容：曲线拟合与插值
一 实验目的
用最小二乘法实现多项式拟合；3次样条函数的应用
二 预备知识
(1)熟悉一般的曲线拟合的最小二乘法原则
(2)熟悉正规方程、差分表、均差表的概念
(3)熟悉“\”、polyfit、polyval、interp1、spline等Matlab命令
三 实验内容与要求
下表给出了氨蒸汽的一组温度和压力数据。那么能否从所列的数据中计算75
氨蒸汽的压力？
表氨蒸汽的压力和温度
（1）用polyfit、polyval命令求解氨蒸汽问题
（2）用解正规方程的方法令求解氨蒸汽问题
（3）已知某平原地区的一条公路经过如下坐标点，请用不同的插值方法绘出这条公路（不考虑公路的宽度）。对于上表给出的数据，编程计算三次样条插值估计的公路长度。
第二篇：东南大学数学建模与实验+实验报告
数学建模与实验
实验报告
授课教师
计算机科学与工程学院
目录
实验1—3.4节“企业利润合理使用”例题的求解………………………………1实验2—Hill密码加密、解密……………………………………………………2实验3—习题5.3“样条差值法绘制公路”求解…………………………………3实验4—Volterra方程组求解（改进欧拉公式与龙格-库塔公式比较）……5实验5—习题6.8“饮酒驾车的药物注射模型”求解……………………………7实验6—银行贷款利息的计算……………………………………………………9
实验1——3.4节“企业利润合理使用”例题的求解
实验目的
运用幂法求解矩阵的最大特征值、特征向量、一致性指标和随机一致性比率。
实验原理
1、判断矩阵A只有一个最大特征根λmax，且λmax=yk+1
。2、可以使用迭代法求解：
（1）任取一个初始向量x0≠0（2）u0=x0
k=1，2，…
vk=Auk?1
mk=max{vk}uvkk=
mk
代码实现
%Eigen.m函数定义
function[w,m,CI,CR]=Eigen(AC)
n=numel(AC)^0.5;%获得矩阵行、列数w=AC(:,1:1);%初始向量x0
flag=0;whileflag~=n
flag=0;m=0;%λ
wc=w;%上一轮计算的ukv=AC*w;%vk=Auk-1
fork=1:n
m=m+v(k);endw=v/m;
%uk=vk/mk
l=abs(wc-w)./w;%uk变化速度(相对误差)
fork=1:n
ifl(k)<=0.0001
flag=flag+1;
%flag=n说明所有元素都趋于稳定endendend
CI=(m-n)/(n-1);
yk
3、一致性指标和随机一致性比率
CI=
λmax?nn?1CR=
CIRIswitchn
%RI表
case1
RI=0;case2
RI=0;case3
RI=0.58;
……case11
RI=1.51;otherwise
RI=2;
endCR=CI/RI;
%执行文件AC=[1,1/5,1/3;5,1,3;3,1/3,1];[w,lmd,CI,CR]=Eigen(AC)C1P=[1,3;1/3,1];
[w,lmd,CI,CR]=Eigen(C1P)C2P=[1,1/5;5,1];
[w,lmd,CI,CR]=Eigen(C2P)C3P=[1,2;1/2,1];
[w,lmd,CI,CR]=Eigen(C3P)
执行结果
A—C判断矩阵：w=
0.10470.63700.2583lmd=
3.0385CI=
0.0193CR=
0.0332
—P判断矩阵：C1C1—w=
0.75000.2500lmd=
2CI=0CR=
NaN
—P判断矩阵：C2C2—w=
0.16670.8333lmd=
2CI=0CR=
NaN
—P判断矩阵：C3C3—w=
0.66670.3333lmd=
2CI=0CR=
NaN
实验2——Hill密码加密、解密
实验目的
利用Hill2加密原理将特定文本加密，并将其解密。
实验原理
加密：将原文本按字母值表转化为数字，并将此数字阵列α两两配对与加密矩阵相乘Aα得到密码阵列β。再次按字母值表反向翻译即可得到密文。
解密：首先计算出解密矩阵A-1(mod26)，按与加密相同的方法将密文转化为数字阵列β，与解密矩阵相乘即可得到原文数字阵列α=A-1β(mod26)，反向翻译可得到原文。
加密矩阵要求模26可逆。
代码实现
%原文=xxxxxxxxxxxxxxxxxxxxCode='';
text='xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'%原文Answer='';
textSize=size(text);ifmod(textSize(2),2)==1
text=[text,text(textSize(2))];end%如果原文长度为奇则末尾补位A=[1,8;0,9]
%加密矩阵
%字母表
A=mod(inv(A)*27,26)
%模26逆矩阵
codeArray=mod(abs(Code)-64,26);%字母表codeSize=size(codeArray);fori=1:2:codeSize(2)-1
b=[codeArray(i),codeArray(i+1)];
a=mod(A*b',26);%模26ifa(1)==0
a(1)=26;endifa(2)==0
a(2)=26;end
a1=setstr(a(1)+64);a2=setstr(a(2)+64);Answer=[Answer,a1,a2];endAnswer
%加密过程%模26%转为ASCII码
%Z转化为26以便输出
%转为ASCII码
textArray=mod(abs(text)-64,26);textSize=size(text);fori=1:2:textSize(2)-1
b=[textArray(i),textArray(i+1)];a=mod(A*b',26);c1=setstr(a(1)+64);c2=setstr(a(2)+64);Code=[Code,c1,c2];end
Code
执行结果
原文：（xxxxxxxxxxxxxxxxxxxxx）text=
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
加密矩阵：A=
10
89
解密矩阵：A=
10
23
密文：Code=
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
解密结果：Answer=
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
实验3——习题5.3“样条差值法绘制公路”求解
实验目的
利用样条差值法，根据已知坐标，绘制整条曲线。
实验原理
估测公路函数满足三次样条差值条件。
公路在(478，296)处折返，因而整条曲线不是函数曲线，故将公路在折点处分段。由于g(x)未知，根据表中数值估测g’(x)的值，其中
g’(x0+)=g’(0+)=-16/30g’(x28-)=g’(478-)=16/48g’(x28+)=g’(478+)=-12/30g’(x38-)=g’(200-)=-10/40估测公路长度时，以x轴的一米为间隔微分公路函数，求其总和。
代码实现
%全部采样点
X=[0,30,50,70,80,90,120,148,170,180,202,212,230,248,268,271,280,290,300,312,320,340,360,372,382,390,416,430,478,440,420,380,360,340,320,314,280,240,200];
Y=[80,64,47,42,48,66,80,120,121,138,160,182,200,208,212,210,200,196,188,186,200,184,188,200,202,240,246,280,296,308,334,328,334,346,356,360,392,390,400];%低值第一段采样点n1=29;
X1=[0,30,50,70,80,90,120,148,170,180,202,212,230,248,268,271,280,290,300,312,320,340,360,372,382,390,416,430,478];
Y1=[80,64,47,42,48,66,80,120,121,138,160,182,200,208,212,210,200,196,188,186,200,184,188,200,202,240,246,280,296];Z1=0:1:478;m1=length(Z1);gbar1=[-8/15,1/3];lmd1(1)=1;mu1(n1)=1;
%估算g’(x)
end
M1=inv(A1)*d1';
%A1*M1=d1
EndA1(1,1)=2;A1(1,2)=1;A1(n1,n1-1)=1;A1(n1,n1)=2;fori=2:n1-1
A1(i,i-1)=mu1(i);A1(i,i)=2;
A1(i,i+1)=lmd1(i);end
d1(1)=6*((Y1(2)-Y1(1))/h1(1)-gbar1(1))/h1(1);d1(n1)=6*(gbar1(2)-(Y1(n1)-Y1(n1-1))/h1(n1-1))/h1(n1-1);fori=2:n1-1
lmd1(i)=h1(i)/(h1(i-1)+h1(i));mu1(i)=1-lmd1(i);
d1(i)=6*((Y1(i+1)-Y1(i))/h1(i)-(Y1(i)-Y1(i-1))/h1(i-1))/(h1(i-1)+h1(i));
%计算hj,μj,λj,dj
fori=1:n1-1
h1(i)=X1(i+1)-X1(i);
fork=1:m1
fori=1:n1-1
ifZ1(k)>=X1(i)&Z1(k)<=X1(i+1)
S1(k)=M1(i)*(X1(i+1)-Z1(k))^3/(6*h1(i))+M1(i+1)*(Z1(k)-X1(i))^3/(6*h1(i))+(Y1(i)-M1(i)*h1(i)^2/6)*(X1(i+1)-Z1(k))/h1(i)+(Y1(i+1)-M1(i+1)*h1(i)^2/6)*(Z1(k)-X1(i))/h1(i);
breakendendendn2=11;
X2=[200,240,280,314,320,340,360,380,420,440,478];
Y2=[400,390,392,360,356,346,334,328,334,308,296];
Z2=200:1:478;m2=length(Z2);gbar2=[-1/4,-6/19];lmd2(1)=1;mu2(n2)=1;fori=1:n2-1
h2(i)=X2(i+1)-X2(i);end
d2(1)=6*((Y2(2)-Y2(1))/h2(1)-gbar2(1))/h2(1);d2(n2)=6*(gbar2(2)-(Y2(n2)-Y2(n2-1))/h2(n2-1))/h2(n2-1);fori=2:n2-1
lmd2(i)=h2(i)/(h2(i-1)+h2(i));mu2(i)=1-lmd2(i);
d2(i)=6*((Y2(i+1)-Y2(i))/h2(i)-(Y2(i)-Y2(i-1))/h2(i-1))/(h2(i-1)+h2(i));endA2(1,1)=2;A2(1,2)=1;
%计算hj,μj,λj,dj
%估算g’(x)
%获得S1(x)各点值
%高值第二段采样点
A2(n2,n2-1)=1;A2(n2,n2)=2;fori=2:n2-1
A2(i,i-1)=mu2(i);A2(i,i)=2;
A2(i,i+1)=lmd2(i);end
M2=inv(A2)*d2';fork=1:m2
fori=1:n2-1
ifZ2(k)>=X2(i)&Z2(k)<=X2(i+1)
S2(k)=M2(i)*(X2(i+1)-Z2(k))^3/(6*h2(i))+M2(i+1)*(Z2(k)-X2(i))^3/(6*h2(i))+(Y2(i)-M2(i)*h2(i)^2/6)*(X2(i+1)-Z2(k))/h2(i)+(Y2(i+1)-M2(i+1)*h2(i)^2/6)*(Z2(k)-X2(i))/h2(i);
breakendendend
%获得S1(x)各点值
%绘图
plot(Z1,S1,Z2,S2,X,Y,'o')%估算公路长度L=0;fort=1:477
L=L+((Z1(t)-Z1(t+1))^2+(S1(t)-S1(t+1))^2)^0.5;end
fort=1:277
L=L+((Z2(t)-Z2(t+1))^2+(S2(t)-S2(t+1))^2)^0.5;end
%A2*M2=d2
执行结果
公路绘制：（见右图）公路长度约为1016.3米L=
1.0163e+003
实验4——Volterra方程组求解(改进欧拉公式与龙格-库塔公式比较)
实验目的
使用“改进的欧拉公式”和“4阶龙格-库塔公式”分别对Volterra方程求解，绘制解曲线、相轨线，并将结果进行比较。
实验原理
改进的欧拉公式局部截断误差为O(h3)，而4阶龙格-库塔公式截断误差达到O(h5)。在本例中将两公式推广到解微分方程组。
代码实现
%Volterra.m函数定义function[as]=Volterra(t,x)as(1)=x(1)*(1-0.1*x(2));as(2)=x(2)*(-0.5+0.02*x(1));%执行文件h=0.1;
t=0:h:15;
%改进的欧拉公式x1=[25,2];fori=1:150
k1=Volterra(t(i),x1(i:i,:));
k2=Volterra(t(i+1),x1(i:i,:)+h*k1);x1(i+1:i+1,:)=x1(i:i,:)+h*(k1+k2)/2;end
%4阶龙格-库塔方法x2=[25,2];fori=1:150
k1=Volterra(t(i),x2(i:i,:));
k2=Volterra(t(i)+h/2,x2(i:i,:)+h*k1/2);k3=Volterra(t(i)+h/2,x2(i:i,:)+h*k2/2);k4=Volterra(t(i)+h,x2(i:i,:)+h*k3);x2(i+1:i+1,:)=x2(i:i,:)+h*(k1+2*k2+2*k3+k4)/6;
end
%绘制解曲线和相轨线subplot(2,2,1)plot(t,x1)subplot(2,2,2)plot(t,x2)subplot(2,2,3)
plot(x1(1:110,1:1),x1(1:110,2:2))subplot(2,2,4)
plot(x2(1:110,1:1),x2(1:110,2:2))
执行结果
曲线绘制：
改进的欧拉公式的解曲线、相轨线4阶龙格-库塔方法的解曲线、相轨线
计算结果比较：
t
00.10.20.30.40.50.60.70.80.911.11.21.31.41.51.61.71.81.922.12.22.32.42.52.62.72.82.933.1……1010.110.210.310.410.510.610.710.8
X1改进的欧拉公式25.000027.080029.330731.762934.387537.214640.253443.511846.995050.705154.639658.789463.137167.653872.295977.001181.683386.227890.485894.271497.360799.4982100.412899.847197.601593.587687.877280.726072.553863.878255.220547.0189……13.703314.817316.028917.346218.777720.332522.020323.851025.8351
X2
4阶龙格-库塔方法25.000027.081829.334531.769034.396037.225740.267343.528347.014150.726554.662758.813263.160267.674172.310577.006281.673886.197390.426994.175597.219499.3046100.163999.546597.262093.231587.533180.422472.311163.703355.106846.9513……13.714714.831116.045517.365818.800720.359422.051423.886925.8762
2.00002.00402.01672.03902.07202.11702.17572.25022.34292.45702.59652.76612.97193.22143.52403.89124.33714.87905.53726.33527.29868.45339.821611.416613.235015.248417.398419.596621.733623.696425.387126.7391……2.20122.15432.11342.07842.04972.02732.01172.00332.0025
X1改进的欧拉公式
2.00002.00412.01702.03942.07262.11782.17682.25152.34462.45912.59912.76932.97593.22633.53013.89874.34634.89025.55076.35107.31668.47289.840911.433213.245015.247317.381919.562121.681323.629725.311926.6618……2.19772.15102.11032.07552.04702.02482.00952.00132.0008
X2
4阶龙格-库塔方法
实验5——习题6.8“饮酒驾车的药物注射模型”求解
实验目的
运用药物注射模型，使用曲线拟合方法，解释饮酒驾车的一些实际问题。
实验原理
由于酒精不需要进入肠道即可被吸收，且胃对其吸收速率也非常快，本题应采用“快速静脉注射模型”。酒精主要存在于血液中，故本例应计算吸收室的血药浓度c1(t)=A1e-αt+B1e-βt
相关系数可以通过拟合法求解。
代码实现、执行结果及分析
formatshortg
%题中提供的某人喝了两瓶啤酒后血液酒精浓度随时间变化表
t=[0.25;0.5;0.75;1;1.5;2;2.5;3;3.5;4;4.5;5;6;7;8;9;10;11;12;13;14;15;16];c=[30;68;75;82;84;77;70;68;58;51;50;41;38;35;28;25;18;15;12;10;7;7;4];%根据此变化表拟合求解相关系数
ft=fittype('A1*exp(-a*x)+B1*exp(-b*x)');
options=fitoptions('Method','NonlinearLeastSquares');
options.StartPoint=[0-100000];
cfit=fit(t,c,ft,options);
plot(cfit,t,c,'o');
A1=cfit.A1
B1=cfit.B1
a=cfit.a
b=cfit.b
由此解得：(拟合曲线见右图)
A1=110.55
B1=-151.46
a=0.17949
b=2.8243
%---1---%
%问题：某人中午12点喝了一瓶啤酒，下午6点检
查合格，晚饭又喝一瓶，次日凌晨2点检查未通过，请对此情况做出解释。
t11_1=6;
c11_1=(A1*exp(-a*t11_1)+B1*exp(-b*t11_1))/2
%下午6点酒精浓度除以2是因为此人只喝了一瓶，浓度减半
t11_2=13.2;
c11_2=(A1*exp(-a*t11_2)+B1*exp(-b*t11_2))/2
t12=7.2;
c12=(A1*exp(-a*t12)+B1*exp(-b*t12))/2
c12A=c12+c11_2
由此解得：
c11_1=18.829
c11_2=5.1712
c12=15.181
c12A=20.352%次日凌晨2点总残留酒精浓度>20，检测不合格%下午6点酒精浓度<20，故检测合格%晚饭喝酒在次日凌晨2点残留酒精浓度%次日凌晨2点总残留酒精浓度%中午喝酒在次日凌晨2点残留酒精浓度
本例中实际计算的时间其实未到达次日凌晨2点，按照这组数据得到次日凌晨2点的酒精浓度将低于20，原因在于题目中提供数据的饮酒者耗散酒精的能力可能比本例中的人要强。数据仅仅体现一种趋势，实际数字因人而异。他第二次检测不合格的根本原因在于中午摄入的酒精仍有残余，提高了整体酒精浓度。
%---2---%
%问题：短时间内喝啤酒3瓶多长时间之后才能驾车？t2=0.2:0.1:24;
fori=1:239
c2(i)=(A1*exp(-a*t2(i))+B1*exp(-
b*t2(i)))*3/2;
end
plot(t2,c2,t2,20);
从右图绘制曲线的交点可以看出，喝啤酒3瓶后约12小时之后血液酒精浓度才会低于20，符合驾车标准。%---3---%
%问题：怎样估计血液中的酒精含量在什么时候最高？通过前面的两张图我们已经可以看出酒精含量的变化趋势，血液中酒精含量大约在1.5小时之后达到最高。%---4---%
%问题：根据模型论证，如果天天喝酒，是否还能开车？t4=0.2:0.1:72;
fori=1:length(t4)
c41(i)=(A1*exp(-a*t4(i))+B1*exp(-b*t4(i)))*1.5;end
fori=121:length(t4)
c42(i)=(A1*exp(-a*t4(i-120))+B1*exp(-b*t4(i-120)))*1.5;
end
fori=241:length(t4)
c43(i)=(A1*exp(-a*t4(i-240))+B1*exp(-b*t4(i-240)))*1.5;
end
fori=361:length(t4)
c44(i)=(A1*exp(-a*t4(i-360))+B1*exp(-b*t4(i-360)))*1.5;
end
fori=481:length(t4)
c45(i)=(A1*exp(-a*t4(i-480))+B1*exp(-b*t4(i-480)))*1.5;
end
fori=1:length(t4)
c4A1(i)=c41(i)+c42(i)+c43(i)+c44(i)+c45(i);%每隔12小时饮啤酒3瓶
c4A2(i)=c41(i)+c43(i)+c45(i);
%每隔24小时饮啤酒3瓶
end
plot(t4,c4A1,t4,20);
plot(t4,c4A2,t4,20);
由右图上可以看到，如果某人每隔12小时饮啤酒3瓶，他的血液酒精浓度将永远在20以上，不能驾车。由右图
下则看到，如果某人每隔24小时饮啤酒3瓶，每天在饮酒12后可以开车。
实验6——银行贷款利息的计算
实验目的
运用个人住房抵押贷款模型，解决一些实际问题。
实验原理
设贷款后第k个月时欠款余额为Ak，月还款m元，月利率r，则Ak+1?Ak=rAk?m。差分方程的解为Ak=A0(1+r)?
k
m
(1+r)k?1,k=1,2,?r
t3=3*12;
A(t3)=A0*(1+r)^t3-m*((1+r)^t3-1)/r;A(t3)由此解得：
ans=255269.41
即三年后欠款额为￥255269.41。
%---4---%
%问题：在问题1的情况下，两年后提前还款2万元，月还款额不变，则还款期限缩短多长时间？t4=2*12;
A(t4)=A0*(1+r)^t4-m*((1+r)^t4-1)/r;A(t4)=A(t4)-20000;
t_shot=eval(solve('m2=(A(t4)*(1+r)^t*r)/((1+r)^t-1)','t'))+t4;dt=t-ceil(t_shot)由此解得：
dt=40.00，即还款期限缩短了40个月。
%---5---%
%，%问题：在问题1的情况下，两年后年利率调为5.495.49%月还款额应调整到多少？t5=2*12;
A(t5)=A0*(1+r)^t5-m*((1+r)^t5-1)/r;ry5=0.0549;r5=ry5/12;
m5=(A(t5)*(1+r5)^(t-t5)*r5)/((1+r5)^(t-t5)-1)由此解得：
m5=2432.83，即月还款额应调整到￥2432.83。
[]
代码实现、执行结果及分析
formatbank
A0=300000;
%每月还款，等额本息，贷款30万，年利率4.59%%---1---%
%问题：若贷款15年，每月应还款多少？ry=0.0459;r=ry/12;t=15*12;
m=(A0*(1+r)^t*r)/((1+r)^t-1)由此解得：
m=2308.80，即每月还款￥2308.80。
%---2---%
%问题：若还款能力为每月￥2600-￥2700，贷款期限应为多长？m2=2600;
t_max=eval(solve('m2=(A0*(1+r)^t*r)/((1+r)^t-1)','t'))m2=2700;
t_min=eval(solve('m2=(A0*(1+r)^t*r)/((1+r)^t-1)','t'))由此解得：t_max=152.51t_min=144.95
即贷款期限应在153个月到145个月之间。%---3---%
%问题：在问题1的情况下，三年后欠款额？