数学实验题库 建立模型并写出求解模型的Matlab代码或程序

Matlab 编程和模拟模型1. 编程训练题目一. 有一个向量()n a a a v ,,,21 =，编程求解下列式子： 2232121121u ni i i a a a s a s +++==∑=- ， 其中u 为不超过n 的最大奇数。

自己给出测试例子，并给出运算结果。

二. 模拟一次扑克牌洗牌过程，输出牌的顺序。

%4*13+2paihao=[11:23,31:43,51:63,71:83,98,99];if length(paihao)~=54,error('error in number')endPAIXU=randperm(54);fapai = paihao(PAIXU)三. 模拟1000次投币过程，输出出现正面的频率。

程序1：num_zheng=0;for i=1:1000,cur=rand;if cur<0.5,num_zheng=num_zheng + 1;endendf = num_zheng/1000四.问题：现有一个文本文件（mydata.txt），包含了火车的出发和到达时刻，基本格式如下（实际中该文件超过1000行）19:35 07:4719:14 07:1219:21 07:1918:53 08:2319:07 07:0519:00 06:5821:37 08:0420:40 07:4021:44 07:2221:30 07:2519:28 06:4822:00 07:1019:21 07:1918:53 08:2319:07 07:05……要求把时间读取并存为矩阵数据，并且需要将其以0:00时刻记为初始时刻，把奇数行数据转为分钟。

准备知识：（1）文件打开（fopen），关闭（fclose）；（2）读取各行数据（fgetl）（3）解析各行数据，并转换数据注：处理方法不唯一。

文本文件：fopen，fclose，fgetl。

%先初始化矩阵mat=[];row = 0;fid=fopen('mydata.txt');while 1tline = fgetl(fid);if feof(fid), break, enddisp(tline)cur=str2num(tline);row = row + 1;%将当前行的(行)向量保存到mat中if isempty(mat),mat=cur;elselen=length(cur);mat(row,1:len)=cur;endendfclose(fid);19 35 7 47 8 31 9 4219 14 7 12 0 0 0 019 21 7 19 0 0 0 018 53 8 23 0 0 0 019 7 7 5 0 0 0 019 0 6 58 0 0 0 021 37 8 4 0 0 0 020 40 7 40 0 0 0 021 44 7 22 0 0 0 021 30 7 25 0 0 0 019 28 6 48 0 0 0 0五.（适中）有一个文件mytest.txt，内容如下：请编程将这个文件中从第2行开始的所有数保存为1个列向量（按行先后顺序排列这些数）程序%先初始化矩阵vec=[];tline = fgetl(fid);%跳过第1行字符串，不处理fid=fopen('mytest.txt');while 1tline = fgetl(fid);if feof(fid), break, enddisp(tline)cur=str2num(tline);%将当前行的(行)向量保存到vec 中vec=[vec,cur];endfclose(fid);vec=vec'%请将此文件中的数值读取为1个向量1 4 72 5 83 6 9vec =147258369六. （较难）有一个函数235243211x x x x x x y ++++=，另有一个函数)(x f 定义如下。

