1、数学实验四 学 生 实 验 报 告实验课程名称 开课实验室 学 院 年级 专业班 学 生 姓 名 刘蛰 学 号 开 课 时 间 至 学年第 学期总 成 绩教师签名数 学 与 统 计 学 院 制开课学院、实验室： 实验时间 ： 年 月 日课程名称实验项目名 称实验项目类型验证演示综合设计其他指导教师成 绩实验目的1 了解最小二乘拟合的基本原理和方法；2 掌握用MATLAB作最小二乘多项式拟合和曲线拟合的方法；3 通过实例学习如何用拟合方法解决实际问题，注意与插值方法的区别。4 了解各种参数辨识的原理和方法；5 通过范例展现由机理分析确定模型结构，拟合方法辨识参数，误差分析等求解实际问题的过程；基
2、础实验一、实验内容1. 下表给出了近两个世纪某国人口的统计数据，试预测2010年该国的人口。1) 可先用以上数据拟合Multhus人口指数增长模型，根据检验结果进一步讨论马尔萨斯人口模型的改进。2) Malthus 模型的基本假设是：人口的增长率为常数，记为r。记时刻t的人口为x(t），且初始时刻的人口为x0，于是得到如下微分方程：3) Malthus 模型在短期内能比较好的拟合数据，但从长期来看，该模型的增长速度过快导致和实际的常识不相符。为此，一个改进是考虑种群内部的竞争。比如2.某年美国旧车价格的调查资料如下表其中xi表示轿车的使用年数，yi表示相应的平均价格。试分析用什么形式的曲线来拟
3、合上述的数据，并预测使用4.5年后轿车的平均价格大致为多少？二、实验过程（一般应包括实验原理或问题分析，算法设计、程序、计算、图表等， 实验结果及分析）由已知的微分方程模型可以求得人口数与时间的函数关系，运用matlab的polyfit、lsqcurvefit函数，即可利用已知数据求得未知参数，使得该曲线与所给数据误差最小。在求得的人口与时间函数关系中，对数据稍作处理，可以将问题转化为一次线性拟合。第一题第一小题M文件：t=linspace(0,20,21);x=3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76 92 106.5 123
4、.2 131.7 150.7 179.3 204 226.5 251.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)grid onyuce=exp(yuce)运行结果截图：运行结果显示：预测人口数为564.7065，很明显，由于模型过于粗糙，没有考虑种群竞争及环境容纳量，预测值明显偏大。第一题第二小题M文件：xdata=linspace(0,20,21);t=0:0.1:23;ydata=3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.
5、6 50.2 62.9 76 92 106.5 123.2 131.7 150.7 179.3 204 226.5 251.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)grid on名字为funfun的M文件：function f=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)*(1
6、0*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)grid on名字为jia的M文件：function f=j
7、ia(x,xdata)f=x(1).*exp(x(2).*xdata);运行结果截图：预测值为46.3305.曲线与数据吻合的较好。应用实验（或综合实验）一、实验内容平面上的圆可以由三个参数确定，半径和圆心和坐标。通过提取圆上点的坐标也能确定圆的方程，而最少的点数是3个点。为减小误差，在实际中一般会得到多个点的坐标通过拟合的方法得到圆的方程。下面的数据来源于一个圆，使用这些数据确定圆的方程。x=(14.5081 13.7508 13.3520 11.6870 10.3406)y=(1.8410 3.6362 3.9491 4.3287 3.7139).二、问题分析 对于这个问题，我分析了很久。
8、最初直接把y从圆的一般方程里将y解出来。，用此等式直接作为函数，供lsqcurvefit调用。但是运行结果中出现复数。之后我想用参数方程 表示圆，可是t的数据未知，已知的两组数据时x，y的。再后来想到用极坐标，把带入一般方程，将数据x，y处理成但是形式与第一种差不多，个人认为应该还是同样的结果，便没有编程序运行。最终只得在网上寻找一段圆拟合的代码，尝试着学习一下。三、数学模型的建立与求解(一般应包括模型、求解步骤或思路，程序放在后面的附录中) 5号宋体四、实验结果及分析 运行结果截图：五、附录（程序等）M文件：t=0:0.01:2*pi;xdata=14.5081 13.7508 13.352
9、0 11.6870 10.3406;ydata=1.8410 3.6362 3.9491 4.3287 3.7139;A=xdata ydata;Par=CircleFitByTaubin(A)x=Par(1)+Par(3)*cos(t);y=Par(2)+Par(3)*sin(t);plot(xdata,ydata,r+,x,y)axis equalaxis (9 15 -1 5)grid on名字为CircleFitByTaubin的M文件：function Par = CircleFitByTaubin(XY) %-% % Circle fit by Taubin% G. Taubin,
10、 Estimation Of Planar Curves, Surfaces And Nonplanar% Space Curves Defined By Implicit Equations, With % Applications To Edge And Range Image Segmentation,% IEEE Trans. PAMI, Vol. 13, pages 1115-1138, (1991)% Input: XY(n,2) is the array of coordinates of n points x(i)=XY(i,1), y(i)=XY(i,2)% Output:
11、Par = a b R is the fitting circle:% center (a,b) and radius R% Note: this fit does not use built-in matrix functions (except mean),% so it can be easily programmed in any programming language%- n = size(XY,1); % number of data points centroid = mean(XY); % the centroid of the data set % computing mo
12、ments (note: all moments will be normed, i.e. divided by n) Mxx = 0; Myy = 0; Mxy = 0; Mxz = 0; Myz = 0; Mzz = 0; for i=1:n Xi = XY(i,1) - centroid(1); % centering data Yi = XY(i,2) - centroid(2); % centering data Zi = Xi*Xi + Yi*Yi; Mxy = Mxy + Xi*Yi; Mxx = Mxx + Xi*Xi; Myy = Myy + Yi*Yi; Mxz = Mxz
13、 + 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; % computing the coefficients of the characteristic polynomial 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*
14、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; % Newtons method starting at x=0 for iter=1:IterMax yold = ynew; ynew = A0 + xnew*(A1 + xnew*(A2 + xnew*A3); if abs(ynew) ab
15、s(yold) disp(Newton-Taubin goes wrong direction: |ynew| |yold|); xnew = 0; break; end Dy = A1 + xnew*(A22 + xnew*A33); xold = xnew; xnew = xold - ynew/Dy; if (abs(xnew-xold)/xnew) = IterMax) disp(Newton-Taubin will not converge); xnew = 0; end if (xnew0.) fprintf(1,Newton-Taubin negative root: x=%fn
16、,xnew); xnew = 0; endend % computing the circle parameters 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人合作完成，请在实验报告上注明合作者的姓名。
