线性拟合法讲解

摘要
摘要
插值法和曲线拟合是两种来源于实际，同时又广泛应用于实际的重要的数值计算方法。

随着计算机技术的不断发展以及人类计算机水平的逐步提高，他们在国民经济和科学研究中占据了越来越重要的地位。

插值法与曲线拟合结合计算机技术例如MATLAB等编译工具可以用来解决许多实际问题，可以做到高效快捷准确的计算出想要的结果。

本文从MATLAB的功能特点出发，研究了数值计算方法中的插值法和曲线拟合两类问题，比较了这两类问题的特点和不同之处，通过多组实验来进行进一步的研究，即使用MATLAB实现通过拉格朗日插值法和曲线拟合解决实际问题。

本文中实现了通过拉格朗日插值法解决三个实际问题，包括了二氧化硫与传感器电压问题，最近十年93#汽油价格变化问题以及中石油股票月K 线并作图分析。

同时还实现了通过多项式曲线拟合解决两个实际问题，即根据研究氮肥（N）的施肥量与土豆产量的影响所得数据曲线拟合出相关函数关系和铁在冶炼过程中含碳量与时间的关系分析并做图研究得到相应的结论。

关键词：插值法曲线拟合 MATLAB 数值计算
ABSTRACT
ABSTRACT
Interpolation and curve fitting are two methods come from actual , while widely used in actual important numerical . With the continuous development of computer technology and the gradual improvement of human computer skills , they occupy an increasingly important position in the national economy and scientific research . Interpolation and curve-fitting combining with computer technology and any other building tools such as MATLAB can be used to solve many practical problems by which can be solved quickly and efficiently accurately calculate the desired results. In this paper, starting from the function features of MATLAB to study the numerical method of interpolation and curve-fitting problems , comparing the characteristics of these two types of problems and differences , to carry out further research by two experiments that use MATLAB implementation to solve practical problems by Lagrange interpolation method to achieve a number of practical problems solved by Lagrange interpolation method , including sulfur dioxide and sensor voltage problem , the last decade of 93 # gasoline prices and oil stocks change on K line soldiers plot analysis .It also achieved by polynomial curve fitting solve two practical problems , according to a study of nitrogen (N) fertilization resulting impact on the amount of data and the potato yield curve fitting and the correlation function of the carbon content of iron in the smelting process and time relationship analysis and make the corresponding figure study conclusions .
Keywords: Interpolation Curve-Fitting MATLAB Numerical calculation
目录1
目录
第一章绪论 (1)
1.1插值法概述 (1)
1.1.1 插值法的背景 (1)
1.1.2 插值法的思想 (2)
1.2.曲线拟合概述 (3)
1.2.1曲线拟合的背景 (3)
1.2.2曲线拟合的思想 (3)
1.3本文研究内容 (4)
第二章插值法基本理论 (5)
2.1插值法基本定义 (5)
2.1.1唯一定理 (5)
2.1.2几何意义 (5)
2.2 拉格朗日（Lagrange）插值 (6)
2.3牛顿（Newton）插值 (7)
2.4样条函数插值方法 (8)
2.4.1二次样条函数插值 (10)
2.4.2三次样条函数插值 (11)
2.5高次插值的龙格现象 (13)
2.6插值法小结 (14)
第三章曲线拟合基本理论 (15)
3.1最小二乘拟合法 (15)
3.2最小二乘拟合函数的求解 (16)
第四章基于MATLAB的插值法仿真研究 (17)
4.1算法及流程图 (17)
4.2代码解析 (17)
4.3插值法应用 (19)
4.3.1二氧化硫浓度与电压 (19)
2 基于Matlab的插值法与曲线拟合数值计算方法研究
4.3.2 93#汽油价格计算 (21)
4.3.3中国石油股票月K线分析 (23)
4.4插值法应用小结 (24)
第五章基于MATLAB的曲线拟合仿真研究 (25)
5.1基于MATLAB实现曲线拟合 (25)
5.2曲线拟合应用 (25)
5.2.1氮肥的施肥量与土豆产量之间的关系分析 (25)
5.2.2铁的冶炼过程中含碳量和时间关系分析 (27)
5.3曲线拟合应用小结 (29)
第六章总结 (31)
6.1插值法与曲线拟合比较 (31)
6.2全文总结 (31)
致谢 .............................................. 错误！未定义书签。

参考文献 . (35)
第一章绪论1
第一章绪论
在工程实践和科学实验中，经常会遇到计算函数问题，就是数值方法本身而言有很多问题最后也都是转化成了函数数值的计算。

然而，有时候遇到非常复杂的函数关系，甚至没有明确的表达式或者有表达式但却无法知道。

插值与拟合往往就是解决这些问题的方法。

