实验二_插值与拟合
- 格式:doc
- 大小:344.92 KB
- 文档页数:14
数值计算方法插值与拟合数值计算方法在科学计算和工程应用中起着重要的作用,其中插值和拟合是其中两个常用的技术。
插值是指通过已知的离散数据点来构造出连续函数或曲线的过程,拟合则是找到逼近已知数据的函数或曲线。
本文将介绍插值和拟合的基本概念和常见的方法。
一、插值和拟合的基本概念插值和拟合都是通过已知数据点来近似表达未知数据的方法,主要区别在于插值要求通过已知数据点的函数必须经过这些数据点,而拟合则只要求逼近这些数据点。
插值更加精确,但是可能会导致过度拟合;拟合则更加灵活,能够通过调整参数来平衡拟合精度和模型复杂度。
二、插值方法1. 线性插值线性插值是一种简单的插值方法,通过已知数据点构造出线段,然后根据插值点在线段上进行线性插值得到插值结果。
2. 拉格朗日插值拉格朗日插值是一种基于多项式插值的方法,通过已知数据点构造出一个多项式,并根据插值点求解插值多项式来得到插值结果。
3. 分段线性插值分段线性插值是一种更加灵活的插值方法,通过将插值区间分成若干小段,然后在每个小段上进行线性插值。
三、拟合方法1. 最小二乘法拟合最小二乘法是一种常用的拟合方法,通过最小化实际观测点和拟合函数之间的残差平方和来确定拟合函数的参数。
2. 多项式拟合多项式拟合是一种基于多项式函数的拟合方法,通过选择合适的多项式次数来逼近已知数据点。
3. 曲线拟合曲线拟合是一种更加灵活的方法,通过选择合适的曲线函数来逼近已知数据点,常见的曲线包括指数曲线、对数曲线和正弦曲线等。
四、插值与拟合的应用场景插值和拟合在实际应用中具有广泛的应用场景,比如图像处理中的图像重建、信号处理中的滤波器设计、金融中的风险评估等。
五、插值与拟合的性能评价插值和拟合的性能可以通过多种指标进行评价,常见的评价指标包括均方根误差、相关系数和拟合优度等。
六、总结插值和拟合是数值计算方法中常用的技术,通过已知数据点来近似表达未知数据。
插值通过已知数据点构造出连续函数或曲线,拟合则找到逼近已知数据的函数或曲线。
程序设计实验报告(matlab)实验一: 程序设计基础实验目的:初步掌握机器人编程语言Matlab。
实验内容:运用Matlab进行简单的程序设计。
实验方法:基于Matlab环境下的简单程序设计。
实验结果:成功掌握简单的程序设计和Matlab基本编程语法。
实验二:多项式拟合与插值实验目的:学习多项式拟合和插值的方法,并能进行相关计算。
实验内容:在Matlab环境下进行多项式拟合和插值的计算。
实验方法:结合Matlab的插值工具箱,进行相关的计算。
实验结果:深入理解多项式拟合和插值的实现原理,成功掌握Matlab的插值工具箱。
实验三:最小二乘法实验目的:了解最小二乘法的基本原理和算法,并能够通过Matlab进行计算。
实验内容:利用Matlab进行最小二乘法计算。
实验方法:基于Matlab的线性代数计算库,进行最小二乘法的计算。
实验结果:成功掌握最小二乘法的计算方法,并了解其在实际应用中的作用。
实验六:常微分方程实验目的:了解ODE的基本概念和解法,并通过Matlab进行计算。
实验内容:利用Matlab求解ODE的一阶微分方程组、变系数ODE、高阶ODE等问题。
实验方法:基于Matlab的ODE工具箱,进行ODE求解。
实验结果:深入理解ODE的基本概念和解法,掌握多种ODE求解方法,熟练掌握Matlab的ODE求解工具箱的使用方法。
总结在Matlab环境下进行程序设计实验,使我对Matlab有了更深刻的认识和了解,也使我对计算机科学在实践中的应用有了更加深入的了解。
通过这些实验的学习,我能够灵活应用Matlab进行各种计算和数值分析,同时也能够深入理解相关的数学原理和算法。
这些知识和技能对我未来的学习和工作都将有着重要的帮助。
第1篇一、实验目的1. 理解数值计算的基本概念和常用算法;2. 掌握Python编程语言进行数值计算的基本操作;3. 熟悉科学计算库NumPy和SciPy的使用;4. 分析算法的数值稳定性和误差分析。
二、实验内容1. 实验环境操作系统:Windows 10编程语言:Python 3.8科学计算库:NumPy 1.19.2,SciPy 1.5.02. 实验步骤(1)Python编程基础1)变量与数据类型2)运算符与表达式3)控制流4)函数与模块(2)NumPy库1)数组的创建与操作2)数组运算3)矩阵运算(3)SciPy库1)求解线性方程组2)插值与拟合3)数值积分(4)误差分析1)舍入误差2)截断误差3)数值稳定性三、实验结果与分析1. 实验一:Python编程基础(1)变量与数据类型通过实验,掌握了Python中变量与数据类型的定义方法,包括整数、浮点数、字符串、列表、元组、字典和集合等。
(2)运算符与表达式实验验证了Python中的算术运算、关系运算、逻辑运算等运算符,并学习了如何使用表达式进行计算。
(3)控制流实验学习了if-else、for、while等控制流语句,掌握了条件判断、循环控制等编程技巧。
(4)函数与模块实验介绍了Python中函数的定义、调用、参数传递和返回值,并学习了如何使用模块进行代码复用。
2. 实验二:NumPy库(1)数组的创建与操作通过实验,掌握了NumPy数组的基本操作,包括创建数组、索引、切片、排序等。
(2)数组运算实验验证了NumPy数组在数学运算方面的优势,包括加、减、乘、除、幂运算等。
(3)矩阵运算实验学习了NumPy中矩阵的创建、操作和运算,包括矩阵乘法、求逆、行列式等。
3. 实验三:SciPy库(1)求解线性方程组实验使用了SciPy库中的线性代数模块,通过高斯消元法、LU分解等方法求解线性方程组。
(2)插值与拟合实验使用了SciPy库中的插值和拟合模块,实现了对数据的插值和拟合,并分析了拟合效果。
插值与拟合方法在实际中,常常要处理由实验或测量所得到的一批离散数据.插值与拟合方法就是要通过这些数据去确定某一类已知函数的参数或寻找某个近似函数,使所得到的近似函数与已知数据有较高的拟合精度.插值问题:要求这个近似函数(曲线或曲面)经过所已知的所有数据点.通常插值方法一般用于数据较少的情况.数据拟合:不要求近似函数通过所有数据点,而是要求它能较好地反映数据的整体变化趋势。
共同点:插值与拟合都是根据实际中一组已知数据来构造一个能够反映数据变化规律的近似函数的方法,由于对近似要求的准则不同,因此二者在数学方法上有很大的差异.插值问题的一般提法:已知某函数)(x f y =(未知)的一组观测(或试验)数据),,2,1)(,(n i y x ii⋅⋅⋅=,要寻求一个函数)(x φ,使iiy x =)(φ),,2,1(n i ⋅⋅⋅=,则)()(x f x ≈φ.实际中,常常在不知道函数)(x f y =的具体表达式的情况下,对于i x x =有实验测量值iy y =),,2,1,0(n i ⋅⋅⋅=,寻求另一函数)(x φ使满足:)()(i i i x f y x ==φ),,2,1,0(n i ⋅⋅⋅=称此问题为插值问题,并称函数)(x φ为)(x f 的插值函数,nx x x x ,,,,21⋅⋅⋅称为插值节点,),,2,1,0()(n i y x ii⋅⋅⋅==φ称为插值条件,即)()(iiix f y x ==φ),,2,1,0(n i ⋅⋅⋅=,则)()(x f x ≈φ.(1) 拉格朗日(Lagrange )插值设函数)(x f y =在1+n 个相异点nx x x x ,,,,21⋅⋅⋅上的函数值为ny y y y ,,,,21⋅⋅⋅,要求一个次数不超过n 的代数多项式nnnx a x a x a a x P +⋅⋅⋅+++=221)(使在节点i x 上有),,2,1,0()(n i y x P ii n ⋅⋅⋅==成立,称之为n 次代数插值问题,)(x P n称为插值多项式.可以证明n 次代数插值是唯一的.事实上: 可以得到j n j n i i j in y x x xx x P j i ∑∏==⎥⎥⎦⎤⎢⎢⎣⎡⎪⎪⎭⎫⎝⎛--=≠00)()( 当1=n 时,有二点一次(线性)插值多项式:101001011)(y x x x x y x x x x x P --+--=当n =2时,有三点二次(抛物线)插值多项式: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)牛顿(Newton ) 插值牛顿插值的基本思想:由于)(x f y =关于二节点10,x x 的线性插值为)()()()()()()()()(00101000010101x x x x x f x f x p x x x x x f x f x f x p ---+=---+= 假设满足插值条件)2,1,0()()(2===i x p y x f iii的二次插值多项式一般形式为))(()()(1212x x x x c x x c c x p --+-+= 由插值条件可得⎪⎩⎪⎨⎧=--+-+=-+=)())(()()()()(21202202101011000x f x x x x c x x c c x f x x c c x f c 可以解出⎪⎪⎪⎩⎪⎪⎪⎨⎧------=--==020101121220101100)()()()()()(),(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所以))(()())(()()(10211020102x x x x c x p x x x x c x x c c x p --+=--+-+=类似的方法,可以得到三次插值多项式等,按这种思想可以得到一般的牛顿插值公式.函数的差商及其性质对于给定的函数)(x f ,用),,,(10n x x x f ⋅⋅⋅表示关于节点nx x x ,,,1⋅⋅⋅的n 阶差商,则有一阶差商:01011)()(),(x x x f x f x x f --=,121221)()(),(x x x f x f x x f --= 二阶差商:021021210),(),(),,(x x x x f x x f xx x f --=n 阶差商:0110211),,,(),,,(),,,(x x x x x f x x x f x x x f n n n n -⋅⋅⋅-⋅⋅⋅=⋅⋅⋅-差商有下列性质:(1)差商的分加性:∑∏=≠=-=⋅⋅⋅nk nk j j j kk n x xx f xx x f 0)(01)()(),,,(.(2)差商的对称性:在),,,(1nx x x f ⋅⋅⋅中任意调换jix x ,的次序其值不变.牛顿插值公式: 一次插值公式为))(,()()(01001x x x x f x f x p -+=二次插值公式为))()(,,()())()(,,())(,()()(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 --+=--+-+=于是有一般的牛顿插值公式为)())()(,,,()()())()(,,,())()(,,())(,()()(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可以证明:其余项为))(())()(,,,,()(11010n n n x x x x x x x x x x x x f x R --⋅⋅⋅--⋅⋅⋅=-实际上,牛顿插值公式是拉格朗日插值公式的一种变形,二者是等价的.另外还有著名的埃尔米特(Hermite )插值等.(3)样条函数插值方法样条,实质上就是由分段多项式光滑连接而成的函数,一般称为多项式样条.由于样条函数的特殊性质,决定了样条函数在实际中有着重要的应用.样条函数的一般概念定义 设给定区间],[b a 的一个分划b x x x a n=<⋅⋅⋅<<=∆1:,如果函数)(x s 满足条件:(1) 在每个子区间),,2,1](,[1n i x x ii ⋅⋅⋅=-上是k 次多项式; (2) )(x s 及直到k -1阶的导数在],[b a 上连续.则称)(x s 是关于分划△的一个k 次多项式样条函数,nx x x ,,,1⋅⋅⋅称为样条节点,121,,,-⋅⋅⋅n x x x 称为内节点,nx x ,0称为边界节点,这类样条函数的全体记作),(k S P∆,称为k 次样条函数空间.若),()(k S x s P∆∈,则)(x s 是关于分划△的k 次多项式样条函数.k 次多项式样条函数的一般形式为∑∑=-=+-+=ki n j k j jii k x x k i x x s 011)(!!)(βα其中),,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 jj kj kj在实际中最常用的是2=k 和3的情况,即为二次样条函数和三次样条函数. 二次样条函数:对于],[b a 上的分划b x x x a n=<⋅⋅⋅<<=∆1:,则)2,()(!2!2)(11222102∆βαααP n j j jS x x x x x s ∈-+++=∑-=+其中)1,2,1(,0,)()(22-=⎪⎩⎪⎨⎧<≥-=-+n j x x x x x x x x j j j j . 三次样条函数:对于],[b a 上的分划b x x xa n =<⋅⋅⋅<<=∆10:,则)3,()(!3!3!2)(1133322103∆βααααP n j j jS x x x x x x s ∈-++++=∑-=+其中)1,2,1(,0,)()(33-=⎪⎩⎪⎨⎧<≥-=-+n j x x x x x x x x jjj j .1 二次样条函数插值)2,()(2∆∈P S x s 中含有2+n 个待定常数,故应需要2+n 个插值条件,因此,二次样条插值问题可分为两类:问题(1):已知插值节点ix 和相应的函数值),,2,1,0(n i y i⋅⋅⋅=,以及端点0x (或n x )处的导数值0'y (或ny '),求)2,()(2∆∈PS x s 使得⎩⎨⎧'=''='⋅⋅⋅==))(()(),,2,1,0()(20022n n i i y x s y x s n i y x s 或(5.1)问题(2):已知插值节点ix 和相应的导数值),,2,1,0(n i y i⋅⋅⋅=',以及端点0x (或n x )处的函数值0y (或ny ),求)2,()(2∆∈P S x s 使得⎩⎨⎧==⋅⋅⋅='='))(()(),,2,1,0()(20022n n i i y x s y x s n i y x s 或(5.2)事实上,可以证明这两类插值问题都是唯一可解的.对于问题(1),由条件(5.1)⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧'=+='==-+++==++==++=∑-=00210211222102121211112020201002)(,,3,2,)(2121)(21)(21)(y x x s n j y x x x x x s yx x x s y x x x s j j i i j i jj j ααβααααααααα 引入记号T n ),,,,,(11210-=ββααα X 为未知向量,T nn y y y y ),,,,(10'= C 为已知向量, ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡---=-0010)(21)(21211)(212110211211021212212222211200x x x x x x x x x x x x x x x n n n n n A 于是,问题转化为求方程组C AX =的解Tn ),,,,,(1121-=ββααα X 的问题,即可得到二次样条函数的)(2x s 的表达式.对于问题(2)的情况类似.2.三次样条函数插值由于)3,()(3∆∈P S x s 中含有3+n 个待定系数,故应需要3+n 个插值条件,因此可将三次样条插值问题分为三类: 问题(1):已知插值节点jx 和相应的函数值),,2,1,0(n j y j⋅⋅⋅=,以及两个端点0x ,n x 处的导数值0'y ,ny ',求)3,()(3∆∈PS x s 使满足条件⎪⎩⎪⎨⎧='='⋅⋅⋅==),0()(),,1,0()(33n j y x s n j y x s j j j j(5.3)问题(2):已知插值节点jx 和相应的函数值),,2,1,0(n j y j⋅⋅⋅=,以及两个端点0x ,nx 处的二阶导数值0y '',n y '',求)3,()(3∆∈PS x s 使满足条件⎪⎩⎪⎨⎧=''=''⋅⋅⋅==),0()(),,1,0()(33n j y x s n j y x s j j j j(5.4)问题(3):类似地,求)3,()(3∆∈PSx 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(5.5)这三类插值问题的条件都是3+n 个,可以证明其解都是唯一的〔8〕.一般的求解方法可以仿照二次样条的情况处理方法,在这里给出一种更简单的方法.仅依问题(1)为例,问题(2)和问题(3)的情况类似处理.由于在)3,()(3∆PS x s ∈区间],[b a 上是一个分段光滑,且具有二阶连续导数的三次多项式,则在子区间],[1+j jx x 上)(3x s ''是线性函数,记),,,1,0)((3n j x s d jj =''=为待定常数.由拉格朗日插值公式可得nj x x h h x x d h x x d x s j j j jj j jj j ,,1,0,,)(1113=-=-+-=''+++显然jjj h d dx s -='''+13)(在],[1+j jx x上为常数.于是在],[1+j j x x 上有31233)(6)(2))(()(j jjj j j j j j x x h d d x x d x x x s y x s --+-+-'+=+(5.6)则当1+=j x x 时,由(5.6)式和问题(1)的条件得121231362)()(+++=-++'+=j j jj j j j j j j y h d d h d h x s y x s故可解得)2(6)(113+++--='j j j jjj j d d h h y y x s(5.7)将(5.7)式代入(5.6)式得)1,,1,0](,[,)(6)(2)()2(6)(1312113-=∈--+-+-⎥⎥⎦⎤⎢⎢⎣⎡+--+=++++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 jj j j jj j j j j j j j(5.8) 在],[1j j x x-上同样的有),,2,1](,[,)(6)(2)()2(6)(131112111111113n 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 =∈--+-+-⎥⎥⎦⎤⎢⎢⎣⎡+--+=------------(5.9) 根据)(3x s的一阶导数连续性,由(5.9)式得)()2(6)0(311113j j j j j j j j x s d d h h y y x s '=++-=-'---- 结合(5.7)式整理得⎪⎪⎭⎫ ⎝⎛---+=++++--+-+----11111111162j 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 引入记号⎪⎪⎭⎫ ⎝⎛---+=+=--+--111116,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(5.10)再由边界条件:nny x s y x s '=''=')(,)(33得⎪⎪⎩⎪⎪⎨⎧⎪⎪⎭⎫ ⎝⎛--'=+⎪⎪⎭⎫ ⎝⎛'--=+----111100010106262n n n n n n n h y y y h d d y h y y h d d(5.11)联立(5.10),(5.11)式得方程组C D A =⋅(5.12)其中⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡---=----2121212112112200n n n n a a a a a aA ,⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=-n n d d d d 110 D ,⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛--'⎪⎪⎭⎫ ⎝⎛'--=----111110001066n n n n n n hy y y h c c y h y y h C 由方程组(6.12)可以唯一解出),,1,0(n j d j=,代入(5.8)式就可以得三次样条函数)(3x s 的表达式.B样条函数插值方法磨光函数实际中的许多问题,往往是既要求近似函数(曲线或曲面)有足够的光滑性,又要求与实际函数有相同的凹凸性,一般插值函数和样条函数都不具有这种性质.如果对于一个特殊函数进行磨光处理生成磨光函数(多项式),则用磨光函数构造出样条函数作为插值函数,既有足够的光滑性,而且也具有较好的保凹凸性,因此磨光函数在一维插值(曲线)和二维插值(曲面)问题中有着广泛的应用.由积分理论可知,对于可积函数通过积分会提高函数的光滑度,因此,我们可以利用积分方法对函数进行磨光处理.定义 若)(x f 为可积函数,对于0>h ,则称积分⎰+-=22,1)(1)(hx h x h dt t f h x f为)(x f 的一次磨光函数,h 称为磨光宽度.同样的,可以定义)(x f 的k 次磨光函数为)1()(1)(22,1,>=⎰+--k dt t f h x f hx h x h k h k事实上,磨光函数)(,x f h k 比)(x f 的光滑程度要高,且当磨光宽度h 很小时)(,x f h k 很接近于)(x f .等距B样条函数对于任意的函数)(x f ,定义其步长为1的中心差分算子δ如下:⎪⎭⎫ ⎝⎛--⎪⎭⎫ ⎝⎛+=2121)(x f x f x f δ在此取0)(+=x x f ,则002121+++⎪⎭⎫ ⎝⎛--⎪⎭⎫ ⎝⎛+=x x x δ是一个单位方波函数(如图5-1),记0)(+=Ωx x δ.并取1=h ,对)(0x Ω进行一次磨光得++++-+++-+++--+-+=-=⎥⎥⎦⎤⎢⎢⎣⎡⎪⎭⎫ ⎝⎛--⎪⎭⎫ ⎝⎛+==⎰⎰⎰⎰)1(2)1(2121)()(11212100212101x x x dt t dt t dt t t dt t x x xx x x x x x ΩΩ显然)(1x Ω是连续的(如图5-2).)(1x Ωo1-1/2 0 1/2 x -1 0 1 x 图5-1图5-2类似地可得到k 次磨光函数为kk j jk j k j k x k C x ++=+⎪⎭⎫ ⎝⎛-++-=Ω∑21!)1()(11 实际上,可以证明:)(x kΩ是分段k 次多项式,且具有1-k 阶连续导数,其k 阶导数有2+k个间断点,记为)1,,2,1,0(21+⋅⋅⋅=+-=k j k j x j.从而可知)(x kΩ是对应于分划+∞<<⋅⋅⋅<<<-∞∆+110:k x x x 的k 次多项式样条函数,称之为基本样条函数,简称为k 次B样条.由于样条节点为)1,,2,1,0(21+⋅⋅⋅=+-=k j k j xj是等距的,故)(x k Ω又称为k 次等距B样条函数.对于任意函数)(x f 的k 次磨光函数,由归纳法可以得到 [4,8] :⎪⎭⎫⎝⎛+≤≤--Ω=⎰∞+∞--22)()(1)(1,h x t h x dt t f htx h x f k h k 特别地,当1)(=x f 时,有1)(11⎰+∞∞--=-dt htx hk Ω,从而1)(⎰+∞∞-=dx x k Ω,且当k ≥1时有递推关系⎥⎦⎤⎢⎣⎡⎪⎭⎫ ⎝⎛-Ω⎪⎭⎫ ⎝⎛---⎪⎭⎫ ⎝⎛+Ω⎪⎭⎫ ⎝⎛++=Ω--212121211)(11x x k x k x k x k k k一维等距B样条函数插值等距B样条函数与通常的样条如下的关系: 定理设有区间],[b a 的均匀分划nab h n j jh x x j -=⋅⋅⋅=+=),,,1,0(:0∆,则对任意 k 次样条函数),()(k S x S p k ∆∈都可以表示为B样条函数族1021-=-=⎭⎬⎫⎩⎨⎧⎪⎭⎫⎝⎛+---n j k j k k j h x x Ω的线性组合[14].根据定理 5.1,如果已知曲线上一组点()jjy x ,,其中),,1,0,0(0n j h jh x x j⋅⋅⋅=>+=,则可以构造出一条样条磨光曲线(即为B样条函数族的线性组合)⎪⎭⎫⎝⎛--=∑--=j h x x c x S n kj k j k 01)(Ω 其中)1,,1,(-⋅⋅⋅+--=n k k j c j为待定常数.用它来逼近曲线,既有较好的精度,又有良好的保凸性.实际中,最常用的是3=k 的情况,即一般形式为⎪⎭⎫ ⎝⎛--=∑+-=j h x x c x S n j j 01133)(Ω 其中3+n 个待定系数)1,,0,1(+⋅⋅⋅-=n j c j可以由三类插值条件确定.由插值条件(5.3)得()()()⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧'=-'='==-='=-'='∑∑∑+-=+-=+-=n n j j n i n j j i n j j y j n c h x S ni y j i c x S y j c h x S 113311330113031)(,,1,0,)(1)(ΩΩΩ(5.13)注意到)(3x Ω的局部非零性及其函数值:61)1(,32)0(33=±=ΩΩ,当2≥x 时0)(3=x Ω;且由)21()21()(223--+='x x x ΩΩΩ知,21)1(,0)0(33=±'='ΩΩ,当2≥x 时0)(3='x Ω.则(5.13)中的每一个方程中只有三个非零系数,具体的为⎪⎩⎪⎨⎧'=+-==++'=+-+-+--n n n i i i i y h c c n i y c c c y h c c 2,,1,0,6421111011(5.14)由方程组(5.14)容易求解出)1,,0,1(+⋅⋅⋅-=n j c j,即可得到三次样条函数)(3x S 表达式.类似地,由插值条件(5.4)得待定系数的)1,,0,1(+⋅⋅⋅-=n j c j所满足的方程组为⎪⎩⎪⎨⎧''=+-==++''=+-+-+--nn n n i i i i y h c c c n i y c c c y h c c c 21111021012,,1,0,642(5.15)由插值条件(5.5)得待定系数的)1,,0,1(+⋅⋅⋅-=n j cj所满足的方程组为⎪⎪⎩⎪⎪⎨⎧==++=-+---=-++-=-+-+-+-+--+--+--ni y c c c c c c c c c c c c c c c c c c c i i i i n n n n n n n n ,,1,0,640)()(2)(0)(0)(0)()(4)(1111011111111011(5.16)方程组(5.15),(5.16)也都是容易求解的.注:上述等距B样条插值公式也适用于近似等距的情形,但在端点0x 和n x 处误差可能较大,实际应用时,为了提高在端点0x 和nx 处的精度,可以适当向左右延拓几个节点.二维等距B样条函数插值设有空间曲面),(y x f z =(未知),如果已知二维等距节点()()τj y ih x y x ji++=0,,)0,(>τh 上的值为),,2,1,0;,,2,1,0(m j n i z ij⋅⋅⋅=⋅⋅⋅=,则相应的B样条磨光曲面的一般形式为⎪⎭⎫ ⎝⎛--⎪⎭⎫⎝⎛--=∑∑--=--=j y y i h x x c y x s l m lj k ij n ki τΩΩ0011),( 其中),,2,1,0;,,2,1,0(m j n i c ij⋅⋅⋅=⋅⋅⋅=为待定常数,l k ,可以取不同值,常用的也是2,=l k 和3的情形.这是一种具有良好保凸性的光滑曲面(函数),在工程设计中是常用的,但只能使用于均匀分划或近似均匀分划的情况.(4) 最小二乘拟合方法最小二乘拟合方法的思想:由于一般插值问题并不总是可解的(即当插值条件多于待定系数的个数时,其问题无解),同时,问题的插值条件本身一般是近似的,为此,只要求在节点上近似地满足插值条件,并使它们的整体误差最小,这就是最小二乘拟合法.最小二乘拟合方法可以分为线性最小二乘拟合方法和非线性最小二乘拟合方法.线性最小二乘拟合方法设{}m k kx 0)(=φ是一个线性无关的函数系,则称线性组合∑==mk k k x a x 0)()(φφ为广义多项式.如三角多项式:∑∑==+=mk k mk kkx b kx ax 0sin cos )(φ.设由给定的一组测量数据),(iiy x 和一组正数),,2,1(n i w i⋅⋅⋅=,求一个广义多项式∑==mk k k x a x 0)()(φφ使得目标函数[]21)(∑=-=ni i i i y x w S φ(5.17)达到最小,则称函数)(x φ为数据),,2,1)(,(n i y x ii⋅⋅⋅=关于权系数),,2,1(n i w i⋅⋅⋅=的最小二乘拟合函数,由于)(x φ关于待定系数ia 是线性的,故此问题又称为线性最小二乘问题. 注意:这里{}m k kx 0)(=φ可根据实际来选择,权系数iw 的选取更是灵活多变的,有时可选取1=i w ,或nw i 1=,对于nw i1=,则相应问题称为均方差的极小化问题.最小二乘拟合函数的求解要使最小二乘问题的目标函数(5.17)达到最小,则由多元函数取得极值的必要条件得),,2,1,0(0m k a Sk==∂∂ 即),,2,1,0(0)()(10m k x y x a w i k ni i m k i k k i ⋅⋅⋅⋅==⎥⎦⎤⎢⎣⎡-∑∑==φφ 亦即),,2,1,0()()()(001m k x y w a x x w n i i k i i j mj n i i k i j i ⋅⋅⋅⋅==⎥⎦⎤⎢⎣⎡∑∑∑===φφφ(5.18)是未知量为ma a a a ,,,,21⋅⋅⋅的线性方程组,称(5.18)式为正规方程组.实际中可适当选择函数系{}m k kx 0)(=φ,由正规方程组解出ma a a a ,,,,210⋅⋅⋅,于是可得最小二乘拟合函数∑==mk kk x a x 0)()(φφ.一般线性最小二乘拟合方法将上面一元函数的最小二乘拟合问题推广到多元函数,即为多维线性最小二乘拟合问题.假设已知多元函数),,,(21nx x x f y ⋅⋅⋅=的一组测量数据);,,,(21iniiiy x x x ⋅⋅⋅),,2,1(m i ⋅⋅⋅=和一组线性无关的函数系{}N k nk x x x 021),,,(=⋅⋅⋅φ,求函数∑=⋅⋅⋅=⋅⋅⋅Nk n k k n x x x a x xx 02121),,,(),,,(φφ对于一组正数mw w w ,,,21⋅⋅⋅,使得目标函数[]2121),,,(∑=⋅⋅⋅-=mi ni i i i i x x x y w S φ达到最小.其中待定系数N a a a a,,,,210⋅⋅⋅由正规方程组),,2,1,0(),(),(0N k y a Nj k j k j⋅⋅⋅==∑=φφφ确定,此处ini i i k mi i k ni i i k mi ni i i j i k j y x x x w y x x x x x x w ),,,(),(),,,(),,,(),(21121121⋅⋅⋅=⋅⋅⋅⋅⋅⋅=∑∑==φφφφφφ注:上面的函数φ关于ia 都是线性的,这就是线性最小二乘拟合问题,对于这类问题的正规组总是容易求解的.如果φ关于ia 是非线性的,则相应的问题称为非线性最小二乘拟合问题.非线性最小二乘拟合方法假设已知多元函数),,,(21nx x x f y ⋅⋅⋅=的一组测量数据);,,,(21iniiiy x x x ⋅⋅⋅),,2,1(m i ⋅⋅⋅=,要求一个关于参数),,2,1,0(N j a j⋅⋅⋅=是非线性的函数),,,;,,,(1021Nn a a a x x x ⋅⋅⋅⋅⋅⋅=φφ对一组正数mw w w ,,,21⋅⋅⋅使得目标函数[]21102110),,,;,,,(),,,(∑=⋅⋅⋅⋅⋅⋅-=⋅⋅⋅mi N ni i i i i N a a a x x x y w a a a S φ达到最小,则称之为非线性最小二乘问题.这类问题属于无约束的最优化问题,一般问题的求解是很复杂的,通常情况下,可以采用共轭梯度法、最速下降法、拟牛顿法和变尺度法等方法求解.实例:黄河小浪底调水调沙问题问题的提出2004年6月至7月黄河进行了第三次调水调沙试验,特别是首次由小浪底、三门峡和万家寨三大水库联合调度,采用接力式防洪预泄放水,形成人造洪峰进行调沙试验获得成功.整个试验期为20多天,小浪底从6月19日开始预泄放水,直到7月13日恢复正常供水结束.小浪底水利工程按设计拦沙量为75.5亿立方米,在这之前,小浪底共积泥沙达14.15亿吨.这次调水调试验一个重要目的就是由小浪底上游的三门峡和万家寨水库泄洪,在小浪底形成人造洪峰,冲刷小浪底库区沉积的泥沙.在小浪底水库开闸泄洪以后,从6月27日开始三门峡水库和万家寨水库陆续开闸放水,人造洪峰于29日先后到达小浪底,7月3日达到最大流量2700立方米/每秒,使小浪底水库的排沙量也不断地增加.下面是由小浪底观测站从6月29日到7月10日检测到的试验数据:表5-1: 试验观测数据单位:水流为立方米/每秒,含沙量为公斤/立方米·84··85·注:以上数据主要是根据媒体公开报道的结果整理而成的,不一定与真实数据完全相符.现在,根据试验数据建立数学模型研究下面的问题:(1) 给出估算任意时刻的排沙量及总排沙量的方法;(2) 确定排沙量与水流量的变化关系.模型的建立与求解对于问题(1),根据所给问题的试验数据,要计算任意时刻的排沙量,就要确定出排沙量随时间变化的规律,可以通过插值来实现.考虑到实际中排沙量应该是随时间连续变化的,为了提高精度,我们采用三次B样条函数进行插值.下面构造三次B样条函数)(x S y =.由试验数据,时间是每天的早8点和晚8点,间隔都是12个小时,共24个点)24,,2,1(⋅⋅⋅=i t i.为了计算方便,令)23,,,1,0(122128⋅⋅⋅=+⎥⎦⎤⎢⎣⎡⋅+-=i i t x i i(5.19)则it 对应于)23,,1,0(1⋅⋅⋅=+=i i x i.于是以)23,,1,0(⋅⋅⋅=i x i为插值节点(等距),步长1=h .其相应的排沙量为)23,,1,0(⋅⋅⋅=i y i 对应关系如下表:·86·表5-2: 插值数据对应关系单位:排沙量为公斤函数)(x S y =所满足的条件为 (1)23,,1,0,)(⋅⋅⋅==i y x S ii;(2) 3500)(,56400)(2223222323231212-=--≈'='=--≈'='x x y y x S y x xy yx S y .取)(x S 的三次B样条函数一般形式为∑-=⎪⎭⎫⎝⎛--=24103)(j j j h x x c x S Ω·87·其中)24,,1,0,1(⋅⋅⋅-=j cj为待定常数,1=h .在这里⎪⎪⎪⎩⎪⎪⎪⎨⎧≥<<+-+-≤+-=Ω2,021,342611,3221)(23233x x x x x x x x x且易知⎪⎪⎪⎩⎪⎪⎪⎨⎧≥±===Ω2,01,610,32)(3x x x x和⎪⎪⎩⎪⎪⎨⎧≥±===Ω'2,01,210,0)(3x x x x 根据B样条函数的性质,)(x S ''在[]23,x x 上连续,则有()∑-=--'='='2413)(j jj xx c x S y Ω由插值条件(1),(2)可得到下列方程组()()()⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧'=-'=''=-'='⋅⋅⋅==-=∑∑∑-=-=-=23241323024130241323)()(23,,1,0,)(y j c x S y j c x S i y j i c x S j j j j i j j i ΩΩΩ 即⎪⎩⎪⎨⎧'=+-'=+-⋅⋅⋅==++-+-23242311112223,,1,0,64y c c y c c i y c c c i i i i 将232324112,2y c c y c c '+='-=-代入前24个方程中的第一个和最后一个,便可得到方程组F AC =,其中·88·⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⋅⋅⋅⋅⋅⋅=⨯232102424,421410141014124c c c c C A ,⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡'-'+=3400048000684000458400266626232322100 y y y y y y F显然A 为满秩阵,方程组F AC =一定有解,用消元法求解可得问题的解为56044.39830=c , 4117111.2031=c , 2159510.7882=c , 9189845.6433=c ,1203106.6364=c , 8239727.8115=c ,8249182.1166=c , 1263543.7217=c ,9287842.9988=c , 2302284.2839=c ,4317419.86810=c , 1304836.24311=c ,3307635.15912=c ,6305423.11913=c ,2270672.36214=c ,4240287.43115=c ,0154177.91216=c ,4103000.92017=c ,99818.406218=c , 43725.454719=c ,49279.775020=c ,32155.445221=c , 2098.444222=c ,7450.777923=c ,-450.777924311.2034,2232324011='+=='-=-y c c y c c . 将)24,,1,0,1(⋅⋅⋅-=j c j代入()∑-=--==24131)(j jj x c x S y Ω(5.20)即得排沙量的变化规律.由(5.19)和(5.20)式可得到第i 时间段(12小时为一段)内,任意时刻]12,0[∈t 的排沙量.则总的排沙量为()dt j t c dx x S Y j j⎰∑⎰-=--Ω==284824132411)(经计算可得1110844.1⨯=Y 吨,即从6月29日至7月10日小浪底水库排沙总量大约为1.844亿吨,此与媒体报道的排沙量基本相符.对于问题(2),研究排沙量与水量的关系,从试验数据可以看出,开始排沙量是随着水流量的增加而增长,而后是随着水流量的减少而减少.显然,变化规律并非是线性的关系,为此,我们问题分为两部分,从开始水流量增加到最大值2720立方米/每秒(即增长的过程)为一段,从水流量的最大值到结束为第二段,分别来研究水流量与排沙量的关系.具体数据如表5-3和5-4.表5-3: 第一阶段试验观测数据 单位:水流为立方米/每秒,含沙量为公斤/立方米表5-4: 第二阶段试验观测数据单位:水流为立方米/每秒,含沙量为公斤/立方米对于第一阶段,由表5-3用Matlab作图(如图5-3)可以看出其变化趋势,我们用多项式作最小二乘拟合.·90··91·图5-3设拟合函数为∑==mk kk x a x 1)(φ确定待定常数),,1,0(m k ak=使得211111102])([∑∑∑===⎥⎦⎤⎢⎣⎡-=-=i i i m k k i k i i y x a y x S φ有最小值.于是可以得到正规方程组为m k x y a x mj i k i i j i j k i ,,1,0,0111111 ==⎪⎭⎫⎝⎛∑∑∑===+ 当3=m 时,即取三次多项式拟合,则3,2,1,0,1113111321112111110111==⎪⎭⎫⎝⎛+⎪⎭⎫ ⎝⎛+⎪⎭⎫ ⎝⎛+⎪⎭⎫ ⎝⎛∑∑∑∑∑==+=+=+=k x y a x a x a x a x i k i i i k i i k i i k i i k i求解可得73321108423.1,103172.1,3.1784,-2492.9318--⨯=⨯-===a a a a .于是可得拟合多项式为332213)(x a x a x a a x +++=φ,最小误差为847.72=S ,拟合效果如图所示.·92·图:三次拟合效果,带*号的为拟合曲线.类似地,当4=m 时,即取四次多项式拟合,则正规方程组为4,3,2,1,0111411143111321112111110111==⎪⎭⎫⎝⎛+⎪⎭⎫ ⎝⎛+⎪⎭⎫ ⎝⎛+⎪⎭⎫ ⎝⎛+⎪⎭⎫ ⎝⎛∑∑∑∑∑∑==+=+=+=+=k x y a x a x a x a x a x i ki i i k i i k i i k i i k i i k i求解可得104633210109312.1,1094.1,102626.7,12.0624,-7434.6557---⨯-=⨯=⨯-===a a a a a 于是可得拟合多项式为443322104)(x a x a x a x a a x ++++=φ,最小误差为102.66=S ,拟合效果如图5-5所示.图5-5:四次拟合效果,带*号的为拟合曲线.从上面的三次多项式拟合和四次多项拟合效果来看,差别不大.基本可以看出排沙量与水流量的关系.图5-6:第二段三·93··94· 次多项式拟合效果对于第二阶段,由表5-4可以类似地处理.我们用线性最小二乘法作三次和四多项式拟合.拟合效果如图5-6和5-7所示,最小误差分别为5.459=S 和1.236=S . 从拟合效果来看,显然四次多项式拟合要比三次多项式拟合好的多.图5-7:第二段四次多项式拟合效果。
插值法和拟合实验报告一、实验目的1.通过实验了解插值法和拟合法在数值计算中的应用;2.掌握拉格朗日插值法、牛顿插值法和分段线性插值法的原理和使用方法;3.学会使用最小二乘法进行数据拟合。
二、实验仪器和材料1.一台计算机;2. Matlab或其他适合的计算软件。
三、实验原理1.插值法插值法是一种在给定的数据点之间“插值”的方法,即根据已知的数据点,求一些点的函数值。
常用的插值法有拉格朗日插值法、牛顿插值法和分段线性插值法。
-拉格朗日插值法:通过一个n次多项式,将给定的n+1个数据点连起来,构造出一个插值函数。
-牛顿插值法:通过递推公式,将给定的n+1个数据点连起来,构造出一个插值函数。
-分段线性插值法:通过将给定的n+1个数据点的连线延长,将整个区间分为多个小区间,在每个小区间上进行线性插值,构造出一个插值函数。
2.拟合法拟合法是一种通过一个函数,逼近已知的数据点的方法。
常用的拟合法有最小二乘法。
-最小二乘法:通过最小化实际观测值与拟合函数的差距,找到最优的参数,使得拟合函数与数据点尽可能接近。
四、实验步骤1.插值法的实验步骤:-根据实验提供的数据点,利用拉格朗日插值法、牛顿插值法、分段线性插值法,分别求出要插值的点的函数值;-比较三种插值法的插值结果,评价其精度和适用性。
2.拟合法的实验步骤:-根据实验提供的数据点,利用最小二乘法,拟合出一个合适的函数;-比较拟合函数与实际数据点的差距,评价拟合效果。
五、实验结果与分析1.插值法的结果分析:-比较三种插值法的插值结果,评价其精度和适用性。
根据实验数据和插值函数的图形,可以判断插值函数是否能较好地逼近实际的曲线。
-比较不同插值方法的计算时间和计算复杂度,评价其使用的效率和适用范围。
2.拟合法的结果分析:-比较拟合函数与实际数据点的差距,评价拟合效果。
可以使用均方根误差(RMSE)等指标来进行评价。
-根据实际数据点和拟合函数的图形,可以判断拟合函数是否能较好地描述实际的数据趋势。
插值与拟合算法分析在数学与计算机科学领域,插值与拟合算法是两种常用的数据处理技术。
插值算法通过已知数据点之间的内插来估算未知数据点的值,而拟合算法则通过求取最佳拟合曲线或函数来逼近已知数据点。
本文将对插值与拟合算法进行详细分析,并比较它们在不同应用中的优缺点。
一、插值算法插值算法主要用于通过已知数据点之间的内插来估算未知数据点的值。
常用的插值算法包括拉格朗日插值、牛顿插值、样条插值等。
这些算法根据插值函数的不同特点,适用于不同类型的数据处理。
1. 拉格朗日插值拉格朗日插值是一种基于代数多项式的插值方法。
它通过构造一个全局多项式函数来拟合已知数据点,并推导出未知数据点的估算值。
拉格朗日插值算法具有简单易懂、计算效率高等优点,但在处理大量数据点时可能会出现龙格现象,导致插值结果有一定误差。
2. 牛顿插值牛顿插值是一种基于差商的插值方法。
它通过计算差商的递推关系,构造一个分段多项式函数来拟合已知数据点。
相比于拉格朗日插值,牛顿插值算法具有更高的数值稳定性和精度,并且可以方便地进行动态插值。
3. 样条插值样条插值是一种基于分段函数的插值方法。
它将整个数据区间划分为若干小段,并使用不同的插值函数对每一段进行插值。
样条插值算法通过要求插值函数的高阶导数连续,能够更好地逼近原始数据的曲线特征,因此在光滑性较强的数据处理中常被使用。
二、拟合算法拟合算法主要用于通过最佳拟合曲线或函数来逼近已知数据点。
常用的拟合算法包括最小二乘拟合、多项式拟合、非线性拟合等。
这些算法可以使拟合曲线与已知数据点尽可能地接近,从而进行更精确的数据分析和预测。
1. 最小二乘拟合最小二乘拟合是一种通过最小化残差平方和来求取最佳拟合曲线的方法。
它利用数据点与拟合曲线的差异来评估拟合效果,并通过求取最小残差平方和的参数值来确定拟合曲线的形状。
最小二乘拟合算法广泛应用于线性回归和曲线拟合等领域。
2. 多项式拟合多项式拟合是一种通过多项式函数来逼近已知数据点的方法。
试验二 插值法一、实验目的(1) 学会Lagrange 插值和牛顿插值等基本插值方法; (2) 讨论插值的Runge 现象,掌握分段线性插值方法(3) 学会Matlab 提供的插值函数的使用方法,会用这些函数解决实际问题。
二、实验要求(1) 按照题目要求完成实验内容; (2) 写出相应的Matlab 程序;(3) 给出实验结果(可以用表格展示实验结果); (4) 分析和讨论实验结果并提出可能的优化实验。
(5) 写出实验报告。
三、实验步骤1、用编好的Lagrange 插值法程序计算书本P66 的例1、用牛顿插值法计算P77的例1。
2、已知函数在下列各点的值为:试用4次牛顿插值多项式4()P x 对数据进行插值,根据{(,),0.20.08,0,1,2,,10i i i x y x i i =+=L },画出图形。
3、在区间[-1,1]上分别取10,20n =用两组等距节点对龙格函数21(),(11)125f x x x =-≤≤+作多项式插值,对不同n 值,分别画出插值函数及()f x 的图形。
4、下列数据点的插值可以得到平方根函数的近似,在区间[0,64]上作图。
(1) 用这9个点作8次多项式插值8()L x 。
5、在区间[0,10]画出sin()y x =的曲线,取插值节点,0,1,,10k x k k ==L 和节点处的函数值sin()k k y x =,作分段线性插值,并画出相应的折线图,将两图形绘在一张图上。
6、搜索Matlab 提供的三个数据拟合函数polyfit()、lsqcurvefit()、nlinfit()的使用方法。
列出它们的典型应用例子。
一、实验目的1、学会Lagrange 插值、牛顿插值和 分段线性插值等基本插值方法;2、讨论插值的Runge 现象,掌握分段线性插值方法3、学会Matlab 提供的插值函数、数据拟合函数的使用方法,并用这些函数解决实际问题。
二、实验原理1、拉格朗日插值多项式2、牛顿插值多项式3、分段线性插值三、实验步骤1、用MA TLAB 编写独立的拉格朗日插值多项式函数2、用MA TLAB 编写独立的牛顿插值多项式函数3、利用编写好的函数计算本章P66例1、P77例1的结果并比较。
实验二 插值与拟合实验名称:插值与拟合 实验类型: 验证性实验 学 时:23.1 实验环境① 操作系统:WindowsXP/Win7 ② 编程环境:自定3.2 实验目的① 掌握多项式插值法的基本思路和步骤; ② 了解整体插值的局限性及分段插值的基本思想。
③ 掌握最小二乘法拟合的基本原理和方法; ④ 培养运用计算机模拟解决问题的能力。
3.3 实验原理和方法3.3.1 多项式插值若n 次多项式() (0,1,,)k l x k n =在n +1个插值节点01...n x x x <<<上满足插值条件:1()()(,0,1,,)0()k i ik i k l x i k n i k δ=⎧===⎨≠⎩则称这n +1个n 次多项式01(), (), ... , ()n l x l x l x 为插值节点01, , ... , n x x x 上的n 次插值基函数。
由于i k ≠时,()=0k i l x ,故0111, , ... , , , ... , k k n x x x x x -+为()k l x 的零点,从而可以设0111()()()()()()k k k k n l x A x x x x x x x x x x -+=-----由()=1k k l x 可得 01111()()()()()k k k k k k k k n A x x x x x x x x x x -+=-----所以可得01110111()()()()()()(0,1,,)()()()()()k k n k k k k k k k k n x x x x x x x x x x l x k n x x x x x x x x x x -+-+-----==-----。
若记10()()nn i i x x x ω+==-∏,则有1'1()()()()n k k n k x l x x x x ωω++=-,从而有1'01()()()()()()()n nn n i i i i i n i x P x l x f x f x x x x ωω+==+==-∑∑。
上式是次数不超过n 的多项式,且满足所有的插值条件,称之为Lagrange 插值多项式。
3.3.2 数据拟合 最小二乘法基本原理 已知数据对()(),1,2,,j jx y j m =,求多项式0()()nii i p x a x n m ==<∑使得20110(,,,)mn i n i j j j i a a a a x y ==⎛⎫Φ=- ⎪⎝⎭∑∑为最小,这就是多项式拟合的最小二乘法。
最小二乘法的算法描述线性函数()p x a bx =+为例,拟合给定数据(),,1,2,,i i x y i m =。
算法描述:步骤1:输入m 值,及(),,1,2,,i i x y i m =。
步骤2:建立法方程组,并求解。
步骤3:输出()p x a bx =+。
3.4 实验内容和步骤3.4.1 多项式插值1. 给定sin110.190809,sin120.207912,sin130.224951,ooo===构造插值多项式计算'sin1130o 。
(1) 编程实现拉格朗日插值,并计算结果。
(2) 将计算结果和查表结果进行比较。
2. 区间[]5,5-作等距划分:101(0,1,,),k x kh k n h n=-+== ,以k x (0,1,,k n =)为节点对函数21()1f x x=+进行插值逼近。
(分别取5,10,20n =) (1) 用多项式插值对()f x 进行逼近,并在同一坐标系下作出函数的图形,进行比较。
写出插值函数对()f x 的逼近程度与节点个数的关系,并分析原因。
(2) 试用分段插值(任意选取)对()f x 进行逼近,在同一坐标下画出图形,观察分段插值函数对()f x 的逼近程度与节点个数的关系。
解:1.(1)拉格朗日插值程序代码:function f = Language(x,y,x0)syms t l;if(length(x) == length(y))n = length(x);elsedisp('x和y的维数不相等!');return; %检错endh=sym(0);for (i=1:n)l=sym(y(i));for(j=1:i-1)l=l*(t-x(j))/(x(i)-x(j));end;for(j=i+1:n)l=l*(t-x(j))/(x(i)-x(j));end;h=h+l;endsimplify(h);if(nargin == 3)f = subs (h,'t',x0); %计算插值点的函数值elsef=collect(h);f = vpa(f,6); %将插值多项式的系数化成6位精度的小数endx=[11,12,13];y=[0.190809,0.207912,0.224951];f=Language(x,y,11.5);vpa(f,8)程序结果:(2)查表所得的精确解:0.199367934,所以插值所得到的结果精确到小数点后五位。
解:2.(1)程序代码:function y=lagrange(x0,y0,x);n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;endx1=[-5:2:5];y=1./(1+x1.^2);x=[-5:0.001:5];y1=L(x1,y);x2=[-5:1:5];y=1./(1+x2.^2);y2=L(x2,y);x3=[-5:0.5:5];y=1./(1+x3.^2);y3=L(x3,y);y=1./(1+x3.^2);plot(x,subs(y1,'t',x),'r',x,subs(y2,'t',x),'g',x,subs(y3,'t',x),'b')hold onplot(x,1./(1+x.*x),'m')axis([-5 5 -0.5 2])xlabel('x')ylabel('y')legend('n=5插值后的图','n=10插值后的图','n=20插值后的图','原图',0)结果与原因分析:通过上图发现,高次插值稳定性差,而低次插值对于较大区间逼近精度又不够,而且,随着插值节点数目的增多,区域中间的确越来越逼近原函数,然而区域两端发生了较大震荡,即发生了Runge现象,这是由叉子余项的表达式所决定的。
解决这一矛盾的有效方法就是采用分段低次代数插值。
(2)分段插值函数程序代码:function output=L1syms f x lx;f=1/(1+x^2);N=input('请输入插值节点数N=');xx=-5:10/N:5;ff=zeros(1,length(xx));for i=1:(N+1)x=xx(i);ff(i)=eval(f);endM = -5:0.01:5;output=zeros(1,length(M));n = 1;for i=2:N+1for x=-5:0.01:5if x<xx(i) && x>=xx(i-1)lx(1)=ff(i-1)*(x-xx(i))/(xx(i-1)-xx(i)); lx(2)=ff(i)*(x-xx(i-1))/(xx(i)-xx(i-1)); output(n) = lx(1)+lx(2);n = n+1;endendendA =-5:0.01:5;Endoutput=L1;plot(A,output,'r');hold onoutput=L1;plot(A,output,'g');hold onoutput=L1;plot(A,output,'b');hold onezplot(f,[-5,5])xlabel('x')ylabel('y')legend('n=5分段差值后的图','n=10分段差值后的图','n=20分段插值后的图','原图',0)请输入插值节点数N=5请输入插值节点数N=10请输入插值节点数N=20程序结果:结果分析:由图可以发现,插值节点数目越多,插值函数越接近原函数,函数没有出现Runge现象,这是因为分段插值将区间分为了n份。
3.4.2 数据拟合1x 1 2 3 4 5i(1) (2) 求出其平方误差。
(1) 求以上数据形如2012()y x c c x c x =++的拟合曲线,及其平方误差。
(2) 求以上数据形如()bx y x ae =的拟合曲线,及其平方误差。
通过画出(1)(2)的图形,观察结果并结合其平方误差,写出你对数据拟合的认识。
解:1.(1)最小二乘法程序代码x=[1 2 3 4 5]; y=[4,4.5,6,8,8.5]; x1=0; x12=0; y1=0; xy=0; n=5; for i=1:n x1=x1+x(i); x12=x12+x(i)^2; y1=y1+y(i); xy=xy+x(i)*y(i); enda0=(y1*x12-x1*xy)/(n*x12-x1*x1); a1=(n*xy-x1*y1)/(n*x12-x1*x1); yy =a0+a1*x plot(x,y,'*') hold on plot(yy,'r') axis([0,6,0,9])结果:(2)平方误差:2.(1)MATLAB程序代码:xx=0;yy=0;x2=0;x3=0;x4=0;xy=0;y2=0;n=11; x1=[2,3,4,7,8,10,11,14,16,18,19];y1=[106.42,108.2,109.5,110,109.93,110.49,110.59,110.6,110.76,111,111.2]; for i=1:nxx=xx+x1(i);x2=x2+x1(i)^2;x3=x3+x1(i)^3;x4=x4+x1(i)^4;yy=yy+y1(i);xy=xy+x1(i)*y1(i);y2=y2+x1(i)*x1(i)*y1(i);endA=[n xx x2;xx x2 x3;x2 x3 x4];B=[yy;xy;y2];C=A\B;y2=C(1)+C(2)*x1+C(3)*x1.*x1;plot(x1,y1,'*')hold onplot(x1,y2,'r')结果:平方误差:程序代码:syms a bx=[2,3,4,7,8,10,11,14,16,18,19];fi=a+b./x;y1=[106.42,108.2,109.5,110,109.93,110.49,110.59,110.6,110.76,111,111.2]; y=log(y1);fy=fi-y;fy2=fy.^2;J=sum(fy.^2);Ja=diff(J,a); Jb=diff(J,b); Ja1=simple(Ja);Jb1=simple(Jb); [a,b]=solve(Ja1,Jb1);f=exp(a).*exp(b./x);plot(x,y1,'*')hold onplot(x,f,'r')xlabel('x')ylabel('y')legend('数值点','拟合后的图',0)>> s=0;for i=1:1:11s=s+(f(i)-y1(i)).^2;endvpa(s)结果平方误差拟合误差平方和为:0.47193209,通过计算比较发现,通过对数变换之后得到的拟合结果比进行二次函数拟合的结果要好很多。