数学实验四.docx
《数学实验四.docx》由会员分享，可在线阅读，更多相关《数学实验四.docx（13页珍藏版）》请在冰点文库上搜索。
数学实验四
学生实验报告
实验课程名称
开课实验室
学院年级专业班
学生姓名刘蛰学号
开课时间至学年第学期
总成绩
教师签名
数学与统计学院制
开课学院、实验室：
实验时间：
年月日
课程
名称
实验项目
名称
实验项目类型
验证
演示
综合
设计
其他
指导
教师
成绩
实验目的
[1]了解最小二乘拟合的基本原理和方法；
[2]掌握用MATLAB作最小二乘多项式拟合和曲线拟合的方法；
[3]通过实例学习如何用拟合方法解决实际问题，注意与插值方法的区别。
[4]了解各种参数辨识的原理和方法；
[5]通过范例展现由机理分析确定模型结构，拟合方法辨识参数，误差分析等求解实际问题的过程；
基础实验
一、实验内容
1.下表给出了近两个世纪某国人口的统计数据，试预测2010年该国的人口。
1）可先用以上数据拟合Multhus人口指数增长模型，根据检验结果进一步讨论马尔萨斯人口模型的改进。
2）Malthus模型的基本假设是：
人口的增长率为常数，记为r。
记时刻t的人口为x（t），且初始时刻的人口为x0，于是得到如下微分方程：
3）Malthus模型在短期内能比较好的拟合数据，但从长期来看，该模型的增长速度过快导致和实际的常识不相符。
为此，一个改进是考虑种群内部的竞争。
比如
2.某年美国旧车价格的调查资料如下表其中xi表示轿车的使用年数，yi表示相应的平均价格。
试分析用什么形式的曲线来拟合上述的数据，并预测使用4.5年后轿车的平均价格大致为多少？
二、实验过程（一般应包括实验原理或问题分析，算法设计、程序、计算、图表等，实验结果及分析）
由已知的微分方程模型可以求得人口数与时间的函数关系，运用matlab的polyfit、lsqcurvefit函数，即可利用已知数据求得未知参数，使得该曲线与所给数据误差最小。
在求得的人口与时间函数关系中，对数据稍作处理，可以将问题转化为一次线性拟合。
第一题第一小题M文件：
t=linspace（0,20,21）;
x=[3.95.37.29.612.917.123.231.438.650.262.97692106.5123.2131.7150.7179.3204226.5251.4];
x=log（x）;
a=polyfit（t,x,1）
y=polyval（a,t）;
yuce=polyval（a,22）;
plot（t,x,'k+',t,y,22,yuce,'ro'）
gridon
yuce=exp（yuce）
运行结果截图：
运行结果显示：
预测人口数为564.7065，很明显，由于模型过于粗糙，没有考虑种群竞争及环境容纳量，预测值明显偏大。
第一题第二小题M文件：
xdata=linspace（0,20,21）;
t=0:
0.1:
23;
ydata=[3.95.37.29.612.917.123.231.438.650.262.97692106.5123.2131.7150.7179.3204226.5251.4];
a0=[.5,1];
[a]=lsqcurvefit（'funfun',a0,xdata,ydata）
yuce=funfun（a,22）
plot（xdata,ydata,'r+',t,funfun（a,t）,22,yuce,'ro'）
gridon
名字为funfun的M文件：
functionf=funfun（x,xdata）;
%f=dsolve（'Dy=x
（1）*y*（1-y/x
（2））','y（0）=3.9','xdata'）
f=x
（2）./（1+1/39.*exp（-x
（1）.*xdata）*（10*x
（2）-39））;
运行结果截图：
M文件中的注释部分为求解微分方程的代码。
微分方程求解过后，将运算符改为.*,./，即可写入函数。
预测结果为267.1928，与实际值比较接近。
第二题M文件：
a0=[12-10];
t=0:
.1:
15;
xdata=1:
10;
ydata=[2615,1943,1494,1087,765,538,484,290,226,204];
[a]=lsqcurvefit（'jia',a0,xdata,ydata）
yuce=jia（a,14.5）
plot（xdata,ydata,'r+',t,jia（a,t）,14.5,yuce,'ro'）
gridon
名字为jia的M文件：
functionf=jia（x,xdata）
f=x
（1）.*exp（x
（2）.*xdata）;
运行结果截图：
预测值为46.3305.曲线与数据吻合的较好。
应用实验（或综合实验）
一、实验内容
平面上的圆可以由三个参数确定，半径和圆心和坐标。
通过提取圆上点的坐标也能确定圆的方程，而最少的点数是3个点。
为减小误差，在实际中一般会得到多个点的坐标通过拟合的方法得到圆的方程。
下面的数据来源于一个圆，使用这些数据确定圆的方程。
x=（14.508113.750813.352011.687010.3406）
y=（1.84103.63623.94914.32873.7139）.
二、问题分析
对于这个问题，我分析了很久。
最初直接把y从圆的一般方程
里将y解出来。
，用此等式直接作为函数，供lsqcurvefit调用。
但是运行结果中出现复数。
之后我想用参数方程
表示圆，可是t的数据未知，已知的两组数据时x，y的。
再后来想到用极坐标，把
带入一般方程，将数据x，y处理成
但是形式与第一种差不多，个人认为应该还是同样的结果，便没有编程序运行。
最终只得在网上寻找一段圆拟合的代码，尝试着学习一下。
三、数学模型的建立与求解（一般应包括模型、求解步骤或思路，程序放在后面的附录中）
5号宋体
四、实验结果及分析
运行结果截图：
五、附录（程序等）
M文件：
t=0:
0.01:
2*pi;
xdata=[14.508113.750813.352011.687010.3406]';
ydata=[1.84103.63623.94914.32873.7139]';
A=[xdataydata];
Par=CircleFitByTaubin（A）
x=Par
（1）+Par（3）*cos（t）;
y=Par
（2）+Par（3）*sin（t）;
plot（xdata,ydata,'r+',x,y）
axisequal
axis（[915-15]）
gridon
名字为CircleFitByTaubin的M文件：
functionPar=CircleFitByTaubin（XY）
%--------------------------------------------------------------------------
%
%CirclefitbyTaubin
%G.Taubin,"EstimationOfPlanarCurves,SurfacesAndNonplanar
%SpaceCurvesDefinedByImplicitEquations,With
%ApplicationsToEdgeAndRangeImageSegmentation",
%IEEETrans.PAMI,Vol.13,pages1115-1138,（1991）
%
%Input:
XY（n,2）isthearrayofcoordinatesofnpointsx（i）=XY（i,1）,y（i）=XY（i,2）
%
%Output:
Par=[abR]isthefittingcircle:
%center（a,b）andradiusR
%
%Note:
thisfitdoesnotusebuilt-inmatrixfunctions（except"mean"）,
%soitcanbeeasilyprogrammedinanyprogramminglanguage
%
%--------------------------------------------------------------------------
n=size（XY,1）;%numberofdatapoints
centroid=mean（XY）;%thecentroidofthedataset
%computingmoments（note:
allmomentswillbenormed,i.e.dividedbyn）
Mxx=0;Myy=0;Mxy=0;Mxz=0;Myz=0;Mzz=0;
fori=1:
n
Xi=XY（i,1）-centroid
（1）;%centeringdata
Yi=XY（i,2）-centroid
（2）;%centeringdata
Zi=Xi*Xi+Yi*Yi;
Mxy=Mxy+Xi*Yi;
Mxx=Mxx+Xi*Xi;
Myy=Myy+Yi*Yi;
Mxz=Mxz+Xi*Zi;
Myz=Myz+Yi*Zi;
Mzz=Mzz+Zi*Zi;
end
Mxx=Mxx/n;
Myy=Myy/n;
Mxy=Mxy/n;
Mxz=Mxz/n;
Myz=Myz/n;
Mzz=Mzz/n;
%computingthecoefficientsofthecharacteristicpolynomial
Mz=Mxx+Myy;
Cov_xy=Mxx*Myy-Mxy*Mxy;
A3=4*Mz;
A2=-3*Mz*Mz-Mzz;
A1=Mzz*Mz+4*Cov_xy*Mz-Mxz*Mxz-Myz*Myz-Mz*Mz*Mz;
A0=Mxz*Mxz*Myy+Myz*Myz*Mxx-Mzz*Cov_xy-2*Mxz*Myz*Mxy+Mz*Mz*Cov_xy;
A22=A2+A2;
A33=A3+A3+A3;
xnew=0;
ynew=1e+20;
epsilon=1e-12;
IterMax=20;
%Newton'smethodstartingatx=0
foriter=1:
IterMax
yold=ynew;
ynew=A0+xnew*（A1+xnew*（A2+xnew*A3））;
ifabs（ynew）>abs（yold）
disp（'Newton-Taubingoeswrongdirection:
|ynew|>|yold|'）;
xnew=0;
break;
end
Dy=A1+xnew*（A22+xnew*A33）;
xold=xnew;
xnew=xold-ynew/Dy;
if（abs（（xnew-xold）/xnew）if（iter>=IterMax）
disp（'Newton-Taubinwillnotconverge'）;
xnew=0;
end
if（xnew<0.）
fprintf（1,'Newton-Taubinnegativeroot:
x=%f\n',xnew）;
xnew=0;
end
end
%computingthecircleparameters
DET=xnew*xnew-xnew*Mz+Cov_xy;
Center=[Mxz*（Myy-xnew）-Myz*Mxy,Myz*（Mxx-xnew）-Mxz*Mxy]/DET/2;
Par=[Center+centroid,sqrt（Center*Center'+Mz）];
end%CircleFitByTaubin
总结与体会
教师签名
年月日
备注：
1、同一章的实验作为一个实验项目，每个实验做完后提交电子稿到服务器的“全校任选课数学实验作业提交”文件夹，文件名为“学院学号姓名实验几”，如“机械20073159张新实验一”。
2、提交的纸质稿要求双面打印，中途提交批改不需要封面，但最后一次需将该课程所有实验项目内页与封面一起装订成册提交。
3、综合实验要求3人合作完成，请在实验报告上注明合作者的姓名。
升级会员 