数学实验四.docx
- 文档编号:14798829
- 上传时间:2023-06-27
- 格式:DOCX
- 页数:13
- 大小:108.97KB
数学实验四.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人合作完成,请在实验报告上注明合作者的姓名。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数学 实验