下载文档原格式
内蒙古科技大学本科生毕业设计说明书（毕业论文）题目：最小二乘曲线拟合及MATLAB实现学生姓名：李亚学号：0972143230专业：测绘工程班级：2009测绘2班指导教师：孙同贺讲师最小二乘曲线拟合及MATLAB实现摘要介绍曲线拟合的基本理论，对最小二乘原理进行了全方位的理论阐述，同时也阐述了曲线拟合的基本原理及多项式曲线拟合模型的建立。
详细的解答了曲线拟合中的最小二乘法，并介绍了部分的正交最小二乘法理论。
重点讲解多项式拟合的具体步骤，同时也介绍了非线性方程的最小二乘拟合，在建立理论的基础上对最小二乘曲线拟合法的MATLAB实现方法进行研究，利用MATLAB2012b的平台对测量数据进行最小二乘曲线拟合，介绍MATLAB的具体构造和曲线拟合工具。
利用MATLAB中的ployfit函数对实测数据进行多项式曲线拟合，并给出曲线拟合MATLAB实现的源程序，给出拟合曲线，并评定拟合的精度证明该方法是行之有效的。
关键词：最小二乘法，曲线拟合，MATLAB，测量数据Curve Fitting in Least-Square Methodand Its Realization with MatlabAbstractTo introduce the basic theory of curve fitting and discuss the least squares principle in this paper, what’s more, we also discuss the basic principle of curve fitting and the establishment of polynomial curve fitting model. Meanwhile, we also introduce the least-square method of curve fitting in detail and part of the theory of orthogonal least square method. We mainly discuss the specific steps of polynomial fitting, and also introduces the nonlinear equation of the least squares fitting at the same time, which established on the theory of least squares curve fitting in MATLAB in order to realize the method to do research. Using MATLAB2012b platform to achieve the goal of measuring data and introducing the special structure of MATLAB and curve fitting tool. We can use ployfit function in MATLAB to polynomial curve fitting of experimental data, and get the MATLAB source program about curve fitting and the fitting curve. Finally, we need to prove the method of assessing the precision of the fitting is effective.Key words: least square method; curve fitting; MATLAB, metrical data最小二乘曲线拟合及MATLAB实现 (I)摘要 (I)CURVE FITTING IN LEAST-SQUARE METHOD AND ITS REALIZATION WITH MATLAB (II)ABSTRACT (II)第一章引言 (1)1.1研究背景 (1)1.1.1 历史理论原理 (1)1.1.2 现代研究 (1)1.2问题定义 (2)1.2.1 曲线拟合的思想 (2)1.2.2 多项式拟合 (3)1.2.3 利用Matlab的polyfit函数进行多项式拟合 (3)1.3论文结构 (3)第二章数据曲线拟合 (4)2.1测量数据 (4)2.2拟合模型 (4)2.3最小二乘原理 (5)2.3.1最小二乘法 (5)2.3.2最小二乘估计与极大似然估计 (7)2.4.1曲线拟合理论 (9)2.4.2最小二乘法线性拟合原理 (10)2.4.3最小二乘非线性拟合 (12)2.4.4正交多项式 (13)2.4.5正交最小二乘曲线拟合 (15)2.5曲线拟合精度评定 (17)第三章MATLAB (19)3.1MATLAB概述 (19)3.1.1MATLAB简介 (19)3.1.2MATLAB的主要组成部分 (21)3.2MATLAB2012B的运行简介 (23)3.2.1启动和退出MATLAB2012b (23)3.2.2MATLAB2012b桌面系统 (24)3.2.3MATLAB函数调用系统 (26)3.2.4MATLAB2012b的帮助系统 (27)3.2.5附件管理系统 (28)3.2.6数据交换系统 (28)3.2.7MATLAB 中的其他系统 (29)3.3最小二乘曲线拟合法的MATLAB实现 (30)第四章最小二乘法曲线拟合的MATLAB实现 (32)4.1使用POLYFIT函数实现多项式拟合 (32)4.3三次多项式的曲线拟合 (34)4.4四次多项式曲线拟合 (35)4.5数据处理和精度评定 (36)第五章总结 (40)参考文献 (41)附录1： (43)MATLAB语言编程源代码 (43)附录2： (45)各次拟合的拟合曲线方程 (45)致谢 (46)外文翻译 (47)外文部分 (47)翻译部分 (54)第一章 引 言1.1研究背景1.1.1 历史理论原理Weierstrass 第一逼近定理[1]对任意函数()f C[a,b]x ∈和任意给定的0ε>,都存在n 次代数多项式()n p x ，满足()()[,]||<max nx a b f x p x ε∈- (1-1-1) Bernstein 多项式（bernstein polynomial ）[1]前苏联数学家Bernstein 曾经给出这样的多项式序列：0(,)(1)n k n k n k n k B f x f x x k n -=⎛⎫⎛⎫=- ⎪ ⎪⎝⎭⎝⎭∑ (1-1-2)在整体上一致逼近()f x ，但它的收敛缓慢，要达到一定的精度，则n 要取很大，计算量大，所以研究如何在给定的精度下，对()f x 进行整体逼近，成为逼近论中的一个重要问题。
曲线拟合（curve-fitting ）：工程实践中，用测量到的一些离散的数据},...2,1,0),,{(m i y x i i =求一个近似的函数)(x ϕ来拟合这组数据，要求所得的拟合曲线能最好的反映数据的基本趋势（即使)(x ϕ最好地逼近()x f ，而不必满足插值原则。
因此没必要取)(i x ϕ=i y ，只要使i i i y x -=)(ϕδ尽可能地小）。
原理：给定数据点},...2,1,0),,{(m i y x i i =。
求近似曲线)(x ϕ。
并且使得近似曲线与()x f 的偏差最小。
近似曲线在该点处的偏差i i i y x -=)(ϕδ，i=1,2,...,m 。
常见的曲线拟合方法:1.使偏差绝对值之和最小2.使偏差绝对值最大的最小3.使偏差平方和最小最小二乘法：按偏差平方和最小的原则选取拟合曲线，并且采取二项式方程为拟合曲线的方法,称为最小二乘法。
推导过程：1. 设拟合多项式为：2. 各点到这条曲线的距离之和，即偏差平方和如下：3. 问题转化为求待定系数0a ...k a 对等式右边求i a 偏导数，因而我们得到了： .......4、 把这些等式化简并表示成矩阵的形式，就可以得到下面的矩阵：5. 将这个范德蒙得矩阵化简后可得到:6. 也就是说X*A=Y ，那么A = (X'*X)-1*X'*Y ，便得到了系数矩阵A ，同时，我们也就得到了拟合曲线。
MATLAB 实现：MATLAB 提供了polyfit （）函数命令进行最小二乘曲线拟合。
调用格式：p=polyfit(x,y,n)[p,s]= polyfit(x,y,n)[p,s,mu]=polyfit(x,y,n)x,y 为数据点，n 为多项式阶数，返回p 为幂次从高到低的多项式系数向量p 。
x 必须是单调的。
矩阵s 包括R （对x 进行QR 分解的三角元素）、df(自由度）、normr(残差）用于生成预测值的误差估计。
matlab最小二乘法曲线拟合最小二乘法是一种常用的曲线拟合方法，它通过最小化实际观测值与拟合曲线之间的平方误差来确定最佳拟合曲线的参数。
给定一组实际观测数据点(xi, yi)，我们希望找到一个拟合曲线y=f(x;θ)，其中θ表示曲线的参数。
最小二乘法的目标是使误差的平方和最小化，即使得下述损失函数最小化：L(θ) = ∑(yi - f(xi;θ))^2其中，∑表示求和运算，L(θ)是损失函数，yi是第i个观测数据点的输出值，f(xi;θ)是根据参数θ计算得到的拟合曲线在第i个观测点的预测值。
为了找到最佳的参数θ，我们通过最小化损失函数来求解优化问题：minimize L(θ)这个问题可以通过求解等式∂L/∂θ = 0 来得到最优参数θ的闭式解。
具体的求解方法，可以通过矩阵和向量的运算来实现。
在Matlab中，可以使用“polyfit”函数进行最小二乘法的曲线拟合。
该函数可以拟合出一条多项式曲线，通过指定最佳拟合的次数，即多项式的阶数。
拟合结果包括最佳参数和拟合误差等信息。
使用方法如下：```% 输入观测数据x = [x1, x2, x3, ...]';y = [y1, y2, y3, ...]';% 拟合曲线order = 1; % 最佳拟合的次数（如一次线性拟合）p = polyfit(x, y, order);% 最佳参数coefficients = p;% 拟合曲线curve = polyval(p, x);% 绘制拟合曲线和观测数据plot(x, y, 'o', x, curve)```这样，就可以使用Matlab的最小二乘法曲线拟合方法来得到最佳的拟合曲线。
曲线拟合（curve-fitting ）：工程实践中，用测量到的一些离散的数据},...2,1,0),,{(m i y x i i =求一个近似的函数)(x ϕ来拟合这组数据，要求所得的拟合曲线能最好的反映数据的基本趋势（即使)(x ϕ最好地逼近()x f ，而不必满足插值原则。
因此没必要取)(i x ϕ=i y ，只要使i i i y x -=)(ϕδ尽可能地小）。
原理：给定数据点},...2,1,0),,{(m i y x i i =。
求近似曲线)(x ϕ。
并且使得近似曲线与()x f 的偏差最小。
近似曲线在该点处的偏差i i i y x -=)(ϕδ，i=1,2,...,m 。
常见的曲线拟合方法:1.使偏差绝对值之和最小2.使偏差绝对值最大的最小3.使偏差平方和最小最小二乘法：按偏差平方和最小的原则选取拟合曲线，并且采取二项式方程为拟合曲线的方法,称为最小二乘法。
推导过程：1. 设拟合多项式为：kk x a x a a x +++=...)(10ϕ2. 各点到这条曲线的距离之和，即偏差平方和如下：3. 问题转化为求待定系数0a ...k a 对等式右边求i a 偏导数，因而我们得到了：.......4、 把这些等式化简并表示成矩阵的形式，就可以得到下面的矩阵：5. 将这个范德蒙得矩阵化简后可得到:6. 也就是说X*A=Y ，那么A = (X'*X)-1*X'*Y ，便得到了系数矩阵A ，同时，我们也就得到了拟合曲线。
MATLAB实现：MATLAB提供了polyfit（）函数命令进行最小二乘曲线拟合。
调用格式：p=polyfit(x,y,n)[p,s]= polyfit(x,y,n)[p,s,mu]=polyfit(x,y,n)x,y为数据点，n为多项式阶数，返回p为幂次从高到低的多项式系数向量p。
x 必须是单调的。
矩阵s包括R（对x进行QR分解的三角元素）、df(自由度）、normr(残差）用于生成预测值的误差估计。
matlab最小二乘法拟合曲线Matlab最小二乘法拟合曲线是一种应用于数据拟合的有效的工具，它的作用是使用最小二乘法来估计未知参数并获得适合拟合的最优拟合曲线，以下是Matlab最小二乘法拟合曲线的具体用法：一、Matlab最小二乘法拟合模型：1、首先，根据需要拟合的数据，定义未知参数的类型、数量和频率；2、接下来，定义未知参数的初始值，以及用于确定参数最优拟合曲线的搜索算法；3、然后，调用最小二乘法函数，使用最小二乘法函数计算拟合参数θ；4、最后，用优化到的θ值生成最优曲线，即得到拟合曲线。
二、Matlab最小二乘法拟合曲线的特点：1、精度高：最小二乘法在误差估计上是最佳的，能控制估计偏差，通过求解思维运算完成最小二乘拟合；2、可以处理多元数据：最小二乘法可以处理多个变量进行统计拟合，有多个自变量时，仍然能生成反映变量之间关系的拟合曲线；3、计算量小：最小二乘法只需计算发生一次，消耗计算量较小，计算正确率高；4、反应速度快：最小二乘法反应速度快，可以很好的拟合多项式，某一特定点的拟合能力强，它具有很高的拟合度。
三、Matlab最小二乘法拟合曲线的应用：1、最小二乘法拟合曲线可以用于多元统计拟合，研究变量之间的关系，可用于实验数据处理和建模；2、最小二乘法拟合曲线也可以用于经济学，可以通过估计最小二乘回归系数进行广义线性模型的预测；3、最小二乘法拟合曲线可以用于工程曲线拟合，如机械设计的几何拟合等，以及测量仪器的校正等；4、最小二乘法拟合曲线也可以用于生物学研究，可以通过进化分类树及类群的状态估计其特征变化趋势；5、最小二乘法拟合曲线还可以用于物理和化学实验中，以及天气、气候等领域。
四、Matlab最小二乘法拟合曲线的优缺点：优点：1、计算量小，计算消耗较小；2、可对多元数据进行拟合，处理变量之间的关系；3、拟合精度高，控制估计偏差；4、反应速度快，容错性强。
缺点：1、处理误差较大的数据时，拟合效果不佳；2、对曲线的凸性要求，不能处理异常数据；3、无法处理变量间的非线性关系，拟合结果也会出现偏差。
一、概述Matlab是一种强大的数学软件，它提供了许多用于数学建模、数据分析、图形可视化等功能的工具。
对数最小二乘拟合是其中一个重要的功能，它可以帮助研究人员处理实验数据，找出数据之间的相关性，从而进行有效的数据分析和建模。
在本文中，我们将讨论Matlab中对数最小二乘拟合的原理、方法和应用。
二、对数最小二乘拟合的原理对数最小二乘拟合是一种特殊的非线性拟合方法，它适用于当实验数据呈现出指数增长或指数衰减的趋势时。
对数最小二乘拟合的原理是将实验数据取对数变换，然后使用最小二乘法进行拟合。
最小二乘法是一种常用的数值优化方法，它通过最小化实际观测值和模型预测值之间的残差平方和来确定模型参数，从而实现拟合。
三、Matlab中对数最小二乘拟合的方法在Matlab中，可以使用“lsqcurvefit”函数进行对数最小二乘拟合。
该函数可以通过最小二乘法拟合非线性方程，并返回拟合参数和拟合效果。
使用该函数时，需要提供拟合的非线性方程、初始参数估计值、实验数据及其权重等信息，以便进行拟合。
在拟合结束后，可以通过绘制拟合曲线和残差图来评估拟合效果。
四、对数最小二乘拟合的应用对数最小二乘拟合在实际应用中具有广泛的意义。
在生物学领域，用对数最小二乘拟合可以研究物种种裙的增长趋势；在经济学领域，可以用来分析经济指标的增长规律等。
通过对数最小二乘拟合，研究人员可以更加准确地描述实验数据的变化规律，从而做出更有力的数据分析和预测。
五、结论对数最小二乘拟合是Matlab中的重要功能之一，它可以帮助研究人员处理实验数据、分析数据规律，并进行数学建模。
本文讨论了对数最小二乘拟合的原理、方法和应用，希望可以为对数最小二乘拟合的研究和应用提供一些有益的参考。
在实际应用中，研究人员可以根据具体的问题和实验数据，灵活运用Matlab中的数学工具，进行数据分析和建模工作。
六、对数最小二乘拟合的优缺点尽管对数最小二乘拟合在处理指数增长或指数衰减的数据时具有一定优势，但同样也存在一些局限性。
学院：光电信息学院 姓名：赵海峰 学号：1001一、曲线拟合的最小二乘法原理：由已知的离散数据点选择与实验点误差最小的曲线)(...)()()(1100x a x a x a x S n n ϕϕϕ+++=称为曲线拟合的最小二乘法。
若记),()()(),(0i k i j mi i k j x x x ϕϕωϕϕ∑==k i k i mi i k d x x f x f ≡=∑=)()()(),(0ϕωϕ上式可改写为),...,1,0(;),(n k d a k j noj j k -=∑=ϕϕ这个方程成为法方程，可写成距阵形式d Ga =其中,),...,,(,),...,,(1010T n T n d d d d a a a a ==⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=),(),(),()(),(),(),(),(),(101110101000n n n n n n G ϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕΛM M M ΛΛ。
它的平方误差为：.)]()([)(||||2022i i mi i x f x S x -=∑=ωδ二、数值实例：下面给定的是乌鲁木齐最近1个月早晨7:00左右(新疆时间)的天气预报所得到的温度数据表，按照数据找出任意次曲线拟合方程和它的图像。
下面应用Matlab编程对上述数据进行最小二乘拟合三、Matlab程序代码：x=[1:1:30];y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9,7,6,5,3,1 ];a1=polyfit(x,y,3) %三次多项式拟合%a2= polyfit(x,y,9) %九次多项式拟合%a3= polyfit(x,y,15) %十五次多项式拟合%b1=polyval(a1,x)b2=polyval(a2,x)b3=polyval(a3,x)r1= sum((y-b1).^2) %三次多项式误差平方和%r2= sum((y-b2).^2) %九次次多项式误差平方和%r3= sum((y-b3).^2) %十五次多项式误差平方和%plot(x,y,'*') %用*画出x,y图像%hold onplot(x,b1, 'r') %用红色线画出x,b1图像%hold onplot(x,b2, 'g') %用绿色线画出x,b2图像%hold onplot(x,b3, 'b:o') %用蓝色o线画出x,b3图像%四、数值结果：不同次数多项式拟和误差平方和为：r1 =r2 =r3 =r1、r2、r3分别表示三次、九次、十五次多项式误差平方和。
1 曲线拟合的线性最小二乘法及其MATLAB 程序例7.2.1 给出一组数据点),(i i y x 列入表7–2中，试用线性最小二乘法求拟合曲线，并用（7.2），（7.3）和（7.4）式估计其误差，作出拟合曲线.表7–2 例7.2.1的一组数据),(y x解 （1）在MATLAB 工作窗口输入程序>> x=[-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6];y=[-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04];plot(x,y,'r*'),legend('实验数据(xi,yi)')xlabel('x'), ylabel('y'),title('例7.2.1的数据点(xi,yi)的散点图')运行后屏幕显示数据的散点图（略）.（3）编写下列MA TLAB 程序计算)(x f 在),(i i y x 处的函数值，即输入程序>> syms a1 a2 a3 a4x=[-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6];fi=a1.*x.^3+ a2.*x.^2+ a3.*x+ a4运行后屏幕显示关于a 1,a 2, a 3和a 4的线性方程组fi =[ -125/8*a1+25/4*a2-5/2*a3+a4,-4913/1000*a1+289/100*a2-17/10*a3+a4,-1331/1000*a1+121/100*a2-11/10*a3+a4,-64/125*a1+16/25*a2-4/5*a3+a4,a4, 1/1000*a1+1/100*a2+1/10*a3+a4,27/8*a1+9/4*a2+3/2*a3+a4, 19683/1000*a1+729/100*a2+27/10*a3+a4, 5832/125*a1+324/25*a2+18/5*a3+a4]编写构造误差平方和的MATLAB 程序>> y=[-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.5068.04];fi=[-125/8*a1+25/4*a2-5/2*a3+a4,-4913/1000*a1+289/100*a2-17/10*a3+a4,-1331/1000*a1+121/100*a2-11/10*a3+a4,-64/125*a1+16/25*a2-4/5*a3+a4, a4, 1/1000*a1+1/100*a2+1/10*a3+a4,27/8*a1+9/4*a2+3/2*a3+a4,19683/1000*a1+729/100*a2+27/10*a3+a4,5832/125*a1+324/25*a2+18/5*a3+a4];fy=fi-y; fy2=fy.^2; J=sum(fy.^2)运行后屏幕显示误差平方和如下J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)^2+(-4913/1000*a1+289/100*a2-17/10*a3+a4+171/2)^2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20)^2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)^2+(a4+91/10)^2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)^2+(27/8*a1+9/4*a 2+3/2*a3+a4+328/25)^2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/2)^2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)^2为求4321,,,a a a a 使J 达到最小，只需利用极值的必要条件0=∂∂ka J )4,3,2,1(=k ，得到关于4321,,,a a a a 的线性方程组，这可以由下面的MA TLAB 程序完成，即输入程序>> syms a1 a2 a3 a4J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)^2+(-4913/1000*a1+289/100*a2-17/10*a3+a4...+171/2)^2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20)^2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)^2+(a 4+91/10)^2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)^2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)^2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/2)^2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)^2;Ja1=diff(J,a1); Ja2=diff(J,a2); Ja3=diff(J,a3); Ja4=diff(J,a4);Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3), Ja41=simple(Ja4),运行后屏幕显示J 分别对a 1, a 2 ,a 3 ,a 4的偏导数如下Ja11=56918107/10000*a1+32097579/25000*a2+1377283/2500*a3+23667/250*a4-8442429/625Ja21 =32097579/25000*a1+1377283/2500*a2+23667/250*a3+67*a4+767319/625Ja31 =1377283/2500*a1+23667/250*a2+67*a3+18/5*a4-232638/125Ja41 =23667/250*a1+67*a2+18/5*a3+18*a4+14859/25解线性方程组Ja 11 =0，Ja 21 =0，Ja 31 =0，Ja 41 =0，输入下列程序>>A=[56918107/10000, 32097579/25000, 1377283/2500, 23667/250; 32097579/25000, 1377283/2500, 23667/250, 67; 1377283/2500, 23667/250, 67, 18/5; 23667/250, 67, 18/5, 18];B=[8442429/625, -767319/625, 232638/125, -14859/25];C=B/A, f=poly2sym(C)运行后屏幕显示拟合函数f 及其系数C 如下C = 5.0911 -14.1905 6.4102 -8.2574f=716503695845759/140737488355328*x^3-7988544102557579/562949953421312*x^2+1804307491277693/281474976710656*x-4648521160813215/562949953421312故所求的拟合曲线为8.25746.410214.19055.0911)(23-+-=x x x x f .（4）编写下面的MA TLAB 程序估计其误差，并作出拟合曲线和数据的图形.输入程序>> xi=[-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6];y=[-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04];n=length(xi);f=5.0911.*xi.^3-14.1905.*xi.^2+6.4102.*xi -8.2574;x=-2.5:0.01: 3.6;F=5.0911.*x.^3-14.1905.*x.^2+6.4102.*x -8.2574;fy=abs(f-y); fy2=fy.^2; Ew=max(fy),E1=sum(fy)/n, E2=sqrt((sum(fy2))/n)plot(xi,y,'r*'), hold on, plot(x,F,'b-'), hold offlegend('数据点(xi,yi)','拟合曲线y=f(x)'),xlabel('x'), ylabel('y'),title('例7.2.1的数据点(xi,yi)和拟合曲线y=f(x)的图形')运行后屏幕显示数据),(i i y x 与拟合函数f 的最大误差E w ，平均误差E 1和均方根误差E 2及其数据点),(i i y x 和拟合曲线y =f (x )的图形（略）.Ew = E1 = E2 =3.105 4 0.903 4 1.240 97.3 函数)(x r k 的选取及其MATLAB 程序例7.3.1 给出一组实验数据点),(i i y x 的横坐标向量为x =(-8.5,-8.7,-7.1,-6.8,-5.10,-4.5,-3.6,-3.4,-2.6,-2.5, -2.1,-1.5, -2.7,-3.6)，纵横坐标向量为y =(459.26,52.81,198.27,165.60,59.17,41.66,25.92, 22.37,13.47, 12.87, 11.87,6.69,14.87,24.22)，试用线性最小二乘法求拟合曲线，并用（7.2），（7.3）和（7.4）式估计其误差，作出拟合曲线.解 （1）在MATLAB 工作窗口输入程序>>x=[-8.5,-8.7,-7.1,-6.8,-5.10,-4.5,-3.6,-3.4,-2.6,-2.5,-2.1,-1.5, -2.7,-3.6];y=[459.26,52.81,198.27,165.60,59.17,41.66,25.92,22.37,13.47, 12.87, 11.87,6.69,14.87,24.22];plot(x,y,'r*'),legend('实验数据(xi,yi)')xlabel('x'), ylabel('y'),title('例7.3.1的数据点(xi,yi)的散点图')运行后屏幕显示数据的散点图（略）.（3）编写下列MA TLAB 程序计算)(x f 在),(i i y x 处的函数值，即输入程序>> syms a bx=[-8.5,-8.7,-7.1,-6.8,-5.10,-4.5,-3.6,-3.4,-2.6,-2.5,-2.1,-1.5,-2.7,-3.6]; fi=a.*exp(-b.*x)运行后屏幕显示关于a 和b 的线性方程组fi =[ a*exp(17/2*b), a*exp(87/10*b), a*exp(71/10*b),a*exp(34/5*b), a*exp(51/10*b), a*exp(9/2*b), a*exp(18/5*b), a*exp(17/5*b), a*exp(13/5*b), a*exp(5/2*b), a*exp(21/10*b), a*exp(3/2*b), a*exp(27/10*b), a*exp(18/5*b)]编写构造误差平方和的MATLAB 程序如下>>y=[459.26,52.81,198.27,165.60,59.17,41.66,25.92,22.37,13.47,12.87, 11.87, 6.69,14.87,24.22];fi =[ a*exp(17/2*b), a*exp(87/10*b), a*exp(71/10*b), a*exp(34/5*b), a*exp(51/10*b), a*exp(9/2*b), a*exp(18/5*b), a*exp(17/5*b), a*exp(13/5*b), a*exp(5/2*b), a*exp(21/10*b), a*exp(3/2*b), a*exp(27/10*b), a*exp(18/5*b)];fy=fi-y;fy2=fy.^2;J=sum(fy.^2)运行后屏幕显示误差平方和如下J =(a*exp(17/2*b)-22963/50)^2+(a*exp(87/10*b)-5281/100)^2+(a*exp(71/10*b)-19827/100)^2+(a*exp(34/5*b)-828/5)^2+(a*exp(51/10*b)-5917/100)^2+(a*exp(9/2*b)-2083/50)^2+(a*exp(18/5*b)-648/25)^2+(a*exp(17/5*b)-2237/100)^2+(a*exp(13/5*b)-1347/100)^2+(a*ex p(5/2*b)-1287/100)^2+(a*exp(21/10*b)-1187/100)^2+(a*exp(3/2*b)-669/100)^2+(a*exp(27/10*b)-1487/100)^2+(a*exp(18/5*b)-1211/50)^2为求b a ,使J 达到最小，只需利用极值的必要条件，得到关于b a ,的线性方程组，这可以由下面的MA TLAB 程序完成，即输入程序>> syms a bJ=(a*exp(17/2*b)-22963/50)^2+(a*exp(87/10*b)-5281/100)^2+(a*exp(71/10*b)-19827/100)^2+(a*exp(34/5*b)-828/5)^2+(a*exp(51/10*b)-5917/100)^2+(a*exp(9/2*b)-2083/50)^2+(a*exp(18/5*b)-648/25)^2+(a*exp(17/5*b)-2237/100)^2+(a*exp(13/5*b)-1347/100)^2+(a*exp(5/2*b)-1287/100)^2+(a*exp(21/10*b)-1187/100)^2+(a*exp(3/2*b )-669/100)^2+(a*exp(27/10*b)-1487/100)^2+(a*exp(18/5*b)-1211/50)^2;Ja=diff(J,a); Jb=diff(J,b);Ja1=simple(Ja), Jb1=simple(Jb),运行后屏幕显示J 分别对b a ,的偏导数如下Ja1 =2*a*exp(3*b)+2*a*exp(17*b)+2*a*exp(87/5*b)+2*exp(68/5*b)*a+2*exp(9*b)*a+2*a*exp(34/5*b)-669/50*exp(3/2*b)-1487/50*exp(27/10*b)-2507/25*exp(18/5*b)-22963/25*exp(17/2*b)-5281/50*exp(87/10*b)-19827/50*exp(71/10*b)-2237/50*exp(17/5*b)-1656/5*exp(34/5*b)-1347/50*exp(13/5*b)-5917/50*exp(51/10*b)-1287/50*exp(5/2*b )-2083/25*exp(9/2*b)-1187/50*exp(21/10*b)+4*a*exp(36/5*b)+2*a*e xp(26/5*b)+2*a*exp(71/5*b)+2*a*exp(51/5*b)+2*a*exp(5*b)+2*a*exp (21/5*b)+2*a*exp(27/5*b)Jb1 =1/500*a*(2100*a*exp(21/10*b)^2+8500*a*exp(17/2*b)^2+6800*a*exp(34/5*b)^2-10035*exp(3/2*b)-40149*exp(27/10*b)-180504*exp (18/5*b)-3903710*exp(17/2*b)-459447*exp(87/10*b)-1407717*exp(71/10*b)-76058*exp(17/5*b)-1126080*exp(34/5*b)-35022*exp(13/5*b)-301767*exp(51/10*b)-32175*exp(5/2*b)-187470*exp(9/2*b)-24927*ex p(21/10*b)+7100*a*exp(71/10*b)^2+5100*a*exp(51/10*b)^2+4500*a*e xp(9/2*b)^2+7200*a*exp(18/5*b)^2+3400*a*exp(17/5*b)^2+2600*a*ex p(13/5*b)^2+2500*a*exp(5/2*b)^2+1500*a*exp(3/2*b)^2+2700*a*exp(27/10*b)^2+8700*a*exp(87/10*b)^2)用解二元非线性方程组的牛顿法的MATLAB 程序求解线性方程组J a1 =0，J b1 =0，得a = b=2.811 0 0.581 6故所求的拟合曲线（7.13）为0811.2)(=x f e x 5816.0-. （7.14）（4）根据（7.2），（7.3），（7.4）和（7.14）式编写下面的MA TLAB 程序估计其误差，并做出拟合曲线和数据的图形.输入程序>> xi=[-8.5 -8.7 -7.1 -6.8 -5.10 -4.5 -3.6 -3.4 -2.6 -2.5-2.1 -1.5 -2.7 -3.6];y=[459.26 52.81 198.27 165.60 59.17 41.66 25.92 22.3713.47 12.87 11.87 6.69 14.87 24.22];n=length(xi); f=2.8110.*exp(-0.5816.*xi); x=-9:0.01: -1;F=2.8110.*exp(-0.5816.*x); fy=abs(f-y); fy2=fy.^2;Ew=max(fy),E1=sum(fy)/n, E2=sqrt((sum(fy2))/n), plot(xi,y,'r*'), hold on plot(x,F,'b-'), hold off,legend('数据点(xi,yi)','拟合曲线y=f(x)')xlabel('x'), ylabel('y'),title('例7.3.1的数据点(xi,yi)和拟合曲线y=f(x)的图形')运行后屏幕显示数据),(i i y x 与拟合函数f 的最大误差E w = 390.141 5，平均误差E 1=36.942 2和均方根误差E 2=106.031 7及其数据点),(i i y x 和拟合曲线y =f (x )的图形（略）.7.4 多项式拟合及其MATLAB 程序例7.4.1 给出一组数据点),(i i y x 列入表7–3中，试用线性最小二乘法求拟合曲线，并用（7.2），（7.3）和（7.4）式估计其误差，作出拟合曲线.表7–3 例7.4.1的一组数据),(y x解 （1）首先根据表7–3给出的数据点i i ，用下列MATLAB 程序画出散点图.在MATLAB 工作窗口输入程序>> x=[-2.9 -1.9 -1.1 -0.8 0 0.1 1.5 2.7 3.6];y=[53.94 33.68 20.88 16.92 8.79 8.98 4.17 9.1219.88];plot(x,y,'r*'), legend('数据点(xi,yi)')xlabel('x'), ylabel('y'),title('例7.4.1的数据点(xi,yi)的散点图')运行后屏幕显示数据的散点图（略）.（3）用作线性最小二乘拟合的多项式拟合的MATLAB 程序求待定系数k a )3,2,1(=k .输入程序>> a=polyfit(x,y,2)运行后输出（7.16）式的系数a =2.8302 -7.3721 9.1382故拟合多项式为2138.91372.72830.2)(2+-=x x x f .（4）编写下面的MATLAB 程序估计其误差，并做出拟合曲线和数据的图形.输入程序>> xi=[-2.9 -1.9 -1.1 -0.8 0 0.1 1.5 2.7 3.6];y=[53.94 33.68 20.88 16.92 8.79 8.98 4.17 9.12 19.88];n=length(xi); f=2.8302.*xi.^2-7.3721.*xi+9.1382x=-2.9:0.001:3.6;F=2.8302.*x.^2-7.3721.*x+8.79;fy=abs(f-y); fy2=fy.^2; Ew=max(fy), E1=sum(fy)/n,E2=sqrt((sum(fy2))/n), plot(xi,y,'r*', x,F,'b-'),legend('数据点(xi,yi)','拟合曲线y=f(x)')xlabel('x'), ylabel('y'),title('例7.4.1 的数据点(xi,yi)和拟合曲线y=f(x)的图形')运行后屏幕显示数据),(i i y x 与拟合函数f 的最大误差E w ，平均误差E1和均方根误差E 2及其数据点(x i ,y i )和拟合曲线y =f (x )的图形（略）.Ew = E1 = E2 =0.745 7, 0.389 2, 0.436 37.5 拟合曲线的线性变换及其MATLAB 程序例7.5.1 给出一组实验数据点),(i i y x 的横坐标向量为x =（7.5 6.8 5.10 4.53.6 3.4 2.6 2.5 2.1 1.5 2.7 3.6），纵横坐标向量为y =（359.26 165.60 59.17 41.66 25.92 22.37 13.47 12.87 11.87 6.69 14.87 24.22），试用线性变换和线性最小二乘法求拟合曲线，并用（7.2），（7.3）和（7.4）式估计其误差，作出拟合曲线.解 （1）首先根据给出的数据点),(i i y x ，用下列MATLAB 程序画出散点图.在MATLAB 工作窗口输入程序>> x=[7.5 6.8 5.10 4.5 3.6 3.4 2.6 2.5 2.1 1.5 2.73.6];y=[359.26 165.60 59.17 41.66 25.92 22.37 13.47 12.87 11.87 6.69 14.87 24.22];plot(x,y,'r*'), legend('数据点(xi,yi)')xlabel('x'), ylabel('y'),title('例7.5.1的数据点(xi,yi)的散点图')运行后屏幕显示数据的散点图（略）.（2）根据数据散点图，取拟合曲线为a y =e bx )0,0(≠>b a , （7.19）其中b a ,是待定系数.令b B a A y Y ===,ln ,ln ，则（7.19）化为Bx A Y +=.在MATLAB 工作窗口输入程序>> x=[7.5 6.8 5.10 4.5 3.6 3.4 2.6 2.5 2.1 1.5 2.73.6];y=[359.26 165.60 59.17 41.66 25.92 22.37 13.47 12.87 11.87 6.69 14.87 24.22];Y=log(y); a=polyfit(x,Y,1); B=a(1);A=a(2); b=B,a=exp(A)n=length(x); X=8:-0.01:1; Y=a*exp(b.*X); f=a*exp(b.*x);plot(x,y,'r*',X,Y,'b-'), xlabel('x'),ylabel('y')legend('数据点(xi,yi)','拟合曲线y=f(x)')title('例7.5.1 的数据点(xi,yi)和拟合曲线y=f(x)的图形')fy=abs(f-y); fy2=fy.^2; Ew=max(fy), E1=sum(fy)/n,E2=sqrt((sum(fy2))/n)运行后屏幕显示a y =e bx 的系数b =0.624 1，a =2.703 9,数据),(i i y x 与拟合函数f 的最大误差Ew =67.641 9，平均误差E 1=8.677 6和均方根误差E 2=20.711 3及其数据点),(i i y x 和拟合曲线9703.2)(=x f e x 1624.0的图形（略）.7.6 函数逼近及其MATLAB 程序最佳均方逼近的MATLAB 主程序function [yy1,a,WE]=zjjfbj(f,X,Y,xx)m=size(f);n=length(X);m=m(1);b=zeros(m,m); c=zeros(m,1);if n~=length(Y)error('X 和Y 的维数应该相同')endfor j=1:mfor k=1:mb(j,k)=0;for i=1:nb(j,k)=b(j,k)+feval(f(j,:),X(i))*feval(f(k,:),X(i));endendc(j)=0;for i=1:nc(j)=c(j)+feval(f(j,:),X(i))*Y(i);endenda=b\c;WE=0;for i=1:nff=0;for j=1:mff=ff+a(j)*feval(f(j,:),X(i));endWE=WE+(Y(i)-ff)*(Y(i)-ff);endif nargin==3return ;endyy=[];for i=1:ml=[];for j=1:length(xx)l=[l,feval(f(i,:),xx(j))];endyy=[yy l'];endyy=yy*a; yy1=yy'; a=a';WE;例7.6.1 对数据X 和Y , 用函数2,,1x y x y y ===进行逼近，用所得到的逼近函数计算在 6.5=x 处的函数值，并估计误差.其中X =(1 3 4 5 6 7 8 9); Y =(-11 -13 -11 -7 -1 7 17 29).解 在MATLAB 工作窗口输入程序>> X=[ 1 3 4 5 6 7 8 9]; Y=[-11 -13 -11 -7 -1 7 17 29];f=['fun0';'fun1';'fun2']; [yy,a,WE]=zjjfbj(f,X,Y,6.5)运行后屏幕显示如下yy =2.75000000000003a =-7.00000000000010 -4.99999999999995 1.00000000000000WE =7.172323350269439e-027例7.6.2 对数据X 和Y ，用函数2,,1x y x y y ===，x y cos =，=y e x ，x y sin =进行逼近，其中X =（0 0.50 1.00 1.50 2.00 2.50 3.00），Y =（0 0.4794 0.8415 0.9815 0.9126 0.5985 0.1645）.解 在MATLAB 工作窗口输入程序>> X=[ 0 0.50 1.00 1.50 2.00 2.50 3.00];Y=[0 0.4794 0.8415 0.9815 0.9126 0.5985 0.1645];f=['fun0';'fun1';'fun2';'fun3';'fun4';'fun5'];xx=0:0.2:3;[yy,a,WE]=zjjfbj(f,X,Y, xx), plot(X,Y,'ro',xx,yy,'b-')运行后屏幕显示如下（图略）yy = Columns 1 through 7-0.0005 0.2037 0.3939 0.5656 0.7141 0.8348 0.9236Columns 8 through 140.9771 0.9926 0.9691 0.9069 0.8080 0.6766 0.5191Columns 15 through 160.3444 0.1642a = 0.3828 0.4070 -0.3901 0.0765 -0.4598 0.5653 WE = 1.5769e-004即，最佳逼近函数为y=0.3828+0.4070*x-0.3901*x^2+0.0765*exp(x) -0.4598*cos(x) +0.5653*sin(x). 7.7 三角多项式逼近及其MATLAB 程序计算三角多项式的MATLAB 主程序function [A,B,Y1,Rm]=sanjiao(X,Y,X1,m)n= length(X)-1;max1=fix((n-1)/2);if m > max1m=max1;endA=zeros(1,m+1);B=zeros(1,m+1);Ym=(Y(1)+Y(n+1))/2; Y(1)=Ym; Y(n+1)=Ym; A(1)=2*sum(Y)/n;for i=1:mB(i+1)=sin(i*X)*Y'; A(i+1)=cos(i*X)*Y';endA=2*A/n; B=2*B/n; A(1)=A(1)/2;Y1=A(1);for k=1:mY1=Y1+A(k+1)*cos(k*X1)+ B(k+1)*sin(k*X1);Tm=A(1)+A(k+1).*cos(k*X)+ B(k+1).*sin(k*X); k=k+1;endY;Tm; Rm=(sum(Y-Tm).^2)/n;例7.7.1 根据],[ππ-上的350,60,13=n 个等距横坐标点ni x i π+π-=2),,2,1,0(n i =和函数3sin 2)(x x f =. （1）求)(x f 的6阶三角多项式逼近，计算均方误差;（2）将这三个三角多项式分别与)(x f 的傅里叶级数nx n n x f n n sin 19)1(318)(121∑∞=+--π= 的前6项进行比较;（3）利用三角多项式分别计算X i = -2, 2.5的值;（4）在同一坐标系中，画出函数)(x f ，350,60,13=n 的三角多项式和数据点的图形.解 （1）输入程序>> X1=-pi:2*pi/13:pi;Y1=2*sin(X1/3);X1i=[-2,2.5];[A1,B1,Y11,Rm1]=sanjiao(X1,Y1,X1i,6),X2=-pi:2*pi/60:pi;Y2=2*sin(X2/3);[A2,B2,Y12,Rm2]=sanjiao(X2,Y2,X1i,6)X3=-pi:2*pi/350:pi;Y3=2*sin(X3/3);[A3,B3,Y13,Rm3]=sanjiao(X3,Y3,X1i,6)X1i=[-2,2.5];Y1=2*sin(X1i/3)for n=1:6bi=(-1)^(n+1)*18*sqrt(3)*n/(pi*(9*n^2-1))end（2）画图，输入程序>>X1=-pi:2*pi/13:pi;Y1=2*sin(X1/3);Xi=-pi:0.001:pi; f=2*sin(Xi/3);[A1,B1,Y1i,R1m]=sanjiao(X1,Y1,Xi,6);X2=-pi:2*pi/60:pi; Y2=2*sin(X2/3); X3=-pi:2*pi/350:pi;Y3=2*sin(X3/3);[A2,B2,Y2i,R2m]=sanjiao(X2,Y2,Xi,6);[A3,B3,Y3i,R3m]=sanjiao(X3,Y3,Xi,6);plot(X1,Y1,'r*', Xi, Y1i,'b-',Xi, Y2i,'g--', Xi, Y3i, 'm:',Xi, f, 'k-.')xlabel('x'),ylabel('y')legend('数据点(xi,yi)','n=13的三角多项式','n=60的三角多项式','n=350的三角多项式','函数f(x)')title('例7.7.1 的数据点(xi,yi)、n=13，60，350的三角多项式T3和函数f(x)的图形')运行后图形（略）.7.8 随机数据点上的二元拟合及其MATLAB 程序例7.8.1 设节点(X,Y,Z )中的X 和Y 分别是在区间]3,3[-和]5.3,5.2[-上的50个随机数，Z 是函数Z =7-3x 3e22y - -x 在(X,Y )的值，拟合点(X I ,Y I )中的X I =-3:0.2:3, Y I =-2.5:0.2:3.5.分别用二元拟合方法中最近邻内插法、三角基线性内插法、三角基三次内插法和MATLAB 4网格化坐标方法计算在(X I ,Y I )处的值，作出它们的图形，并与被拟和曲面进行比较.解 （1）最近邻内插法.输入程序>> x=rand(50,1);y=rand(50,1); %生成50个一元均匀分布随机数x 和y ， x ，y .X=-3+(3-(-3))*x;%利用x 生成的随机变量.Y=-2.5+(3.5-(-2.5))*y; %利用y 生成的随机变量.Z=7-3* X.^3 .* exp(-X.^2 - Y.^2); %在每个随机点（X,Y ）处计算Z 的值.X1=-3:0.2:3;Y1=-2.5:0.2:3.5;[XI,YI] = meshgrid(X1,Y1); %将坐标（XI,YI ）网格化.ZI=griddata(X,Y,Z,XI,YI, 'nearest') %计算在每个插值点（XI,YI ）处的插值ZI.mesh(XI,YI, ZI) %作二元拟合图形.xlabel('x'), ylabel('y'), zlabel('z'),title('用最近邻内插法拟合函数z =7-3 x^3 exp(-x^2 - y^2) 的曲面和节点的图形')%legend('拟合曲面','节点(xi,yi,zi)')hold on %在当前图形上添加新图形.plot3(X,Y,Z, 'bo') %用兰色小圆圈画出每个节点(X,Y,Z). hold of %结束在当前图形上添加新图形.运行后屏幕显示用最近邻内插法拟合函数Z =7-3x 3e 22y - -x 在两组不同节点处的曲面及其插值Z I （略）.（2）三角基线性内插法.输入程序>> x=rand(50,1);y=rand(50,1); %生成50个一元均匀分布随机数x 和y ， x ，y .X=-3+(3-(-3))*x;%利用x 生成 上的随机变量.Y=-2.5+(3.5-(-2.5))*y; %利用y 生成 上的随机变量.Z=7-3* X.^3 .* exp(-X.^2 - Y.^2); %在每个随机点（X,Y ）处计算Z 的值.X1=-3:0.2:3;Y1=-2.5:0.2:3.5;[XI,YI] = meshgrid(X1,Y1); %将坐标（XI,YI ）网格化.ZI=griddata(X,Y,Z,XI,YI, 'linear') %计算在每个插值点（XI,YI ）处的插值ZI.mesh(XI,YI, ZI) %作二元拟合图形.xlabel('x'), ylabel('y'), zlabel('z'),title('用三角基线性内插法拟合函数z =7-3 x^3 exp(-x^2 - y^2) 的曲面和节点的图形')%legend('拟合曲面','节点(xi,yi,zi)')hold on %在当前图形上添加新图形.plot3(X,Y,Z, 'bo') %用兰色小圆圈画出每个节点(X,Y,Z). hold of %结束在当前图形上添加新图形.运行后屏幕显示用三角基线性内插法拟合函数Z =7-3x 3e 22y - -x 在两组不同节点处的曲面和节点的图形及其插值Z I （略）.（3）三角基三次内插法.输入程序>> x=rand(50,1);y=rand(50,1); %生成50个一元均匀分布随机数x 和y ， x ，y .X=-3+(3-(-3))*x;%利用x 生成 上的随机变量.Y=-2.5+(3.5-(-2.5))*y; %利用y 生成 上的随机变量.Z=7-3* X.^3 .* exp(-X.^2 - Y.^2); %在每个随机点（X,Y ）处计算Z 的值.X1=-3:0.2:3;Y1=-2.5:0.2:3.5;[XI,YI] = meshgrid(X1,Y1); %将坐标（XI,YI ）网格化.ZI=griddata(X,Y,Z,XI,YI, 'cubic') %计算在每个插值点（XI,YI ）处的插值ZI.mesh(XI,YI, ZI) %作二元拟合图形.xlabel('x'), ylabel('y'), zlabel('z'),title('用三角基三次内插法拟合函数z =7-3 x^3 exp(-x^2 - y^2) 的曲面和节点的图形')%legend('拟合曲面','节点(xi,yi,zi)')hold on %在当前图形上添加新图形.plot3(X,Y,Z, 'bo') %用兰色小圆圈画出每个节点(X,Y,Z). hold of %结束在当前图形上添加新图形.运行后屏幕显示用三角基三次内插法拟合函数Z =7-3x 3e 22y - -x 在两组不同节点处的曲面和节点的图形及其插值Z I （略）.（4）MATLAB 4网格化坐标方法.输入程序>> x=rand(50,1);y=rand(50,1); %生成50个一元均匀分布随机数x 和y ， x ，y .X=-3+(3-(-3))*x;%利用x 生成 上的随机变量.Y=-2.5+(3.5-(-2.5))*y; %利用y 生成 上的随机变量.Z=7-3* X.^3 .* exp(-X.^2 - Y.^2); %在每个随机点（X,Y ）处计算Z 的值.X1=-3:0.2:3; Y1=-2.5:0.2:3.5;[XI,YI] = meshgrid(X1,Y1); %将坐标（XI,YI ）网格化.ZI=griddata(X,Y,Z,XI,YI, 'v4') %计算在每个插值点（XI,YI ）处的插值ZI.mesh(XI,YI, ZI) %作二元拟合图形.xlabel('x'), ylabel('y'), zlabel('z'),title('用MATLAB 4网格化坐标方法拟合函数z =7-3 x^3 exp(-x^2 - y^2) 的曲面和节点的图形')%legend('拟合曲面','节点(xi,yi,zi)')hold on %在当前图形上添加新图形.plot3(X,Y,Z, 'bo') %用兰色小圆圈画出每个节点(X,Y,Z). hold of %结束在当前图形上添加新图形.运行后屏幕显示用MATLAB 4网格化坐标方法拟合函数Z =7-3x 3e22y - -x 在两组不同节点处的曲面和节点的图形及其插值ZI （略）.（5）作被拟合曲面Z =7-3x 3e 22y - -x 和节点的图形.输入程序>> x=rand(50,1);y=rand(50,1); %生成50个一元均匀分布随机数x 和y ， x ，y .X=-3+(3-(-3))*x;%利用x 生成随机变量.Y=-2.5+(3.5-(-2.5))*y; %利用y 生成随机变量.Z=7-3* X.^3 .* exp(-X.^2 - Y.^2); %在每个随机点（X,Y ）处计算Z 的值.X1=-3.:0.1:3.;Y1=-2.5:0.1:3.5;[XI,YI] = meshgrid(X1,Y1); %将坐标（XI,YI ）网格化.ZI=7-3* XI.^3 .* exp(-XI.^2 - YI.^2);mesh(XI,YI, ZI) %作二元拟合图形.xlabel('x'), ylabel('y'), zlabel('z'),title('被拟合函数z =7-3 x^3 exp(-x^2 - y^2) 的曲面和节点的图形')%legend('被拟合函数曲面','节点(xi,yi,zi)')hold on %在当前图形上添加新图形.plot3(X,Y,Z, 'bo') %用兰色小圆圈画出每个节点(X,Y,Z). hold of %结束在当前图形上添加新图形.运行后屏幕显示被拟合函数Z =7-3x 3e22y - -x 的曲面和节点的图形及其函数值ZI （略）.7.9 随机数据点上的n 元拟合及其MATLAB 程序例7.9.1 首先利用MATLAB 函数rand 产生随机数据X 1,Y 1,Z 1,然后用线性变换b at u += (其中5,5-==b a )将随机数据X 1,Y 1 ,Z 1 变换为节点坐标(X,Y,Z ),再用函数)1(373+-=z y x w e222z y x ---生成数据W , 用三元最近邻内插法方法计算函数w 在插值点x i ,10:5.0:3-= y i ,13:5.0:2-=z i =y i 处拟合数据的值，并作其图形.解 输入程序>> X1=-5+5*rand(10,1);Y1=-5+5*rand(10,1);Z1=Y1;[X,Y,Z] = meshgrid(X1,Y1,Z1);W=7-3* X.^3 .* Y.*(Z+1).* exp(-X.^2 - Y.^2- Z.^2);xi=-3:0.5:10;yi=-2:0.5:13;zi=yi;[XI,YI,ZI] = meshgrid(xi,yi,zi);W1=griddata3(X, Y, Z, W, XI, YI, ZI, 'nearest');slice(XI,YI,ZI,W1,[-2 4 9.5],9,[-2 2 9]),%shading flat%lighting flatxlabel('x'), ylabel('y'), zlabel('z'),title('被拟合函数W=7-3X^3Y(Z+1)exp(-X^2 - Y^2- Z^2) '); hold oncolorbar('horiz')view([-30 45])运行后屏幕显示三元线性拟合值及其图形（略）.例7.9.2 设节点(X,Y,Z,W )中的X,Y 和Z 分别是在区间]3,3[-和]5.3,5.2[-,Y =Z 上的15个随机数，W 是函数x w +=2e 222z y x ---在(X,Y,Z )的值，拟合点(x i ,y i ,z i )中的x i =-3:0.2:3， y i =-2.5:0.2:3.5，z i =y i ， 用'linear '方法计算拟合数据的值，并作其图形.解 输入程序>> x=rand(15,1); y=rand(15,1);X1=-3+(3-(-3))*x;Y1=-2.5+(3.5-(-2.5))*y;Z1=Y1;[X,Y,Z] = meshgrid(X1,Y1,Z1);W=2+X.* exp(-X.^2 - Y.^2- Z.^2);xi=-3:0.2:3; yi=-2.5:0.2:3.5; zi=yi;[X2,Y2,Z2]=meshgrid(xi,yi,zi);W1=griddata3(X, Y, Z, W, X2,Y2,Z2,'linear');slice(X2,Y2,Z2,W1,[-1 0 1.5],2,[-2 3]),shading flat,lighting flat,xlabel('x'), ylabel('y'), zlabel('z'),title('被拟合函数W=2+X exp(-X^2 - Y^2- Z^2)');hold on,colorbar('horiz'), view([-3 5])运行后屏幕显示三元线性拟合值及其图形（略）.。
---------------------------------------------------------------最新资料推荐------------------------------------------------------ 最小二乘法拟合圆公式推导及matlab实现20 09 -01 -1 7 | 最小二乘法拟合圆公式推导及 m at lab 实现最小二乘法(least squares analysis)是一种数学优化技术，它通过最小化误差的平方和找到一组数据的最佳函数匹配。
最小二乘法是用最简的方法求得一些绝对不可知的真值，而令误差平方之和为最小。
最小二乘法通常用于曲线拟合 (least squares fitting) 。
这里有拟合圆曲线的公式推导过程和 vc 实现。
m at lab 实现：function [ R,A,B] = circ(x,y,N) x1 = 0; x2 = 0; x3= 0; y1 = 0; y2 = 0; y3 = 0; x1 y1 = 0; x1 y2 =0; x2y1 = 0; for i = 1 : N x1 = x1 + x(i); x2 = x2 + x(i)* x(i); x3 = x3 + x(i) * x(i)* x(i); y1 = y1 + y(i); y2 = y2 + y(i)* y(i); y3 =y3 + y(i) * y(i)* y(i); x1 y1 = x1 y1 + x(i) * y(i) ; x1 y2 = x1 y2 + x(i)* y(i)* y(i); x2y1 = x2y1 +x(i)* x(i)* y(i); end C = N * x2 - x1 * x1 ; D = N * x1y1 - x1 * y1 ; E = N * x3 + N * x1 y2 - (x2 + y2) *x1 ; G = N * y2 - y1 * y1 ; H = N * x2y1 + N * y3- (x2 + y2) * y1 ; a = (H * D - E * G) /(C * G - D * D); b = (H * C - E * D)/(D * D - G * C) ; c = -(a * x11 / 3+ b * y1 + x2 + y2)/N; A = a/(-2) ; % x 坐标 B = b/(-2) ; % y 坐标 R = sqrt(a * a + b * b- 4 * c)/2; void CViewActionImageTool::LeastSquaresFitting() { if(m_nNum3) { return; } int i=0; doubleX1=0; double Y1=0; double X2=0; double Y2=0; double X3=0; double Y3=0; double X1Y1=0; doubleX1Y2=0; double X2Y1=0; for (i=0;im_nNum;i++) { X1 = X1 + m_points[i].x; Y1 = Y1 +m_points[i].y; X2 = X2 + m_points[i].x*m_points[i].x; Y2 = Y2 + m_points[i].y*m_points[i].y; X3 = X3 +m_points[i].x*m_points[i].x*m_points[i].x; Y3 = Y3 +m_points[i].y*m_points[i].y*m_points[i].y; X1Y1 =X1Y1 + m_points[i].x*m_points[i].y; X1Y2 = X1Y2 +m_points[i].x*m_points[i].y*m_points[i].y; X2Y1 =X2Y1 + m_points[i].x*m_points[i].x*m_points[i].y; } double C,D,E,G,H,N; double a,b,c; N = m_nNum; C =N*X2 - X1*X1; D = N*X1Y1 - X1*Y1; E = N*X3 + N*X1Y2 -(X2+Y2)*X1; G = N*Y2 - Y1*Y1; H = N*X2Y1 + N*Y3 -(X2+Y2)*Y1; a = (H*D-E*G)/(C*G-D*D); b =(H*C-E*D)/(D*D-G*C); c = -(a*X1 + b*Y1 + X2 + Y2)/N; double A,B,R; A = a/(-2); B = b/(-2); R =---------------------------------------------------------------最新资料推荐------------------------------------------------------sqrt(a*a+b*b-4*c)/2; m_fCenterX = A; m_fCenterY = B; m_fRadius = R; return;}3 / 3。
曲线拟合( curve-fitting )：工程实践中，用测量到的一些离散的数据{( X, yj,i 0,1,2,...m}求一个近似的函数(x)来拟合这组数据，要求所得的拟合曲线能最好的反映数据的基本趋势(即使(x)最好地逼近f x，而不必满足插值原则。
因此没必要取(X)=y i，只要使i (X i) y i尽可能地小)。
原理：给定数据点{( x i,y i),i 0 ,1 , 2, . . . m} 。
求近似曲线( x) 。
并且使得近似曲线与f x 的偏差最小。
近似曲线在该点处的偏差i(x i ) y i，i=1,2,...,m 。
常见的曲线拟合方法：1. 使偏差绝对值之和最小2. 使偏差绝对值最大的最小3. 使偏差平方和最小最小二乘法：按偏差平方和最小的原则选取拟合曲线，并且采取二项式方程为拟合曲线的方法,称为最小二乘法。
推导过程：1. 设拟合多项式为：2. 各点到这条曲线的距离之和，即偏差平方和如下：3•问题转化为求待定系数a0...a k对等式右边求q偏导数，因而我们得到了：4、把这些等式化简并表示成矩阵的形式，就可以得到下面的矩阵：5. 将这个范德蒙得矩阵化简后可得到：6. 也就是说X*A=Y ，那么A = (X'*X)-1*X'*Y ，便得到了系数矩阵A ，同时，我们也就得到了拟合曲线。
MATLAB 实现：MATLAB 提供了polyfit ()函数命令进行最小二乘曲线拟合。
调用格式：p=polyfit(x,y,n)[p,s]= polyfit(x,y,n)[p,s,mu]=polyfit(x,y,n)x,y 为数据点，n 为多项式阶数，返回p 为幂次从高到低的多项式系数向量p。
x 必须是单调的。
矩阵s包括R (对x进行QR分解的三角元素)、df(自由度)、normr(残差)用于生成预测值的误差估计。
[p,s,mu]=polyfit(x,y,n) 在拟合过程中，首先对x 进行数据标准化处理，以在拟合中消除量纲等影响，mu包含标准化处理过程中使用的x的均值和标准差。
polyval( ) 为多项式曲线求值函数，调用格式：y=polyval(p,x)[y,DELTA]=polyval(p,x,s)y=polyval(p,x)为返回对应自变量x在给定系数P的多项式的值。
[y,DELTA]=polyval(p,x,s) 使用polyfit函数的选项输出s得出误差估计丫 DELTA。
它假设polyfit函数数据输入的误差是独立正态的，并且方差为常数。
则丫 DELTA 将至少包含50%的预测值。
如下给定数据的拟合曲线：x=[0.5,1.0,1.5,2.0,2.5,3.0]，y=[1.75,2.45,3.81,4.80,7.00,8.60]。
解：MATLAB 程序如下：x=[0.5,1.0,1.5,2.0,2.5,3.0];y=[1.75,2.45,3.81,4.80,7.00,8.60];p=polyfit(x,y,2) x1=0.5:0.05:3.0;y1=polyval(p,x1); plot(x,y,'*r',x1,y1,'-b') 运行结果如图1 计算结果为：p =0.5614 0.8287 1.1560 即所得多项式为y=0.5614xA2+0.08287x+1.15560图1 最小二乘法曲线拟合示例对比检验拟合的有效性：例：在[0, n区间上对正弦函数进行拟合，然后在[0,2 n区间画出图形，比较拟合区间和非拟合区间的图形，考察拟合的有效性。
在MATLAB 中输入如下代码：clearx=0:0.1:pi; y=sin(x);[p,mu]=polyfit(x,y,9) x1=0:0.1:2*pi;y1=sin(x1);% 实际曲线y2=polyval(p,x1);%根据由区间0到pi上进行拟合得到的多项式计算0到2pi上的函数值，%需要注意的是polyval ()返回的函数值在pi到2pi上并没有进行拟合plot(x1,y2,'k*',x1,y1,'k-')运行结果：p =0.0000 0.0000 -0.0003 0.0002 0.0080 0.0002 -0.1668 0.0000 1.0000 0.0000 mu =R: [10x10 double]df: 22 normr: 1.6178e-07MATLAB 的最优化工具箱还提供了lsqcurvefit ()函数命令进行最小二乘曲线拟合(Solve non li near curve-fitti ng (data-fitti ng) problems in least-squares sense)。
调用格式：x = lsqcurvefit(fun,x0,xdata,ydata)x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub)x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)x = lsqcurvefit(problem)[x,resnorm] = lsqcurvefit(...)[x,resnorm,residual] = lsqcurvefit(...)[x,resnorm,residual,exitflag] = lsqcurvefit(...)[x,resnorm,residual,exitflag,output] = lsqcurvefit(...) [x,resnorm,residual,exitflag,output,lambda] = lsqcurvefit(...) [x,resnorm,residual,exitflag,output,lambda,jacobian] = x0 为初始解向量;xdata ，ydata 为满足关系ydata=F(x, xdata) 的数据;lb、ub为解向量的下界和上界，若没有指定界，则lb=[ ] , ub=[];options 为指定的优化参数;fun 为拟合函数，其定义方式为:x = lsqcurvefit(@myfun,x0,xdata,ydata) ，其中myfun 已定义为function F = myfun(x,xdata)F =…%计算x处拟合函数值fun的用法与前面相同；resnorm=sum ((fun(x,xdata)-ydataF2)，即在x 处残差的平方和； residual=fun(x,xdata)-ydata ，即在x 处的残差;exitflag 为终止迭代的条件;output 为输出的优化信息;lambda 为解x 处的Lagrange 乘子;jacobian 为解x 处拟合函数fun 的jacobian 矩阵。
例：lsqcurvefit() 优化程序Data = ...[0.0000 5.89550.1000 3.56390.2000 2.51730.3000 1.97900.4000 1.89900.5000 1.39380.6000 1.13590.7000 1.00960.8000 1.03430.9000 0.84351.0000 0.68561.1000 0.61001.2000 0.53921.3000 0.39461.4000 0.39031.5000 0.54741.6000 0.34591.7000 0.13701.8000 0.22111.9000 0.17042.0000 0.2636];t = Data(:,1);y = Data(:,2);% axis([0 2 -0.5 6])plot(t,y,'ro')title('Data points')%We would like to fit the function y = c(1)*exp(-lam(1)*t) + c(2)*exp(-lam(2)*t) to the data %The lsqcurvefit function solves this type of problem easily.%To begin, define the parameters in terms of one variable x:%x(1) = c(1)%x(2) = lam(1)%x(3) = c(2)%x(4) = lam(2)%Then define the curve as a function of the parameters x and the data t:F = @(x,xdata)x(1)*exp(-x(2)*xdata) + x(3)*exp(-x(4)*xdata);x0 = [1 1 1 0];[x,resnorm,~,exitflag,output] = lsqcurvefit(F,x0,t,y)hold onplot(t,F(x,t))hold offFsumsquares = @(x)sum((F(x,t) - y).A2);opts = optimset('LargeScale','off');[xunc,ressquared,eflag,outputu] = ...fminunc(Fsumsquares,x0,opts)fprintf(['There were %d iterations using fminunc,' ...' and %d using lsqcurvefit.\n'], ...outputu.iterations,output.iterations)fprintf(['There were %d function evaluations using fminunc,' ...' and %d using lsqcurvefit.'], ...outputu.funcCount,output.funcCount)type fitvectorx02 = [1 0];F2 = @(x,t) fitvector(x,t,y);[x2,resnorm2,~,exitflag2,output2] = lsqcurvefit(F2,x02,t,y)fprintf(['There were %d function evaluations using the 2-d ' ...'formulation, and %d using the 4-d formulation.'], ... output2.funcCount,output.funcCount) x0bad = [5 1 1 0];[xbad,resnormbad,~,exitflagbad,outputbad] = ...lsqcurvefit(F,x0bad,t,y)hold onplot(t,F(xbad,t),'g')legend('Data','Global fit','Bad local fit','Location','NE')hold offfpri ntf(['The residual norm at the good ending point is %f,' ...'and the residual norm at the bad ending point is %f.'], ...res no rm,res no rmbad)displayE ndOfDemoMessage(mfile name)拟合效果如下：直线的最小二乘拟合：y= a+bx式中有两个待定参数，a代表截距，b代表斜率。
