数学建模与实验实验报告

数学建模与试验试验汇报讲课教师 计算机科学与工程学院计算机科学与技术系oxstar目录试验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只有一种最大特性根，且。 2、可以使用迭代法求解： （1）任取一种初始向量 （2） k=1，2，… 3、一致性指标和随机一致性比率代码实现%Eigen.m函数定义function[w,m,CI,CR]=Eigen(AC)n=numel(AC)^0.5; %获得矩阵行、列数w=AC(:,1:1); %初始向量x0flag=0;whileflag~=nflag=0;m=0; %λwc=w; %上一轮计算的ukv=AC*w; %vk=Auk-1fork=1:nm=m+v(k);endw=v/m; %uk=vk/mkl=abs(wc-w)./w; %uk变化速度(相对误差)fork=1:nifl(k)<=0.0001flag=flag+1;%flag=n阐明所有元素都趋于稳定endendendCI=(m-n)/(n-1);switchn %RI表case1RI=0;case2RI=0;case3RI=0.58;……case11RI=1.51;otherwiseRI=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.0332C1—P判断矩阵：w=0.75000.2500lmd=2CI=0CR=NaNC2—P判断矩阵：w=0.16670.8333lmd=2CI=0CR=NaNC3—P判断矩阵：w=0.66670.3333lmd=2CI=0CR=NaN试验2——Hill密码加密、解密试验目的运用Hill2加密原理将特定文本加密，并将其解密。试验原理 加密：将原文本按字母值表转化为数字，并将此数字阵列α两两配对与加密矩阵相乘Aα得到密码阵列β。再次按字母值表反向翻译即可得到密文。 解密：首先计算出解密矩阵A-1(mod26)，按与加密相似的措施将密文转化为数字阵列β，与解密矩阵相乘即可得到原文数字阵列α=A-1β(mod26)，反向翻译可得到原文。 加密矩阵规定模26可逆。代码实现%原文=我是来自计算机学院三班的学生牛星Code='';text='WOSHILAIZIJISUANJIXUEYUANSANBANDEXUESHENG??????' %原文Answer='';textSize=size(text);ifmod(textSize(2),2)==1text=[text,text(textSize(2))];end %假如原文长度为奇则末尾补位A=[1,8;0,9] %加密矩阵textArray=mod(abs(text)-64,26); %字母表textSize=size(text);fori=1:2:textSize(2)-1b=[textArray(i),textArray(i+1)];a=mod(A*b',26); %模26c1=setstr(a(1)+64); %转为ASCII码c2=setstr(a(2)+64);Code=[Code,c1,c2];end %加密过程CodeA=mod(inv(A)*27,26) %模26逆矩阵codeArray=mod(abs(Code)-64,26); %字母表codeSize=size(codeArray);fori=1:2:codeSize(2)-1b=[codeArray(i),codeArray(i+1)];a=mod(A*b',26); %模26ifa(1)==0a(1)=26;endifa(2)==0a(2)=26;end %Z转化为26以便输出a1=setstr(a(1)+64); %转为ASCII码a2=setstr(a(2)+64);Answer=[Answer,a1,a2];endAnswer执行成果原文：（我是来自计算机学院三班的学生牛星）text=WOSHILAIZIJISUANJIXUEYUANSANBANDEXUESHENG??????加密矩阵：A=1809密文：Code=MEETADUCTCDCEGIVDCJGWQCIJOIVJITJOHISETMVOV??????解密矩阵：A=1203解密成果：Answer=WOSHILAIZIJISUANJIXUEYUANSANBANDEXUESHENG??????试验3——习题5.3“样条差值法绘制公路”求解试验目的运用样条差值法，根据已知坐标，绘制整条曲线。试验原理 估测公路函数满足三次样条差值条件。 公路在(478，296)处折返，因而整条曲线不是函数曲线，故将公路在折点处分段。由于g(x)未知，根据表中数值估测g’(x)的值，其中 g’(x0+)=g’(0+)=-16/30 g’(x28-)=g’(478-)=16/48 g’(x28+)=g’(478+)=-12/30 g’(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]; %估算g’(x)lmd1(1)=1;mu1(n1)=1;fori=1:n1-1h1(i)=X1(i+1)-X1(i);endd1(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-1lmd1(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));End %计算hj,μj,λj,djA1(1,1)=2;A1(1,2)=1;A1(n1,n1-1)=1;A1(n1,n1)=2;fori=2:n1-1A1(i,i-1)=mu1(i);A1(i,i)=2;A1(i,i+1)=lmd1(i);endM1=inv(A1)*d1'; %A1*M1=d1fork=1:m1fori=1:n1-1ifZ1(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);breakendendend %获得S1(x)各点值%高值第二段采样点n2=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]; %估算g’(x)lmd2(1)=1;mu2(n2)=1;fori=1:n2-1h2(i)=X2(i+1)-X2(i);endd2(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-1lmd2(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));end %计算hj,μj,λj,djA2(1,1)=2;A2(1,2)=1;A2(n2,n2-1)=1;A2(n2,n2)=2;fori=2:n2-1A2(i,i-1)=mu2(i);A2(i,i)=2;A2(i,i+1)=lmd2(i);endM2=inv(A2)*d2'; %A2*M2=d2fork=1:m2fori=1:n2-1ifZ2(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:477L=L+((Z1(t)-Z1(t+1))^2+(S1(t)-S1(t+1))^2)^0.5;endfort=1:277L=L+((Z2(t)-Z2(t+1))^2+(S2(t)-S2(t+1))^2)^0.5;end执行成果公路绘制：（见右图）公路长度约为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:150k1=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:150k1=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阶龙格-库塔措施的解曲线、相轨线计算成果比较：tX1X2X1X2改善的欧拉公式4阶龙格-库塔措施改善的欧拉公式4阶龙格-库塔措施025.000025.00002.00002.00000.127.080027.08182.00402.00410.229.330729.33452.01672.01700.331.762931.76902.03902.03940.434.387534.39602.07202.07260.537.214637.22572.11702.11780.640.253440.26732.17572.17680.743.511843.52832.25022.25150.846.995047.01412.34292.34460.950.705150.72652.45702.4591154.639654.66272.59652.59911.158.789458.81322.76612.76931.263.137163.16022.97192.97591.367.653867.67413.22143.22631.472.295972.31053.52403.53011.577.001177.00623.89123.89871.681.683381.67384.33714.34631.786.227886.19734.87904.89021.890.485890.42695.53725.55071.994.271494.17556.33526.3510297.360797.21947.29867.31662.199.498299.30468.45338.47282.2100.4128100.16399.82169.84092.399.847199.546511.416611.43322.497.601597.262013.235013.24502.593.587693.231515.248415.24732.687.877287.533117.398417.38192.780.726080.422419.596619.56212.872.553872.311121.733621.68132.963.878263.703323.696423.6297355.220555.106825.387125.31193.147.018946.951326.739126.6618…………1013.703313.71472.2.197710.114.817314.83112.15432.151010.216.028916.04552.11342.110310.317.346217.36582.07842.075510.418.777718.80072.04972.047010.520.332520.35942.02732.024810.622.020322.05142.01172.009510.723.851023.88692.00332.001310.825.835125.87622.00252.0008试验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.A1B1=cfit.B1a=cfit.ab=cfit.b由此解得：(拟合曲线见右图)A1=110.55B1=-151.46a=0.17949b=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 %中午喝酒在次日凌晨2点残留酒精浓度t12=7.2;c12=(A1*exp(-a*t12)+B1*exp(-b*t12))/2 %晚饭喝酒在次日凌晨2点残留酒精浓度c12A=c12+c11_2 %次日凌晨2点总残留酒精浓度由此解得：c11_1=18.829 %下午6点酒精浓度<20，故检测合格c11_2=5.1712c12=15.181c12A=20.352 %次日凌晨2点总残留酒精浓度>20，检测不合格本例中实际计算的时间其实未抵达次日凌晨2点，按照这组数据得到次日凌晨2点的酒精浓度将低于20，原因在于题目中提供数据的饮酒者耗散酒精的能力也许比本例中的人要强。数据仅仅体现一种趋势，实际数字因人而异。他第二次检测不合格的主线原因在于中午摄入的酒精仍有残存，提高了整体酒精浓度。%---2---%%问题：短时间内喝啤酒3瓶多长时间之后才能驾车？t2=0.2:0.1:24;fori=1:239c2(i)=(A1*exp(-a*t2(i))+B1*exp(-b*t2(i)))*3/2;endplot(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;endfori=121:length(t4)c42(i)=(A1*exp(-a*t4(i-120))+B1*exp(-b*t4(i-120)))*1.5;endfori=241:length(t4)c43(i)=(A1*exp(-a*t4(i-240))+B1*exp(-b*t4(i-240)))*1.5;endfori=361:length(t4)c44(i)=(A1*exp(-a*t4(i-360))+B1*exp(-b*t4(i-360)))*1.5;endfori=481:length(t4)c45(i)=(A1*exp(-a*t4(i-480))+B1*exp(-b*t4(i-480)))*1.5;endfori=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瓶endplot(t4,c4A1,t4,20);plot(t4,c4A2,