matlab在数学建模中的应用——第三章程序代码2第一篇：matlab在数学建模中的应用——第三章程序代码2 clearsyms a b;c=[a b]';A=[174179 183 189 207 234 220.5 256270 285];B=cumsum(A);% 原始数据累加n=length(A);for i=1:(n-1)C(i)=(B(i)+B(i+1))/2;% 生成累加矩阵end% 计算待定参数的值D=A;D(1)=[];D=D';E=[-C;ones(1,n-1)];c=inv(E*E')*E*D;c=c';a=c(1);b=c(2);% 预测后续数据F=[];F(1)=A(1);for i=2:(n+10)F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a;endG=[];G(1)=A(1);for i=2:(n+10)G(i)=F(i)-F(i-1);%得到预测出来的数据endt1=1995:2004;t2=1995:2014;G, a, b % 输出预测值，发展系数和灰色作用量plot(t1,A,'o',t2,G)%原始数据与预测数据的比较第二篇：matlab在数学建模中的应用——第三章程序代码1 clearsyms a b;c=[a b]';A=[89677,99215,109655,120333,135823,159878,182321,209 407,246619,300670];B=cumsum(A);% 原始数据累加n=length(A);for i=1:(n-1)C(i)=(B(i)+B(i+1))/2;% 生成累加矩阵end% 计算待定参数的值D=A;D(1)=[];D=D';E=[-C;ones(1,n-1)];c=inv(E*E')*E*D;c=c';a=c(1);b=c(2);% 预测后续数据F=[];F(1)=A(1);for i=2:(n+10)F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a;endG=[];G(1)=A(1);for i=2:(n+10)G(i)=F(i)-F(i-1);%得到预测出来的数据endt1=1999:2008;t2=1999:2018;Gplot(t1,A,'o',t2,G)%原始数据与预测数据的比较xlabel('年份')ylabel('利润')第三篇：matlab在数学建模中的应用——第一章程序代码clearclc% 读入人口数据（1971－2000年）Y=[338************8345******345093452***345***93452**** **]% 读入时间变量数据（t＝年份－1970）T=[*********2627282930]% 线性化处理for t = 1:30,x(t)=exp(-t);y(t)=1/Y(t);end% 计算，并输出回归系数Bc=zeros(30,1)+1;X=[c,x'];B=inv(X'*X)*X'*y'for i=1:30,% 计算回归拟合值z(i)=B(1,1)+B(2,1)*x(i);% 计算离差s(i)=y(i)-sum(y)/30;% 计算误差w(i)=z(i)-y(i);end% 计算离差平方和SS=s*s';% 回归误差平方和QQ=w*w';% 计算回归平方和UU=S-Q;% 计算，并输出F检验值F=28*U/Q% 计算非线性回归模型的拟合值for j=1:30,Y(j)=1/(B(1,1)+B(2,1)*exp(-j));end% 输出非线性回归模型的拟合曲线（Logisic曲线）plot(T,Y)第四篇：2)线性代数在数学建模中的应用例举8015985.docAct3 总复习【Arrangement】1）模拟题2）线性代数在数学建模中的应用例举3）线性代数在考研中的地位和重要性【Content】模拟题一、填空题（每题4分，共20分）：1、n阶方阵A的行列式，则行列式。

数学实验题库建⽴模型并写出求解模型的Matlab代码或程序数学实验题库实验1 Matlab概述12实验2 函数图形绘图 3实验3 数列极限与函数极限 2 实验4 导数与偏导数的计算 2 实验5 ⽅程近似解的求法 3实验6 定积分的近似计算 312实验7 多元函数的极值问题 31．某化⼯⼚⽣产A 、B 、C 、D 四种产品，每种产品⽣产1吨消耗⼯时和产值如下：要求全⼚年产值为1000万元以上，建⽴使⽣产消耗总⼯时最⼩的数学模型，并求解. 解：设⽣产A 产品1x 吨、B 产品 2x 吨、C 产品3x 吨、D 产品4x 吨，则所⽤⼯时为123410030040075x x x x +++，产值为12345100.5x x x x +++线性规划模型为：1234min 10030040075z x x x x =+++ 1234..5100.510000s t x x x x +++≥MATLAB 代码为：clear;c=[100;300;400;75]; A=[1 5 10 0.5]*(-1); b=[10000]*(-1); Aeq=[]; beq=[]; beq0=[]; lb=0*c;ub=[inf;inf;inf;inf]; digits(5);[x,fval]=linprog(c,A,b,Aeq,beq,lb,ub)2．贝尔⾦属公司要⽣产两种灯，制造⼀盏中国海灯需要耗费黄铜2磅和3个铣床⼩时，⽽3制造⼀盏马坦扎斯海湾灯需要耗费黄铜4磅和1个铣床⼩时，另外每盏中国海灯需要2⼈特制的东⽅灯罩，这种灯罩必须从⾹港进⼝，⽬前每个⽣产周期，由于联邦法的限制，只能进⼝100个。

且下⼀周期公司的黄铜供应量限制为320磅，铣床时间限制为180⼩时，⽽每盏中国海灯的利润为60美元，每盏马坦扎斯海湾灯的利润为30美元，为得到最⼤利润，贝尔公司应该如何安排⽣产？建⽴使利润最⼤的数学模型，并求解. 解：设⽣产中国海灯1x 盏、马坦扎斯海湾灯 2x 盏，则利润为126030x x + 线性规划模型为：12max 6030z x x =+12121212243203180..201000,0x x x x s t x x x x +≤??+≤? ?+?≤??>>?MATLAB 代码为：clear;c=[-60;-30];A=[2 4;3 1; 2 0];b=[320;180;100]; Aeq=[];beq=[];lb=[0;0];ub=[inf;inf];[x,fval]=linprog(c,A,b,Aeq,beq,lb,ub)3．伯恩公司⽣产铝制品的煎锅和焙盘，每个煎锅或焙盘都需要10盎司的铝。

matlab数学建模常用模型及编程摘要：一、引言二、MATLAB 数学建模的基本概念1.矩阵的转置2.矩阵的旋转3.矩阵的左右翻转4.矩阵的上下翻转5.矩阵的逆三、MATLAB 数学建模的常用函数1.绘图函数2.坐标轴边界3.沿曲线绘制误差条4.在图形窗口中保留当前图形5.创建线条对象四、MATLAB 数学建模的实例1.牛顿第二定律2.第一级火箭模型五、结论正文：一、引言数学建模是一种将现实世界中的问题抽象成数学问题，然后通过数学方法来求解的过程。

在数学建模中，MATLAB 作为一种强大的数学软件，被广泛应用于各种数学问题的求解和模拟。

本文将介绍MATLAB 数学建模中的常用模型及编程方法。

二、MATLAB 数学建模的基本概念在使用MATLAB 进行数学建模之前，我们需要了解一些基本的概念，如矩阵的转置、旋转、左右翻转、上下翻转以及矩阵的逆等。

1.矩阵的转置矩阵的转置是指将矩阵的一行和一列互换，得到一个新的矩阵。

矩阵的转置运算符是单撇号(’)。

2.矩阵的旋转利用函数rot90(a,k) 将矩阵a 旋转90 的k 倍，当k 为1 时可省略。

3.矩阵的左右翻转对矩阵实施左右翻转是将原矩阵的第一列和最后一列调换，第二列和倒数第二列调换，依次类推。

matlab 对矩阵a 实施左右翻转的函数是fliplr(a)。

4.矩阵的上下翻转matlab 对矩阵a 实施上下翻转的函数是flipud(a)。

5.矩阵的逆对于一个方阵a，如果存在一个与其同阶的方阵b，使得：a·bb·a=|a|·|b|·I，则称矩阵b 是矩阵a 的逆矩阵。

其中，|a|表示矩阵a 的行列式，I 是单位矩阵。

在MATLAB 中，我们可以使用函数inv(a) 来求解矩阵a 的逆矩阵。

三、MATLAB 数学建模的常用函数在MATLAB 数学建模过程中，我们经常需要使用一些绘图和数据处理函数，如绘图函数、坐标轴边界、沿曲线绘制误差条、在图形窗口中保留当前图形、创建线条对象等。

matlab数学建模程序代码
当进行数学建模时，MATLAB是一个强大的工具，用于实现和测试模型。

下面是一个简单的MATLAB代码示例，演示如何使用MATLAB进行一维线性回归建模：
```matlab
%生成示例数据
x=[1,2,3,4,5];
y=[2.8,3.9,4.8,5.5,6.3];
%进行一维线性回归
coefficients=polyfit(x,y,1);
slope=coefficients(1);
intercept=coefficients(2);
%绘制原始数据和回归线
scatter(x,y,'o','DisplayName','原始数据');
hold on;
plot(x,polyval(coefficients,x),'r-','DisplayName','回归线');
hold off;
%添加标签和图例
xlabel('X轴');
ylabel('Y轴');
title('一维线性回归建模示例');
legend('show');
%输出回归方程的系数
fprintf('回归方程：y=%.2fx+%.2f\n',slope,intercept);
```
此代码生成了一些示例数据，然后使用一维线性回归对数据进行建模。

回归方程的系数将被计算，并且原始数据与回归线将在图上显示。

请注意，这只是一个简单的示例，实际上，你可能需要根据你的具体问题修改代码。

2、已知速度曲线v(t) 上的四个数据点下表所示t=[0.15,0.16,0.17,0.18];v=[3.5,1.5,2.5,2.8];x=0.15:0.001:0.18y=i n t e r p1(t,v,x,'s p l i n e')S=t r a p z(x,y)p=p o l y f i t(x,y,5);d p=p o l y de r(p);d p x=p o l y v a l(d p,0.18)运行结果S=0.0687Dpx=-3、计算图片文件tu.bmp 给出的两个圆A，B 的圆心，和两个圆的两条外公切线和两条内公切线的切点的坐标。

（1）计算A 圆的圆心坐标I=imread('tu.bmp');[m,n]=size(I)BW=im2bw(I)BW(:,200:512)=1;figure, imshow(BW)ed=edge(BW);[y,x]=find(ed);x0=mean(x), y0=mean(y)r1=max(x)-min(x),r2=max(y)-min(y)r=(r1+r2)/4x0 =109.7516y0 =86.7495r1 =162r2 =158r =80（2）B圆的圆心坐标和半径I=imread('tu.bmp');BW=im2bw(I)BW(:,1:200)=1;imshow(BW)ed=edge(BW);[y,x]=find(ed);x0=mean(x), y0=mean(y)r1=max(x)-min(x),r2=max(y)-min(y)r=(r1+r2)/4x0 =334.0943y0 =245.7547r1 =165r2 =158 r = 80.7500外公切线上的切点f=@(x)[(x(1,1)-109.7516)^2+(x(1,2)-86.7495)^2-80.5^2(x(2,1)-334.0943)^2+(x(2,2)-245.7547)^2-80.75^2(x(2,2)-x(1,2))*(x(1,2)-86.7495)+(x(2,1)-x(1,1))*(x(1,1)-109.7516)(x(2,2)-x(1,2))*(x(2,2)-245.7547)+(x(2,1)-x(1,1))*(x(2,1)-334.0943)(x(1,1)-x(2,1))^2+(x(1,2)-x(2,2))^2+0.75^2-(334.0943-109.7516)^2-(245.7 516-86.7495)^2];xy1=fsolve(f,rand(2,2))xy2=fsolve(f,100*rand(2,2))xlswrite('book1.xls',xy1)xlswrite('book1.xls',xy2,'Sheet1','A4')xy1 =156.2419 21.0312380.7270 179.8309xy2 =153.7425 48.4651289.4819 284.38084、求微分方程组的数值解，并画出解曲线dy=@(t,y)[-10*y(1)+10*y(2);28*y(1)-y(2)-y(1)*y(3);-8/3*y(3)+y(1)*y(2)]; [t,y]=ode45(dy,[0,10],[1;0;0])subplot(3,1,1),plot(t,y(:,1),'*')subplot(3,1,2),plot(t,y(:,2),'*')subplot(3,1,3),plot(t,y(:,3),'*')0123456789105、预测2012-2020年美国人口数量。

Matlab 数学实验报告一、实验目的通过以下四组实验,熟悉MATLAB的编程技巧,学会运用MATLAB的一些主要功能、命令，通过建立数学模型解决理论或实际问题。

了解诸如分岔、混沌等概念、学会建立Malthu模型和Logistic 模型、懂得最小二乘法、线性规划等基本思想。

二、实验内容2.1实验题目一2.1.1实验问题Feigenbaum曾对超越函数y=λsin(πｘ)（λ为非负实数）进行了分岔与混沌的研究，试进行迭代格式x k+1=λsin(πx k),做出相应的Feigenbaum图2.1.2程序设计clear;clf;axis([0,4,0,4]);hold onfor r=0:0.3:3.9x=[0.1];for i=2:150x(i)=r*sin(3.14*x(i-1));endpause(0.5)for i=101:150plot(r,x(i),'k.');endtext(r-0.1,max(x(101:150))+0.05,['\it{r}=',num2str(r)]) end加密迭代后clear;clf;axis([0,4,0,4]);hold onfor r=0:0.005:3.9x=[0.1];for i=2:150x(i)=r*sin(3.14*x(i-1));endpause(0.1)for i=101:150plot(r,x(i),'k.');endend运行后得到Feigenbaum图2.2实验题目二2.2.1实验问题某农夫有一个半径10米的圆形牛栏，长满了草。

他要将一头牛拴在牛栏边界的桩栏上，但只让牛吃到一半草，问拴牛鼻子的绳子应为多长？2.2.2问题分析如图所示，E为圆ABD的圆心，AB为拴牛的绳子，圆ABD为草场，区域ABCD为牛能到达的区域。

问题要求区域ABCD等于圆ABC的一半，可以设BC等于x,只要求出∠a和∠b就能求出所求面积。

二，hamiton回路算法提供一种求解最优哈密尔顿的算法---三边交换调整法，要求在运行jiaohuan3（三交换法）之前，给定邻接矩阵C和节点个数N，结果路径存放于R中。

bianquan.m文件给出了一个参数实例，可在命令窗口中输入bianquan,得到邻接矩阵C和节点个数N以及一个任意给出的路径R,，回车后再输入jiaohuan3,得到了最优解。

由于没有经过大量的实验，又是近似算法，对于网络比较复杂的情况，可以尝试多运行几次jiaohuan3,看是否能到进一步的优化结果。

%%%%%%bianquan.m%%%%%%%N=13;for i=1:Nfor j=1:NC(i,j)=inf;endendfor i=1:NC(i,i)=0;endC(1,2)=6.0;C(1,13)=12.9;C(2,3)=5.9;C(2,4)=10.3;C(3,4)=12.2;C(3,5)=17.6;C(4,13)=8.8;C(4,7)=7.4;C(4,5)=11.5;C(5,2)=17.6;C(5,6)=8.2;C(6,9)=14.9;C(6,7)=20.3;C(7,9)=19.0;C(7,8)=7.3;C(8,9)=8.1;C(8,13)=9.2;C(9,10)=10.3;C(10,11)=7.7;C(11,12)=7.2;C(12,13)=7.9;for i=1:Nfor j=1:Nif C(i,j) < infC(j,i)=C(i,j);endendendfor i=1:NC(i,i)=0;endR=[4 7 6 5 3 2 1 13 12 11 10 9 8];<pre name="code" class="plain">%%%%%%%%jiaohuan3.m%%%%%%%%%%n=0;for I=1:(N-2)for J=(I+1):(N-1)for K=(J+1):Nn=n+1;Z(n,:)=[I J K];endendendR=1:Nfor m=1:(N*(N-1)*(N-2)/6)I=Z(m,1);J=Z(m,2);K=Z(m,3); r=R;if J-I~=1&K-J~=1&K-I~=N-1 for q=1:(J-I)r(I+q)=R(J+1-q);endfor q=1:(K-J)r(J+q)=R(K+1-q);endendif J-I==1&K-J==1r(K)=R(J);r(J)=R(K);endif J-I==1&K-J~=1&K-I~=N-1 for q=1:(K-J)r(I+q)=R(I+1+q); endr(K)=R(J);endif K-J==1&J-I~=1&K~=Nfor q=1:(J-I)r(I+1+q)=R(I+q); endr(I+1)=R(K);endif I==1&J==2&K==Nfor q=1:(N-2)r(1+q)=R(2+q);endr(N)=R(2);endif I==1&J==(N-1)&K==Nfor q=1:(N-2)r(q)=R(1+q);endr(N-1)=R(1);endif J-I~=1&K-I==N-1for q=1:(J-1)r(q)=R(1+q);endr(J)=R(1);endif J==(N-1)&K==N&J-I~=1r(J+1)=R(N);for q=1:(N-J-1)r(J+1+q)=R(J+q);endendif cost_sum(r,C,N)<cost_sum(R,C,N)R=rendendfprintf('总长为%f\n',cost_sum(R,C,N))%%%%%%cost_sum.m%%%%%%%%functiony=cost_sum(x,C,N)y=0;for i=1:(N-1)y=y+C(x(i),x(i+1));endy=y+C(x(N),x(1));三，灰色预测代码<pre name="code" class="plain">clearclcX=[136 143 165 152 165 181 204 272 319 491 571 605 665 640 628];x1(1)=X(1);X1=[];for i=1:1:14x1(i+1)=x1(i)+X(i+1);X1=[X1,x1(i)];endX1=[X1,X1(14)+X(15)]for k=3:1:15p(k)=X(k)/X1(k-1);p1(k)=X1(k)/X1(k-1);endp,p1clear kZ=[];for k=2:1:15z(k)=0.5*X1(k)+0.5*X1(k-1);Z=[Z,z(k)];endZB=[-Z',ones(14,1)]Y=[];clear ifor i=2:1:15Y=[Y;X(i)];endYA=inv(B'*B)*B'*Yclear ky1=[];for k=1:1:15y(k)=(X(1)-A(2)/A(1))*exp(-A(1)*(k-1))+A(2)/A(1); y1=[y1;y(k)];endy1clear kX2=[];for k=2:1:15x2(k)=y1(k)-y1(k-1);X2=[X2;x2(k)];endX2=[y1(1);X2]e=X'-X2m=abs(e)./X's=e'*en=sum(m)/13clear ksyms ky=(X(1)-A(2)/A(1))*exp(-A(1)*(k-1))+A(2)/A(1)Y1=[];for j=16:1:21y11=subs(y,k,j)-subs(y,k,j-1);Y1=[Y1;y11];endY1%程序中的变量定义:alpha是包含α、μ值的矩阵;%ago是预测后累加值矩阵;var是预测值矩阵;%error是残差矩阵; c是后验差比值function basicgrey(x,m) %定义函数basicgray(x)if nargin==1 %m为想预测数据的个数，默认为1 m=1;endclc; %清屏,以使计算结果独立显示if length(x(:,1))==1 %对输入矩阵进行判断,如不是一维列矩阵,进行转置变换x=x';endn=length(x); %取输入数据的样本量x1(:,1)=cumsum(x); %计算累加值,并将值赋及矩阵be for i=2:n %对原始数列平行移位 Y(i-1,:)=x(i,:);endfor i=2:n %计算数据矩阵B的第一列数据z(i,1)=0.5*x1(i-1,:)+0.5*x1(i,:);endB=ones(n-1,2); %构造数据矩阵BB(:,1)=-z(2:n,1);alpha=inv(B'*B)*B'*Y; %计算参数α、μ矩阵for i=1:n+m %计算数据估计值的累加数列,如改n+1为n+m可预测后m个值ago(i,:)=(x1(1,:)-alpha(2,:)/alpha(1,:))*exp(-alpha(1, :)*(i-1))+alpha(2,:)/alpha(1,:);endvar(1,:)=ago(1,:);f or i=1:n+m-1 %可预测后m个值var(i+1,:)=ago(i+1,:)-ago(i,:); %估计值的累加数列的还原,并计算出下m个预测值end[P,c,error]=lcheck(x,var); %进行后验差检验[rela]=relations([x';var(1:n)']); %关联度检验ago %显示输出预测值的累加数列alpha %显示输出参数α、μ数列var %显示输出预测值error %显示输出误差P %显示计算小残差概率 c %显示后验差的比值crela %显示关联度judge(P,c,rela) %评价函数显示这个模型是否合格<pre name="code" class="plain">function judge(P,c,rela) %评价指标并显示比较结果if rela>0.6'根据经验关联度检验结果为满意（关联度只是参考主要看后验差的结果）'else'根据经验关联度检验结果为不满意（关联度只是参考主要看后验差的结果）'endif P>0.95&c<0.5'后验差结果显示这个模型评价为“优”'else if P>0.8&c<0.5'后验差结果显示这个模型评价为“合格”'else if P>0.7&c<0.65'后验差结果显示这个模型评价为“勉强合格”' else'后验差结果显示这个模型评价为“不合格”' endendendfunction [P,c,error]=lcheck(x,var)%进行后验差检验n=length(x);for i=1:nerror(i,:)=abs(var(i,:)-x(i,:)); %计算绝对残差c=std(abs(error))/std(x); %调用统计工具箱的标准差函数计算后验差的比值cs0=0.6745*std(x);ek=abs(error-mean(error));pk=0;for i=1:nif ek(i,:)<s0pk=pk+1;endendP=pk/n; %计算小残差概率%附带的质料里有一部分讲了关联度function [rela]=relations(x)%以x(1,:)的参考序列求关联度[m,n]=size(x);for i=1:mfor j=n:-1:2x(i,j)=x(i,j)/x(i,1);endfor i=2:mx(i,:)=abs(x(i,:)-x(1,:)); %求序列差endc=x(2:m,:);Max=max(max(c)); %求两极差Min=min(min(c));p=0.5; %p称为分辨率,0<p<1,一般取p=0.5for i=1:m-1for j=1:nr(i,j)=(Min+p*Max)/(c(i,j)+p*Max); %计算关联系数endendfor i=1:m-1rela(i)=sum(r(i,:))/n; %求关联度end四，非线性拟合function f=example1(c,tdata)f=c(1)*(exp(-c(2)*tdata)-exp(-c(3)*tdata));<pre name="code" class="plain">function f=zhengtai(c,x) f=(1./(sqrt(2.*3.14).*c(1))).*exp(-(x-c(1)).^2./(2.*c( 2)^2));x=1:1:12;y=[01310128212]';c0=[2 8];for i=1:1000c=lsqcurvefit(@zhengtai,c0,x,y);c0=c;endy1=(1./(sqrt(2.*3.14).*c(1))).*exp(-(x-c(1)).^2./(2.*c (2)^2));plot(x,y,'r-',x,y1);legend('实验数据','拟合曲线')x=[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]';y=[30 68 75 82 82 77 68 68 58 51 50 41 38 35 28 25 18 15 12 10 7 7 4]';f=@(c,x)c(1)*(exp(-c(2)*x)-exp(-c(3)*x));c0=[114 0.1 2]';for i=1:50opt=optimset('TolFun',1e-3);[c R]=nlinfit(x,y,f,c0,opt)c0=c;hold onplot(x,c(1)*(exp(-c(2)*x)-exp(-c(3)*x)),'g')endt=[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];y=[30 68 75 82 82 77 68 68 58 51 50 41 38 35 28 25 18 15 12 10 7 7 4];c0=[1 1 1];for i=1:50 c=lsqcurvefit(@example1,c0,t,y);c0=c;endy1=c(1)*(exp(-c(2)*t)-exp(-c(3)*t));plot(t,y,' +',t,y1);legend('实验数据','拟合曲线')五，插值拟合相关知识在生产和科学实验中，自变量及因变量间的函数关系有时不能写出解析表达式，而只能得到函数在若干点的函数值或导数值，或者表达式过于复杂而需要较大的计算量。