例如一些问题需要建立函数关系，即y=f(x)。

虽然从原则上说，它在某个区间[a,b]上是存在的，但通常只能观测到它的部分信息，即只能获取[a,b]上一系列离散点上的值，这些值构成了观测数据。

这就是说，我们只知道的一张关于x与f(x)的观测数据表，如下：
表1.1 关于x与f(x)的观测数据表
xi x1 x2 …xn f(xi) f(x1) f(x2) …f(xn) 然而，当不知道函数在其他点x上的取值时，只能用一个经验函数y=g(x)对真实函数y=f(x)作近似。

下面本文将引入两种两种常用办法来确定经验函数y=g(x)，它们分别是插值法与曲线拟合法。

1.1插值法概述
1.1.1 插值法的背景
插值法是一种古老的数学方法，早在一千多年前的隋唐时期定制历法时就广泛应用了二次插值.二次插值法的创立, 是隋唐数学的一项重大成绩。

插值法是根据两个自变量的已知函数值的近似计算方法。

这种方法是很有使用价值的。

例如, 在天文观测中, 人们不可能每时每刻都进行观测, 因此只能得到日月五星某些时刻在天球上的位置。

利用这些观测记录推算日月五星在其他时刻的位置, 就要用到插值法, 这对于天文计算特别是日月交食的推算是十分重要的。

实际上在《周髀》和《九章》中就已有了一次插值公式。

东汉末天文学家刘洪制定《乾象历》, 为计算月球在近地点后( n+s) 日的共行度数, 采用了一次插值公式
[1]:f(n+s)=f(n)+s△,0<s<1, △=f(n+1)- f(n)其中n 为月球在近地点后运行的整日数, f( n) 为对应的月球位置函数。

此后, 曹魏杨伪、姚秦姜岌、刘宋何承天、南齐祖
2 基于Matlab的插值法与曲线拟合数值计算方法研究
冲之等各家历法计算月行数时也采用了这种算法。

随着天文学的发展和观测精度的提高, 天文学家不仅发现了月球运动的不均匀性, 而且也发现了太阳和五星视运动的不均匀性, 也就是说, 日月五星的视运动并非是时间的一次函数。

为了编制更好的历法, 特别是为了精确计算合朔( 太阳与月亮同时出现) 和日食时刻, 何承天、祖冲之以前所长期采用的一次插值法, 误差太大, 已经不能满足这种要求, 于是中国天算家开始了新的探索。

然而插值理论却是在17世纪微积分产生后才逐步发展起来，Newton插值公式理论是当时的重要成果[1]。

由于计算机的使用以及航空、造船、精密仪器的加工，插值法在理论和实践上都得到进一步发展，获得了广泛的应用。

在计算机广泛使用的今天, 插值法跟MATLAB 等软件的程序结合起来, 发挥其更广泛的作用, 以至可以解决高达几百次甚至上千的插值拟和( 如果计算机允许的话) 。

为此, 大量的数学家和软件编辑者投入到数值分析和数学软件的工作中, 使得插值法得到了空前的发展, 现在的Lagrange 插值, 逐次线性插值, Newton 插值, Hermite 插值,分段低次插值, 还有三次样条插值等等, 就不胜枚举。