数学建模算法的matlab代码句柄图形(图形窗口)二，hamiton回路算法提供一种求解最优哈密尔顿的算法---三边交换调整法，要求在运行jiaohuan3（三交换法）之前，给定邻接矩阵C和节点个数N，结果路径存放于R 中。

bianquan.m文件给出了一个参数实例，可在命令窗口中输入bianquan,得到邻接矩阵C和节点个数N以及一个任意给出的路径R,，回车后再输入jiaohuan3,得到了最优解。

由于没有经过大量的实验，又是近似算法，对于网络比较复杂的情况，可以尝试多运行几次jiaohuan3,看是否能到进一步的优化结果。

%%%%%%bianquan.m%%%%%%%N=13;for i=1:Nfor j=1:NC(i,j)=inf;endendfor i=1:NC(i,i)=0;endC(1,2)=6.0;C(1,13)=12.9;C(2,3)=5.9;C(2,4)=10.3;C(3,4)=12.2;C(3,5)=17.6;C(4,13)=8.8;C(4,7)=7.4;C(4,5)=11.5;C(5,2)=17.6;C(5,6)=8.2;C(6,9)=14.9;C(6,7)=20.3;C(7,9)=19.0;C(7,8)=7.3;C(8,9)=8.1;C(8,13)=9.2;C(9,10)=10.3;C(10,11)=7.7;C(11,12)=7.2;C(12,13)=7.9;for i=1:Nfor j=1:Nif C(i,j) < infC(j,i)=C(i,j);endendendfor i=1:NC(i,i)=0;endR=[4 7 6 5 3 2 1 13 12 11 10 9 8]; %%%%%%%%jiaohuan3.m%%%%%%%%%% n=0; for I=1:(N-2)for J=(I+1):(N-1)for K=(J+1):Nn=n+1;Z(n,:)=[I J K];endendendR=1:Nfor m=1:(N*(N-1)*(N-2)/6)I=Z(m,1);J=Z(m,2);K=Z(m,3);r=R;if J-I~=1&K-J~=1&K-I~=N-1for q=1:(J-I)r(I+q)=R(J+1-q);endfor q=1:(K-J)r(J+q)=R(K+1-q);endendif J-I==1&K-J==1r(K)=R(J);r(J)=R(K);endif J-I==1&K-J~=1&K-I~=N-1 for q=1:(K-J)r(I+q)=R(I+1+q);endr(K)=R(J);endif K-J==1&J-I~=1&K~=N for q=1:(J-I)r(I+1+q)=R(I+q);endr(I+1)=R(K);endif I==1&J==2&K==Nfor q=1:(N-2)r(1+q)=R(2+q);endr(N)=R(2);endif I==1&J==(N-1)&K==N for q=1:(N-2)r(q)=R(1+q);endr(N-1)=R(1);endif J-I~=1&K-I==N-1for q=1:(J-1)r(q)=R(1+q);endr(J)=R(1);endif J==(N-1)&K==N&J-I~=1r(J+1)=R(N);for q=1:(N-J-1)r(J+1+q)=R(J+q);endendif cost_sum(r,C,N)<cost_sum(r,c,n)< p="">R=rendendfprintf('总长为%f\n',cost_sum(R,C,N))%%%%%%cost_sum.m%%%%%%%%functiony=cost_sum(x,C,N)y=0;for i=1:(N-1) y=y+C(x(i),x(i+1));endy=y+C(x(N),x(1));三，灰色预测代码clearclcX=[136 143 165 152 165 181 204 272 319 491 571 605 665 640 628];x1(1)=X(1);X1=[];for i=1:1:14x1(i+1)=x1(i)+X(i+1);X1=[X1,x1(i)];endX1=[X1,X1(14)+X(15)]for k=3:1:15p(k)=X(k)/X1(k-1);p1(k)=X1(k)/X1(k-1);endp,p1clear kZ=[];for k=2:1:15z(k)=0.5*X1(k)+0.5*X1(k-1);Z=[Z,z(k)];endZB=[-Z',ones(14,1)]Y=[];clear ifor i=2:1:15Y=[Y;X(i)];endYA=inv(B'*B)*B'*Yclear ky1=[];for k=1:1:15y(k)=(X(1)-A(2)/A(1))*exp(-A(1)*(k-1))+A(2)/A(1); y1=[y1;y(k)];endy1clear kX2=[];for k=2:1:15x2(k)=y1(k)-y1(k-1);X2=[X2;x2(k)];endX2=[y1(1);X2]e=X'-X2m=abs(e)./X's=e'*en=sum(m)/13clear ksyms ky=(X(1)-A(2)/A(1))*exp(-A(1)*(k-1))+A(2)/A(1)Y1=[];for j=16:1:21y11=subs(y,k,j)-subs(y,k,j-1);Y1=[Y1;y11];endY1%程序中的变量定义:alpha是包含α、μ值的矩阵;%ago是预测后累加值矩阵;var是预测值矩阵;%error是残差矩阵; c是后验差比值function basicgrey(x,m) %定义函数basicgray(x)if nargin==1 %m 为想预测数据的个数，默认为1 m=1;endclc; %清屏,以使计算结果独立显示if length(x(:,1))==1 %对输入矩阵进行判断,如不是一维列矩阵,进行转置变换x=x';endn=length(x); %取输入数据的样本量x1(:,1)=cumsum(x); %计算累加值,并将值赋与矩阵be for i=2:n %对原始数列平行移位Y(i-1,:)=x(i,:);endfor i=2:n %计算数据矩阵B的第一列数据z(i,1)=0.5*x1(i-1,:)+0.5*x1(i,:);endB=ones(n-1,2); %构造数据矩阵BB(:,1)=-z(2:n,1);alpha=inv(B'*B)*B'*Y; %计算参数α、μ矩阵for i=1:n+m %计算数据估计值的累加数列,如改n+1为n+m可预测后m个值ago(i,:)=(x1(1,:)-alpha(2,:)/alpha(1,:))*exp(-alpha(1,:)*(i-1))+alpha(2,:)/alpha(1,:);endvar(1,:)=ago(1,:);for i=1:n+m-1 %可预测后m个值var(i+1,:)=ago(i+1,:)-ago(i,:); %估计值的累加数列的还原,并计算出下m个预测值end[P,c,error]=lcheck(x,var); %进行后验差检验[rela]=relations([x';var(1:n)']); %关联度检验ago %显示输出预测值的累加数列alpha %显示输出参数α、μ数列var %显示输出预测值error %显示输出误差P %显示计算小残差概率 c %显示后验差的比值crela %显示关联度judge(P,c,rela) %评价函数显示这个模型是否合格function judge(P,c,rela)%评价指标并显示比较结果if rela>0.6'根据经验关联度检验结果为满意（关联度只是参考主要看后验差的结果）' else'根据经验关联度检验结果为不满意（关联度只是参考主要看后验差的结果）' endif P>0.95&c<0.5'后验差结果显示这个模型评价为“优”'else if P>0.8&c<0.5'后验差结果显示这个模型评价为“合格”'else if P>0.7&c<0.65'后验差结果显示这个模型评价为“勉强合格”'else'后验差结果显示这个模型评价为“不合格”'endendendfunction [P,c,error]=lcheck(x,var)%进行后验差检验n=length(x);for i=1:nerror(i,:)=abs(var(i,:)-x(i,:)); %计算绝对残差endc=std(abs(error))/std(x); %调用统计工具箱的标准差函数计算后验差的比值cs0=0.6745*std(x);ek=abs(error-mean(error));pk=0;for i=1:nif ek(i,:)<s0< p="">pk=pk+1;endendP=pk/n; %计算小残差概率%附带的质料里有一部分讲了关联度function [rela]=relations(x)%以x(1,:)的参考序列求关联度[m,n]=size(x);for i=1:mfor j=n:-1:2x(i,j)=x(i,j)/x(i,1);endendfor i=2:mx(i,:)=abs(x(i,:)-x(1,:)); %求序列差endc=x(2:m,:);Max=max(max(c)); %求两极差Min=min(min(c));p=0.5; %p称为分辨率,0<p< p="">for i=1:m-1for j=1:nr(i,j)=(Min+p*Max)/(c(i,j)+p*Max); %计算关联系数endendfor i=1:m-1rela(i)=sum(r(i,:))/n; %求关联度end四，非线性拟合function f=example1(c,tdata)f=c(1)*(exp(-c(2)*tdata)-exp(-c(3)*tdata));function f=zhengtai(c,x)f=(1./(sqrt(2.*3.14).*c(1))).*exp(-(x-c(1)).^2./(2.*c(2)^2));x=1:1:12;y=[01310128212]';c0=[2 8];for i=1:1000c=lsqcurvefit(@zhengtai,c0,x,y);c0=c;endy1=(1./(sqrt(2.*3.14).*c(1))).*exp(-(x-c(1)).^2./(2.*c(2)^2));plot(x,y,'r-',x,y1);legend('实验数据','拟合曲线')x=[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]';y=[30 68 75 82 82 77 68 68 58 51 50 41 38 35 28 25 18 15 12 10 7 7 4]';f=@(c,x)c(1)*(exp(-c(2)*x)-exp(-c(3)*x));c0=[114 0.1 2]';for i=1:50opt=optimset('TolFun',1e-3);[c R]=nlinfit(x,y,f,c0,opt)c0=c;hold onplot(x,c(1)*(exp(-c(2)*x)-exp(-c(3)*x)),'g')endt=[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 1314 15 16];y=[30 68 75 82 82 77 68 68 58 51 50 41 38 35 28 25 1815 12 10 7 7 4];c0=[1 1 1];for i=1:50 c=lsqcurvefit(@example1,c0,t,y);c0=c;endy1=c(1)*(exp(-c(2)*t)-exp(-c(3)*t));plot(t,y,'+',t,y1);legend('实验数据','拟合曲线')五，插值拟合相关知识在生产和科学实验中，自变量与因变量间的函数关系有时不能写出解析表达式，而只能得到函数在若干点的函数值或导数值，或者表达式过于复杂而需要较大的计算量。