1.1.2 插值法的思想
已知如下一个关于xi与f(xi)的观测数据表
表1.2 插值法思想观测数据表
xi x1 x2 …xn
f(xi) f(x1) f(x2) …f(xn) 求一个经验函数y=g(x)，使
.
n
…
…
,1
),
(
)
(=
=i
x
f
x
g
i
i式（1-1）
插值的任务就是由已知的观测点(xi,yi)为物理量(未知量),建立一个简单的、连续的解析模型g(x) ，以便能根据该模型推测该物理量在非观测点处的特性[2]。

第一章 绪论
3
图1.1 插值法思想 1.2.曲线拟合概述
1.2.1曲线拟合的背景
实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；
疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。

曲线拟合（curve fitting ）是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系。

1.2.2曲线拟合的思想
曲线拟合的基本思想是根据已知数据表，如下表：
表1.3 曲线拟合数据观测表 xi
x1 x2 … xn f(xi) f(x1) f(x2) … f(xn)
求一个经验函数y=g(x)，使得
min ))()((21=-∑=n i i
i x f x g 式（1-2）
所求的的函数g(x)的拟合函数。

而每个数据点分布在拟合曲线的两侧，拟合曲线不需要经过数据点。

如下图：
o x y ● ●
●
● ● y 0 x 1
x 2 x n y 1
y 2
y n x 0 y=f(x) g(x)
4 基于Matlab的插值法与曲线拟合数值计算方法研究
图1.2 曲线拟合思想
1.3本文研究内容
本文题目是基于MATLAB的插值法与曲线拟合的数值计算方法研究。

各章节内容分别下：
第一章，绪论，即研究学习了解插值法与曲线拟合相关研究背景以及数学思想。

第二章，插值法基本理论，即深入的学习研究了插值法，包括拉格朗日插值，牛顿差值，样条插值等多种插值方法。

第三章，曲线拟合基本理论，即深入学习研究了曲线拟合，最主要的最小二乘法。

以及最小二乘方程的解法。

第四章，基于MATLAB的插值法仿真研究，即使用MATLAB实现了通过拉格朗日插值法解决多个实际问题，包括了二氧化硫与传感器电压问题，最近十年93#
汽油价格变化问题以及中石油股票月K线并作图分析。

第五章，基于MATLAB的曲线拟合仿真研究，即使用MATLAB实现了通过多项式曲线拟合解决两个实际问题。

根据研究氮肥（N）的施肥量与土豆产量的影响所得数据曲线拟合出相关函数关系和铁在冶炼过程中含碳量与时间的关系分析并做图分析。

第六章，总结，即比较插值法与曲线拟合两类问题的特点与不同之处，总结全文研究要点。

第二章 插值法基本理论
5
第二章 插值法基本理论
2.1插值法基本定义
设函数)(x f y =在一个区间],[b a 上连续，且在n+1个不同的点
b x x x a n ≤≤,,,10 上分别取值n y y y ,,,10 ，
在一个性质优良且便于计算的函数类φ 中，求一个简单函数)(x p ，使得
()()n i y x P i i ,1,0== 式（2-1）
而在i x x ≠上，作为)(x f 的近似。

称该区间为插值区间，点 为插值节点，称式（2-1）为 )(x f 的插值条件，并称函数类φ为插值函数类，称 )(x p 为函数在节点 n x x x ,,,10 处的插值函数。

求插值函数 p(x) 的方法就称为插值法。

插值函数类φ的取法不同，所求得的插值函数p(x)逼近f(x)的效果就不同它的选择取决于使用上的需要。

常用的方法有代数多项式、三角多项式和有理函数等[3]。

当选用代数多项式作为插值函数时，那相应的插值问题就称为多项式插值。

在多项式插值中，最常见、最基本的问题是：求一次数不超过n 的代数多项式
()n n n x a x a a x P +++= 10 式（2-2）
使 ()()n i y x P i i n ,,,1,0 == 式（2-3）
其中 n a a a ,,,10 为实数。

满足插值条件(2-3)的多项式(2-2)，称为函数f(x) 在节点处的n 次插值值多项式。

2.1.1唯一定理
满足插值条件的、次数不超过n 的多项式是存在而且是唯一的。

2.1.2几何意义
n 次插值多项式 ()x P n 的几何意义：过曲线y = f(x) 上的n+1个点
6 基于Matlab 的插值法与曲线拟合数值计算方法研究
),,1,0)(,(n i y x i i =作一条n 次代数曲线
()x P n ，作为曲线y = f(x) 的近似，如图
2.1。

图2.1 插值法几何意义
2.2 拉格朗日（Lagrange ）插值
设函数)(x f y =在1+n 个相异点n x x x x ,,,,210⋅⋅⋅上的函数值为n y y y y ,,,,210⋅⋅⋅，要求一个次数不得超过n 的代数多项式
n n n x a x a x a a x P +⋅⋅⋅+++=2210)( 式（2-4）
使得在节点i x 上有),,2,1,0()(n i y x P i i n ⋅⋅⋅==成立，称为n 次代数插值问题，)(x P n 称为插值多项式．则可以证明n 次代数插值是唯一的。

事实上: 可以得到：
j n j n i i j i n y x x x x x P j i ∑∏==⎥⎥⎦
⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛--=≠00)()( 式（2-5） 当1=n 时，有二点一次（线性）插值多项式：
10
1001011)(y x x x x y x x x x x P --+--= 式（2-6） 当n =2时，有三点二次（抛物线）插值多项式：
()x P y n =()x f y =0x 1x n x X a b 0y 1y n
y 0Y
第二章 插值法基本理论
7
2120210121012002010212)
)(())(())(())(())(()
)(()(y x x x x x x x x y x x x x x x x x y x x x x x x x x x P ----+----+----=
式（2-7）
2.3牛顿（Newton ）插值
牛顿插值的基本思想：
由于)(x f y =关于二节点10,x x 的线性插值为：
)()
()()()()()()()(00
101000010101x x x x x f x f x p x x x x x f x f x f x p ---+=---+
=式（2-8）
假设满足插值条件)2,1,0()()(2===i x p y x f i i i 的二次插值多项式的一般形式为
))(()()(1020102x x x x c x x c c x p --+-+= 式（2-9）
由插值条件可得：
⎪⎩⎪
⎨⎧=--+-+=-+=)
())(()()()()(2120220210
1011000x f x x x x c x x c c x f x x c c x f c 式（2-10）
可以解出，得：
⎪⎪⎪⎩
⎪
⎪⎪⎨
⎧------=--==020
101121220
101100
)()()()()()(),(x x x x x f x f x x x f x f c x x x f x f c x f c 式（2-11） 所以：
))(()())(()()(10211020102x x x x c x p x x x x c x x c c x p --+=--+-+=式（2-12） 类似的方法，可以得到三次插值多项式等，如果按这种思想则可以得到一般的牛顿插值公式[4]。

函数的差商及其性质：
对于给定的函数)(x f ，用),,,(10n x x x f ⋅⋅⋅表示关于节点n x x x ,,,10⋅⋅⋅的n 阶差商，则有：
8 基于Matlab 的插值法与曲线拟合数值计算方法研究
一阶差商：010110)()(),(x x x f x f x x f --=
，1
21221)
()(),(x x x f x f x x f --=
二阶差商：0
21021210)
,(),(),,(x x x x f x x f x x x f --=
由此得：n 阶差商：0
1102110)
,,,(),,,(),,,(x x x x x f x x x f x x x f n n n n -⋅⋅⋅-⋅⋅⋅=⋅⋅⋅-
差商有以下性质：
（1）差商的分加性：∑
∏=≠=-=⋅⋅⋅n
k n
k j j j k
k n x x
x f x x x f 0)
(010)
()
(),,,(。

（2）差商的对称性：在),,,(10n x x x f ⋅⋅⋅中任意调换j i x x ,的次序其值不会变。

牛顿插值法公式： 一次插值公式为：
))(,()()(01001x x x x f x f x p -+= 式（2-13）
二次插值公式为：
)
)()(,,()()
)()(,,())(,()()(1021011021001002x x x x x x x f x p x x x x x x x f x x x x f x f x p --+=--+-+=式（2-14）
于是有一般的牛顿插值公式为：
)
())()(,,,()()())()(,,,())()(,,())(,()()(11010111010102100100----⋅⋅⋅--⋅⋅⋅+=-⋅⋅⋅--⋅⋅⋅+⋅
⋅⋅+--+-+=n n n n n n x x x x x x x x x f x p x x x x x x x x x f x x x x x x x f x x x x f x f x p 式（2-15）
可以证明：其余项为：
))(())()(,,,,()(11010n n n x x x x x x x x x x x x f x R --⋅⋅⋅--⋅⋅⋅=- 式（2-16）
实际上，牛顿插值公式是拉格朗日插值公式的一种变形，二者是等价的。

另外还有著名的埃尔米特（Hermite ）插值等。

2.4样条函数插值方法
样条（spline ），在英语中是指富有弹性的细长木条。

样条函数实质上是由分段多项式光滑连接成的函数，一般称之为多项式样条。

由于样条函数的特殊性质，
第二章 插值法基本理论
9
决定了样条函数在实际中有着很多重要的应用。

样条函数的一般概念：
定义：设给定区间],[b a 的一个分划b x x x a n =<⋅⋅⋅<<=∆10:，如果函数)(x s 满足条件：(1) 在每个子区间),,2,1](,[1n i x x i i ⋅⋅⋅=-上是k 次多项式； (2) )(x s 及直到k -1阶的导数在],[b a 上连续。

则称)(x s 是关于分划△的一个k 次多项式样条函数，n x x x ,,,10⋅⋅⋅称为样条节点，
121,,,-⋅⋅⋅n x x x 称为内节点，n x x ,0称为边界节点，这类样条函数的全体记作),(k S P ∆，称为k 次样条函数空间。

若),()(k S x s P ∆∈，则)(x s 是关于分划△的k 次多项式样条函数。

k 次多项式样条函数的一般形式为
∑
∑
=-=+-+=k
i n j k j j
i
i k x x k i x x s 0
1
1
)(!
!
)(βα 式（2-17）
其中),,1,0(k i i =α和)1,,2,1(-=n j j β均为任意常数，而
)1,,2,1(,0,)()(-=⎪⎩⎪⎨⎧<≥-=-+n j x x x x x x x x j
j k
j k
j 式（2-18）
在实际中最常用的是2=k 和３的情况，即为二次样条函数和三次样条函数。

二次样条函数：对于],[b a 上的分划b x x x a n =<⋅⋅⋅<<=∆10:,则
)2,()(!
2!
2)(1
1
22
2
102∆βαααP n j j j
S x x x x x s ∈-++
+=∑
-=+式（2-19）
其中)1,2,1(,0,)()(2
2-=⎪⎩⎪⎨
⎧<≥-=-+
n j x x x x x x x x j
j
j j 。

三次样条函数：对于],[b a 上的分划b x x x a n =<⋅⋅⋅<<=∆10:,则
)3,()(!
3!
3!
2)(1
1
33
3
2
2
103∆βααααP n j j j
S x x x x x x s ∈-++
+
+=∑
-=+式（2-20）
其中)1,2,1(,0,)()(3
3
-=⎪⎩⎪⎨
⎧<≥-=-+
n j x x x x x x x x j
j
j j 。

10 基于Matlab 的插值法与曲线拟合数值计算方法研究
2.4.1二次样条函数插值
)2,()(2∆∈P S x s 中含有2+n 个待定常数，故应需要2+n 个插值条件，因此，二次样条插值问题可分为两类：
问题(1)：已知插值节点i x 和相应的函数值),,2,1,0(n i y i ⋅⋅⋅=，以及端点0x (或n x )处的导数值0'y (或n y ')，求)2,()(2∆∈P S x s 使得
⎩⎨⎧'=''='⋅⋅⋅==))(()(),,2,1,0()(20022n
n i i y x s y x s n i y x s 或 式（2-21） 问题(2)：已知插值节点i x 和相应的导数值),,2,1,0(n i y i ⋅⋅⋅='，以及端点0x (或n x )处的函数值0y (或n y )，求)2,()(2∆∈P S x s 使得
⎩⎨⎧==⋅⋅⋅='='))(()(),,2,1,0()(200
22n n i i y x s y x s n i y x s 或 式（2-22） 事实上，可以证明这两类插值问题都是唯一可解的． 对于问题(1)，由条件(2-21)
⎪⎪⎪⎪⎩⎪⎪⎪
⎪⎨⎧'=+='==-+++==++==++=∑-=0
0210211
2
2210212121111
202
020100
2)(,,3,2,)(2121)(21)(2
1)(y x x s n
j y x x x x x s y
x x x s y x x x s j j i i j i j j j ααβααααααααα 引入记号T n ),,,,,(11210-=ββααα X 为未知向量，T n
n y y y y ),,,,(10'= C 为已知向量，
⎥⎥⎥
⎥⎥⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡---=-00
1
)(21)(2
1
211)(2
1
211
2112110
2
12122122222112
00
x x x x x x x x x x x x x x x n n n n n A 于是，问题转化为求方程组C AX =的解T n ),,,,,(11210-=ββααα X 的问题，
第二章 插值法基本理论
11
即可得到二次样条函数的)(2x s 的表达式。

2.4.2三次样条函数插值
由于)3,()(3∆∈P S x s 中含有3+n 个待定系数，故应需要3+n 个插值条件，因此可将三次样条插值问题分为三类：
问题(1)：已知插值节点j x 和相应的函数值),,2,1,0(n j y j ⋅⋅⋅=，以及两个端点0x ，n x 处的导数值0'y ，n y '，求)3,()(3∆∈P S x s 使满足条件
⎪⎩⎪⎨
⎧='='⋅⋅⋅==),0()()
,,1,0()(33n j y x s n j y x s j j
j j 式（2-23） 问题(2)：已知插值节点
j
x 和相应函数值
)
,,2,1,0(n j y j ⋅⋅⋅=，以及两端点0x ，n x 处
的二阶导数值0
y ''，n y ''，求)3,()(3∆∈P S x s 使满足条件 ⎪⎩⎪⎨
⎧=''=''⋅⋅⋅==),0()()
,,1,0()(33n j y x s n j y x s j j
j j 式（2-24） 问题(3)：类似地，求)3,()(3∆∈P S x s 使满足条件
⎪⎩⎪⎨⎧=+=-==)
2,1,0)(0()0(),,1,0()(0)
(3)(33k x s x s n j y x s k n k j j 式（2-25） 这三类插值问题的条件都是3+n 个，可以证明其解都是唯一的[5]。

一般的求解方法可以仿照二次样条的情况处理方法，在这里给出一种更简单的方法．仅依问题(1)为例，问题(2)和问题(3)的情况类似处理。

由于在)3,()(3∆P S x s ∈区间],[b a 上是一个分段光滑，且具有二阶连续导数的三
次多项式，则在子区间],[1+j j x x 上)(3
x s ''是线性函数，记),,,1,0)((3n j x s d j j =''=为待定常数．由拉格朗日插值公式可得
n j x x h h x x d h x x d x s j j j j
j j j
j j ,,1,0,,)(11
13
=-=-+-=''+++
显然j
j
j h d d x s -='''+13)(在],[1+j j x x 上为常数．于是在],[1+j j x x 上有
31233)(6)(2
))(()(j j
j
j j j j j j x x h d d x x d x x x s y x s --+
-+-'+=+ 式（2-26）
12 基于Matlab 的插值法与曲线拟合数值计算方法研究
则当1+=j x x 时，由式(2-26)和问题(1)的条件得
12123
136
2)()(+++=-+
+'+=j j j
j j j j j j j y h d d h d h x s y x s
故可解得
)2(6
)(113
+++--='j j j j
j
j j d d h h y y x s 式（2-27）
将式(2-27)代入式(2-26)得
)
1,,1,0](,[,)(6)(2)()2(6)(1312
113-=∈--+-+-⎥⎥⎦
⎤⎢⎢⎣⎡+--+=++++n j x x x x x h d d x x d x x d d h h y y y x s j j j j
j
j j j
j j j j j j j j 式（2-28）
在],[1j j x x -上同样的有
),,2,1](,[,)(6)(2)()2(6)(1311
12
11
1111113n j x x x x x h d d x x d x x d d h h y y y x s j j j j j j j j j j j j j j j j =∈--+-+-⎥⎥⎦⎤⎢⎢⎣⎡+--+=------------ 式（2-29）
根据)(3x s 的一阶导数连续性，由式(2-29)得
)()2(6
)0(3
111
1
3
j j j j j j j j x s d d h h y y x s '=++-=-'---- 结合式(2-27)整理得
⎪⎪⎭
⎫ ⎝⎛--
-+=++
++--+-+----11111
1
11
162j j j j j j j j j j j j j j j j j h y y h y y h h d h h h d d h h h 引入记号⎪⎪⎭
⎫ ⎝⎛--
-+=
+=
--+--11111
6,j j j j j j j j j j j j j h y y h y y h h c h h h a ，111--+=-j j j j h h h a ．则 )1,,2,1(,2)1(11-==++-+-n j c d a d d a j j j j j j 式（2-30）
再由边界条件：n n y x s y x s '=''=')(,)(3003
得
第二章 插值法基本理论
13
⎪⎪⎩
⎪⎪⎨
⎧⎪⎪⎭⎫ ⎝⎛--'=+⎪⎪⎭
⎫ ⎝⎛'--=+----1111
000
10106262n n n n n n n h y y y h d d y h y y h d d 式（2-31） 联立式(2-30),(2-31)得方程组
C D A =⋅
其中
⎥⎥⎥⎥⎥⎥⎥
⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡---=----2121212112112
200n n n n a a a a a a
A ，⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=-n n d d d d 110 D ，⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣
⎡⎪⎪⎭⎫ ⎝
⎛--'⎪⎪⎭⎫ ⎝⎛'--=----1
11
110001066n n n n n n h
y y y h c c y h y y h C 式（2-32） 由方程组(2.32)可以唯一解出),,1,0(n j d j =，代入式（2-28）就可以得三次样条函数)(3x s 的表达式．
除此以外还有很多其它的插值方法例如分段线性插值、二次样条插值、三次样条插值等。

2.5高次插值的龙格现象
插值多项式余项公式说明插值节点越多，一般说来误差越小，函数逼近越好，但这也不是绝对的，因为余项的大小既与插值节点的个数有关，也与函数f(x)的高阶导数有关。

换句话说，适当地提高插值多项式的次数，有可能提高计算结果的准确程度，但并非插值多项式的次数越高越好。

当插值节点增多时，不能保证非节点处的插值精度得到改善，有时反而误差更大。

考察函数 55,11)(2
≤≤-+=
x x
x f 右图给出了)(5x P 和)
(10x P 的图像,当n 增大时,)(x P n 在两端会发出激烈的振荡,这就是所谓龙格现象。

该现象表明,在大范围内使用高次插值,逼近的效果往往是不理想的 。

14 基于Matlab的插值法与曲线拟合数值计算方法研究
y
P10(x)
f(x)
P5(x)
-5 0 5 x
图2.2 龙格现象
2.6插值法小结
插值法是实用性很强的方法，它们解决的实际问题虽然各式各样，但抽象为数学问题却有它的共性，即利用已知的数据去寻求某个较为简单的函数P(x)来逼近f(x)。

插值法给出了寻求这种近似函数的原则，以及构造近似函数的几种具体方法。

插值法要求近似函数在已知的数据点必须与f(x)完全一致。

第三章 曲线拟合基本理论
15
第三章 曲线拟合基本理论
插值方法要求插值曲线严格地通过所给的每一个数据点（样点），数据本身存在的不可避免的误差会反映到插值结果中。

如果数据误差很大的话，用插值方法显然就不恰当了，可以用曲线拟合方法。

曲线拟合法: 设法构造一条曲线反映所给数据点的总趋势(但不一定过点），以消除其局部的波动 。

3.1最小二乘拟合法
如果已知函数f(x)在若干点xi(i=1,2,…,n)处的值yi,便可根据插值原理来建立插值多项式作为f(x)的近似。

但是在科学实验和生产实践中，往往会遇到这样一种情况，即节点上的函数值并不是很精确，这些函数值很多是由实验或观测得到的数据，不可避免地带有测量误差，如果要求所得的近似函数曲线精确无误地通过所有的点(xi,yi),就会使曲线保留所有这一切测试误差。

那么当个别数据的误差较大时,插值效果显然是不会理想的。

此外,由实验或观测提供的数据个数往往很多,如果用插值法,势必得到次数较高的插值多项式，这样计算起来会很烦琐。

换句话说:求一条曲线,使数据点均在离此曲线的上方或下方不远处,所求的曲线称为拟合曲线,它既能反映数据的总体分布,又不至于出现局部较大的波动,更能反映被逼近函数的特性,使求得的逼近函数与已知函数从总体上来说其偏差按某种方法度量达到最小,这就是曲线拟合[6]。

与函数插值问题不同,曲线拟合不要求曲线通过所有已知点,而是要求得到的近似函数能反映数据的基本关系。

在某种意义上,曲线拟合更有实用价值。

其中令
i i i y x y -=)(δ在回归分析中称为残差。

一般使用∑==m i i 0
2
22δδ
∑=-=m
i i i y x y 0
2))((作为
衡量)(x y 跟数据点),(i i y x 偏离程度大小的度量标准，称为平方误差。

在对给出的实验(或观测)数据作曲线拟合时,怎样才算拟合得最好呢？一般希望各实验(或观测)数据与拟合曲线的偏差的平方和最小,这就是最小二乘原理。

16 基于Matlab 的插值法与曲线拟合数值计算方法研究
3.2最小二乘拟合函数的求解
要使最小二乘问题的目标函数(5.17)达到最小，则由多元函数取得极值的必要条件得
),,2,1,0(0m k a S
k
==∂∂ 式（3-1） 即
),,2,1,0(0)()(10m k x y x a w i k n
i i m k i k k i ⋅⋅⋅⋅==⎥⎦
⎤⎢⎣⎡-∑∑==φφ 式（3-2） 亦即
),,2,1,0()()()(001m k x y w a x x w n
i i k i i j m
j n i i k i j i ⋅⋅⋅⋅==⎥⎦
⎤
⎢⎣⎡∑∑∑===φφφ 式（3-3）
是未知量为m a a a a ,,,,210⋅⋅⋅的线性方程组，称(3-3)式为正规方程组。

实际中可适当选择函数系{}m
k k x 0)(=φ，由正规方程组解出m a a a a ,,,,210⋅⋅⋅，于是
可得最小二乘拟合函数∑==m
k k k x a x 0
)()(φφ。

第四章 基于MA TLAB 的插值法仿真研究
17
第四章 基于MATLAB 的插值法仿真研究
4.1算法及流程图
我们可以从前文得知，拉格朗日插值多项式为：
j n
j n i i j i
n y x x x
x x P j i ∑∏==⎥⎥⎦
⎤
⎢⎢⎣⎡⎪⎪⎭⎫
⎝⎛--=≠00)
()( 式(4-1) MATLAB 中没有提供现成的函数用于实现拉格朗日插值，我们可以自己编写lang.m 函数来实现拉格朗日插值。

如图4.1，我们画出lang 函数的实现流程图。

从图中我们可以看到，函数刚开始需要校验x ，y 两向量长度是否相同，然后求出拉格朗日基函数，最后经过循环求出插值多项式并求出x0点处对应的函数值。

4.2代码解析
结合以上流程图，我们写出了以下MATLAB 程序，实现了拉格朗日插值法[7]。

通过多次for 循环将基函数运算，求出拉格朗日插值多项式。

function s=Lagrangenew(x,y,x0)
%lagrange 插值，x ，y 为已知的插值点及其函数值 %x0为要求的插值点的x 值 nx=length(x); ny=length(y); if nx~=ny
warning('矢量x 与y 的长度应该相等') return end
m=length(x0);
%按照公式，对要求的插值点矢量x0的每个元素进行计算 for i=1:m
t=0.0;
18 基于Matlab的插值法与曲线拟合数值计算方法研究
for j=1:nx
u=1.0;
for k=1:nx
if k~=j
u=u*(x0(i)-x(k))/(x(j)-x(k));
end
end
t=t+u*y(j);
end
s(m)=t;
end
Return
开始
输入插值点（xi，yi）
输入要求的值x0
初始result=0
求出x向量长度n(x)
求出y向量长度n(y)
判断n(x),n(y)是否相等Error
I=1
临时变量temp=yi
i=j?
Temp=temp*(x-xi) Temp=temp*(xi-xj) Result=result=temp
i=n?
输出结果result
结束J=j+1
J<n
i=i+1
否
是
=
≠
≠
=
图4.1 插值法流程图
第四章基于MA TLAB的插值法仿真研究19
4.3插值法应用
在研究过程，我做了很多的尝试，尝试用插值法解决实际问题，用以体现插值法的高效准确性，本文中将提到我的三次尝试与试验。

首先是二氧化硫传感器中二氧化硫浓度与传感电压之间的关系，其次是最近十年来93#汽油价格变化问题，最后就是中国石油股票（601857）2013年各月月末收盘价变化。

4.3.1二氧化硫浓度与电压
首先，我们通过实验室获取了一组2SH12二氧化硫传感器数据，我们打算根据这十组数据来观察插值法是否适用于这些采集自实验室的数据。

如下：
表4.1 2SH12二氧化硫传感器数据V/PPM,R0=2k，电压信号单位为V.
编号/浓度0ppm 50 100 150 200 250 300 350 400 450 500
VI70601135 0.276 1.46 1.67 1.82 1.97 2.08 2.2 2.3 2.4 2.47 2.53
当我们代入五组以上数据到我们通过MATLAB实现的拉格朗日插值法中时，我们看一下结果：
当浓度为250ppm时，电压值为2.025V，如图：
图4.2 250ppm时5点插值法
20 基于Matlab的插值法与曲线拟合数值计算方法研究
当浓度为400ppm时，传感器电压值为2.588.如图：
图4.3 400ppm时5点插值法
然后，我们增加数据带入个数，即增加插值点个数。

同时，我们将浓度值为100ppm时的电压值设置为未知，根据插值法计算当时的电压值得v=1.5856。

程序运行结果如下：
图4.4 100ppm时10点插值
当把浓度值为250ppm时，计算得传感器电压值为2.0951V。

第四章基于MA TLAB的插值法仿真研究21
图4.5 250ppm浓度10点插值
继续实验，当浓度值为400ppm时，传感器电压值为2.3156V.
图4.6 400ppm浓度时10点插值
根据以上五次试验，我们可以发现，插值法求出的电压值基本符合实际变化趋势，且误差较小。

当二氧化硫传感器的浓度为已知数值时，我们就可以大致求出该浓度下传感器的电压值。

同时我们可以看出，随着插值点的增加，误差会逐步减小。

4.3.2 93#汽油价格计算
查阅中国网上历年来的资料，我们了解到了2004年到2013年历年来我国93#汽油的平均价格[8]，如下，我们想要看看这些来自与市场的数据是否可以使用插值法来求解。

22 基于Matlab 的插值法与曲线拟合数值计算方法研究
表4.2 93#汽油历年油价数据
我们将上述数据带入到我们我们通过MATLAB 实现的拉格朗日插值法程序当中，以便经行我们的实验。

当把2009年的油价作为未知数据时，我们可以通过程序计算得到此时油价为6.321，同时得到运行结果：
图4.7 2009年油价插值
随后我们进行了多组实验，得到如下实验结果：
表4.3 历年油价真实值与计算值比较
年份 2005 2006 2007 2008 2009 2010 2011 实际油价/元 3.987 4.779 4.87 5.66 5.754 6.69 7.50 插值油价/元
11.92
2.795
5.72
5.09
6.321
5.839
9.789
我们可以从以上实验数据中明显发现，2004到2013年以来93#油价的计算值与实际值之间有很大的误差，经过一系列分析查证，我们了解到，插值法所描绘的图像是经过每一点的折线，所以每个点的数据误差都会对对计算造成很大影响，而本次实验的数据采集本身具有局限性同时数据变动也较大，所以插值法不适于用来解决波动较大，误差较大的实际问题。

年份 2004 2005 2006 2007 2008 2009
2010 2011 2012
2013
93#油价/元 3.56 3.987 4.779 4.87 5.66
5.754
6.69
7.50
7.616 7.865
第四章基于MA TLAB的插值法仿真研究23 4.3.3中国石油股票月K线分析
我们在东方财经网上找到了中石油股票（601857）2013年月K线，并读取了各月月末收盘价数据，如下：
图4.8 中国石油2013年月k线
表4.4 中石油股票2013个月月末收盘价
月份 1 2 3 4 5 6 7 8 9 10 11 12
收盘价9.26 9.05 8.69 8.48 8.52 8.35 8.08 7.83 7.84 7.84 7.91 7.71 我们想要看看拉格朗日插值法是否可以计算出2013年某个特定月份的月末收盘价，前提是已知了当月其他月份的相关数据。

于是我们进行了如下实验：
将以上数据带入到程序中，当已知其他月份求8月份收盘价时，我们得出结果：当x=8时，收盘价y=7.903.如下图：
图4.9 8月份油价插值。