郑州网站app开发,网页基础知识,网站建设费计入销售费用的子目,宁夏建网站报价声明#xff1a;以下笔记中的图片以及内容 均整理自“数学建模学习交流”清风老师的课程资料#xff0c;仅用作学习交流使用 文章目录 插值算法定义三个类型插值举例插值多项式分段插值三角插值 一般插值多项式原理拉格朗日插值法龙格现象分段线性插值 牛顿插值法 Hermite埃尔…声明以下笔记中的图片以及内容 均整理自“数学建模学习交流”清风老师的课程资料仅用作学习交流使用 文章目录 插值算法定义三个类型插值举例插值多项式分段插值三角插值 一般插值多项式原理拉格朗日插值法龙格现象分段线性插值 牛顿插值法 Hermite埃尔米特插值原理分段三次埃尔米特插值构造应用 三次样条插值定义应用 三次Hermite插值和三次样条插值的对比n维数据的插值 拟合算法最小二乘法拟合评价Matlab自带拟合工具箱cftool 插值算法
实际上本栏重点只有三次Hermite插值和三次样条插值的两行调用代码其他的全是废原话理
定义
设函数 y f ( x ) y f(x) yf(x)在区间 [ a , b ] [a,b] [a,b]上有定义且已知在点 a ≤ x 0 x 1 ⋯ x n ≤ b a \leq x_0 x_1 \cdots x_n \leq b a≤x0x1⋯xn≤b 上的值分别为 y 0 , y 1 , ⋯ , y n y_0,y_1,\cdots,y_n y0,y1,⋯,yn 若存在一简单函数 P x Px Px使 P ( x i ) y i ( i 0 , 1 , 2 ⋯ , n ) P(x_i) y_i\ (i 0,1,2\cdots,n)\ P(xi)yi (i0,1,2⋯,n) 则称 P ( x ) P(x) P(x)为 f ( x ) f(x) f(x)的插值函数点 x 0 , x 1 , ⋯ , x n x_0,x_1,\cdots,x_n x0,x1,⋯,xn称为插值节点包含插值节点的区间 [ a , b ] [a,b] [a,b]称为插值区间求插值函数 P ( x ) P(x) P(x)的方法称为插值法。
三个类型插值举例
插值多项式
示例已知函数 y f ( x ) y f(x) yf(x)在点 x 0 0 x_0 0 x00 x 1 1 x_1 1 x11 x 2 2 x_2 2 x22上的值分别为 f ( 0 ) 1 f(0)1 f(0)1 f ( 1 ) 2 f(1)2 f(1)2 f ( 2 ) 5 f(2)5 f(2)5。我们可以构造一个二次插值多项式 P ( x ) P(x) P(x)来逼近 f ( x ) f(x) f(x)。 设 P ( x ) a 0 a 1 x a 2 x 2 P(x)a_0 a_1x a_2x^2 P(x)a0a1xa2x2将已知点代入可得方程组 { a 0 a 1 × 0 a 2 × 0 2 1 a 0 a 1 × 1 a 2 × 1 2 2 a 0 a 1 × 2 a 2 × 2 2 5 \begin{cases}a_0 a_1\times0 a_2\times0^2 1\\a_0 a_1\times1 a_2\times1^2 2\\a_0 a_1\times2 a_2\times2^2 5\end{cases} ⎩ ⎨ ⎧a0a1×0a2×021a0a1×1a2×122a0a1×2a2×225 解这个方程组可得 a 0 1 a_0 1 a01 a 1 0 a_1 0 a10 a 2 1 a_2 1 a21所以插值多项式 P ( x ) 1 x 2 P(x)1 x^2 P(x)1x2。在区间 [ 0 , 2 ] [0,2] [0,2]上这个多项式可以用来近似原函数 f ( x ) f(x) f(x)。
分段插值
示例假设我们要对函数 f ( x ) 1 1 25 x 2 f(x)\frac{1}{1 25x^2} f(x)125x21在区间 [ − 1 , 1 ] [-1,1] [−1,1]上进行插值。如果我们只用一个高次多项式进行插值可能会出现龙格现象即在区间端点附近出现剧烈振荡。这时可以采用分段插值比如将区间 [ − 1 , 1 ] [-1,1] [−1,1]等分成若干个子区间如 [ − 1 , − 0.5 ] , [ − 0.5 , 0 ] , [ 0 , 0.5 ] , [ 0.5 , 1 ] [-1,-0.5],[-0.5,0],[0,0.5],[0.5,1] [−1,−0.5],[−0.5,0],[0,0.5],[0.5,1]。在每个子区间上分别进行低次多项式插值如线性插值或二次插值。 在区间 [ − 1 , − 0.5 ] [-1,-0.5] [−1,−0.5]上已知端点值 f ( − 1 ) f(-1) f(−1)和 f ( − 0.5 ) f(-0.5) f(−0.5)通过线性插值得到该区间上的插值函数 P 1 ( x ) P_1(x) P1(x)使得 P 1 ( − 1 ) f ( − 1 ) P_1(-1)f(-1) P1(−1)f(−1) P 1 ( − 0.5 ) f ( − 0.5 ) P_1(-0.5)f(-0.5) P1(−0.5)f(−0.5)并且在该区间内 P 1 ( x ) P_1(x) P1(x)近似 f ( x ) f(x) f(x)。同理在其他子区间上也进行类似的操作得到相应的插值函数 P 2 ( x ) P_2(x) P2(x)、 P 3 ( x ) P_3(x) P3(x)、 P 4 ( x ) P_4(x) P4(x)。这样就构成了整个区间 [ − 1 , 1 ] [-1,1] [−1,1]上的分段插值函数。
三角插值
示例对于一个周期为 2 π 2\pi 2π的函数 f ( x ) f(x) f(x)我们在区间 [ 0 , 2 π ] [0,2\pi] [0,2π]上取 n n n个等距节点 x k 2 k π n x_k \frac{2k\pi}{n} xkn2kπ k 0 , 1 , ⋯ , n − 1 k 0,1,\cdots,n - 1 k0,1,⋯,n−1已知这些节点上的函数值 f ( x k ) f(x_k) f(xk)。 三角插值多项式可以表示为 S n ( x ) a 0 2 ∑ k 1 n ( a k cos k x b k sin k x ) S_n(x)\frac{a_0}{2}\sum_{k 1}^{n}(a_k\cos kx b_k\sin kx) Sn(x)2a0∑k1n(akcoskxbksinkx)其中系数 a k a_k ak和 b k b_k bk通过以下公式计算 a k 1 π ∫ 0 2 π f ( x ) cos k x d x a_k\frac{1}{\pi}\int_{0}^{2\pi}f(x)\cos kx dx akπ1∫02πf(x)coskxdx k 0 , 1 , ⋯ , n k 0,1,\cdots,n k0,1,⋯,n b k 1 π ∫ 0 2 π f ( x ) sin k x d x b_k\frac{1}{\pi}\int_{0}^{2\pi}f(x)\sin kx dx bkπ1∫02πf(x)sinkxdx k 1 , 2 , ⋯ , n k 1,2,\cdots,n k1,2,⋯,n 例如假设 f ( x ) f(x) f(x)是一个简单的周期函数在 [ 0 , 2 π ] [0,2\pi] [0,2π]上 f ( x ) { x , 0 ≤ x π 2 π − x , π ≤ x 2 π f(x)\begin{cases}x, 0\leq x \pi\\2\pi - x, \pi\leq x 2\pi\end{cases} f(x){x,2π−x,0≤xππ≤x2π通过计算上述积分可以得到相应的三角插值多项式 S n ( x ) S_n(x) Sn(x)用它来近似原函数 f ( x ) f(x) f(x)在整个周期上的行为。
这些插值方法在不同的场景下各有优势分段插值可以避免高次多项式的振荡问题插值多项式计算相对简单三角插值适用于周期函数等情况它们都为函数的近似和数据的处理提供了有效的手段。
一般插值多项式
原理
定理 设有 n 1 n 1 n1个互不相同的节点 ( x i , y i ) (x_i,y_i) (xi,yi) ( i 0 , 1 , 2 , ⋯ , n ) (i 0,1,2,\cdots,n) (i0,1,2,⋯,n)则存在唯一的多项式 L n ( x ) a 0 a 1 x a 2 x 2 ⋯ a n x n L_n(x)a_0 a_1x a_2x^2 \cdots a_nx^n Ln(x)a0a1xa2x2⋯anxn使得 L n ( x j ) y j , ( j 0 , 1 , 2 , ⋯ , n ) L_n(x_j)y_j \ \ ,(j 0,1,2,\cdots,n) Ln(xj)yj ,(j0,1,2,⋯,n) 成立 即就是仅存在唯一多项式使其过给定的 n 1 n1 n1个点 证明利用范德蒙行列式 构造方程组 { a 0 a 1 x 0 a 2 x 0 2 ⋯ a n x 0 n y 0 a 0 a 1 x 1 a 2 x 1 2 ⋯ a n x 1 n y 1 ⋯ ⋯ ⋯ a 0 a 1 x n a 2 x n 2 ⋯ a n x n n y n \begin{cases} a_0 a_1x_0 a_2x_0^2 \cdots a_nx_0^n y_0 \\ a_0 a_1x_1 a_2x_1^2 \cdots a_nx_1^n y_1 \\ \cdots\cdots\cdots \\ a_0 a_1x_n a_2x_n^2 \cdots a_nx_n^n y_n \end{cases} ⎩ ⎨ ⎧a0a1x0a2x02⋯anx0ny0a0a1x1a2x12⋯anx1ny1⋯⋯⋯a0a1xna2xn2⋯anxnnyn令 A [ 1 x 0 ⋯ x 0 n 1 x 1 ⋯ x 1 n ⋮ ⋮ ⋯ ⋮ 1 x n ⋯ x n n ] A\begin{bmatrix} 1 x_0 \cdots x_0^n \\ 1 x_1 \cdots x_1^n \\ \vdots \vdots \cdots \vdots \\ 1 x_n \cdots x_n^n \end{bmatrix} A 11⋮1x0x1⋮xn⋯⋯⋯⋯x0nx1n⋮xnn X [ a 0 a 1 ⋯ a n ] T X\begin{bmatrix} \ a_0 \ a_1 \ \cdots \ a_n \end{bmatrix}^T X[ a0 a1 ⋯ an]T Y [ y 0 y 1 ⋯ y n ] T Y\begin{bmatrix} \ y_0 \ y_1 \ \cdots \ y_n \end{bmatrix}^T Y[ y0 y1 ⋯ yn]T方程组的矩阵形式 A X Y AX Y AXY由于 ∣ A ∣ ∏ i 1 n ∏ j 0 i − 1 ( x i − x j ) ≠ 0 |A|\prod_{i1}^{n}\prod_{j0}^{i - 1}(x_i - x_j) \neq 0 ∣A∣∏i1n∏j0i−1(xi−xj)0所以方程组有唯一解从而 L n ( x ) a 0 a 1 x a 2 x 2 ⋯ a n x n L_n(x)a_0 a_1x a_2x^2 \cdots a_nx^n Ln(x)a0a1xa2x2⋯anxn唯一存在. 注 只要 n 1 n 1 n1个节点互异满足上述插值条件的多项式是唯一存在的。如果不限制多项式的次数插值多项式并不唯一。
拉格朗日插值法
在数值分析中在若干个不同的地方得到相应的观测值拉格朗日插值法可以找到一个多项式其恰好在各个观测的点取到观测到的值。 两个点: ( x 0 , y 0 ) , ( x 1 , y 1 ) (x_0,y_0),(x_1,y_1) (x0,y0),(x1,y1) f ( x ) x − x 1 x 0 − x 1 y 0 x − x 0 x 1 − x 0 y 1 f(x)\frac{x - x_1}{x_0 - x_1}y_0\frac{x - x_0}{x_1 - x_0}y_1 f(x)x0−x1x−x1y0x1−x0x−x0y1 三个点: ( x 0 , y 0 ) , ( x 1 , y 1 ) , ( x 2 , y 2 ) (x_0,y_0),(x_1,y_1),(x_2,y_2) (x0,y0),(x1,y1),(x2,y2) f ( x ) ( x − x 1 ) ( x − x 2 ) ( x 0 − x 1 ) ( x 0 − x 2 ) y 0 ( x − x 0 ) ( x − x 2 ) ( x 1 − x 0 ) ( x 1 − x 2 ) y 1 ( x − x 0 ) ( x − x 1 ) ( x 2 − x 0 ) ( x 2 − x 1 ) y 2 f(x)\frac{(x - x_1)(x - x_2)}{(x_0 - x_1)(x_0 - x_2)}y_0\frac{(x - x_0)(x - x_2)}{(x_1 - x_0)(x_1 - x_2)}y_1\frac{(x - x_0)(x - x_1)}{(x_2 - x_0)(x_2 - x_1)}y_2 f(x)(x0−x1)(x0−x2)(x−x1)(x−x2)y0(x1−x0)(x1−x2)(x−x0)(x−x2)y1(x2−x0)(x2−x1)(x−x0)(x−x1)y2 拉格朗日插值多项式一般形式 记 ω n 1 ( x ) ( x − x 0 ) ( x − x 1 ) ⋯ ( x − x n ) \omega_{n 1}(x) (x - x_0)(x - x_1)\cdots(x - x_n) ωn1(x)(x−x0)(x−x1)⋯(x−xn) 则 ω n 1 ′ ( x k ) ( x k − x 0 ) ⋯ ( x k − x k − 1 ) ( x k − x k 1 ) ⋯ ( x k − x n ) . \omega_{n 1}(x_k) (x_k - x_0)\cdots(x_k - x_{k - 1})(x_k - x_{k 1})\cdots(x_k - x_n). ωn1′(xk)(xk−x0)⋯(xk−xk−1)(xk−xk1)⋯(xk−xn). 记 l i ( x ) ( x − x 0 ) ⋯ ( x − x i − 1 ) ( x − x i 1 ) ⋯ ( x − x n ) ( x i − x 0 ) ⋯ ( x i − x i − 1 ) ( x i − x i 1 ) ⋯ ( x i − x n ) l_i(x)\frac{(x - x_0)\cdots(x - x_{i - 1})(x - x_{i 1})\cdots(x - x_n)}{(x_i - x_0)\cdots(x_i - x_{i - 1})(x_i - x_{i 1})\cdots(x_i - x_n)} li(x)(xi−x0)⋯(xi−xi−1)(xi−xi1)⋯(xi−xn)(x−x0)⋯(x−xi−1)(x−xi1)⋯(x−xn) ( x − x 0 ) ⋯ ( x − x i − 1 ) ( x − x i ) ( x − x i 1 ) ⋯ ( x − x n ) ( x − x i ) ( x i − x 0 ) ⋯ ( x i − x i − 1 ) ( x i − x i 1 ) ⋯ ( x i − x n ) \frac{(x - x_0)\cdots(x - x_{i - 1})(x - x_i)(x - x_{i 1})\cdots(x - x_n)}{(x - x_i)(x_i - x_0)\cdots(x_i - x_{i - 1})(x_i - x_{i 1})\cdots(x_i - x_n)} (x−xi)(xi−x0)⋯(xi−xi−1)(xi−xi1)⋯(xi−xn)(x−x0)⋯(x−xi−1)(x−xi)(x−xi1)⋯(x−xn) ω n 1 ( x ) ( x − x i ) ω n 1 ′ ( x i ) \frac{\omega_{n 1}(x)}{(x - x_i)\omega_{n 1}(x_i)} (x−xi)ωn1′(xi)ωn1(x) 则拉格朗日插值多项式为 L n ( x ) ∑ i 0 n y i l i ( x ) L_n(x)\sum_{i 0}^{n}y_il_i(x) Ln(x)∑i0nyili(x) 即 L n ( x ) ∑ k 0 n y k ω n 1 ( x ) ( x − x k ) ω n 1 ′ ( x k ) . L_n(x)\sum_{k 0}^{n}y_k\frac{\omega_{n 1}(x)}{(x - x_k)\omega_{n 1}(x_k)}. Ln(x)∑k0nyk(x−xk)ωn1′(xk)ωn1(x).
但实际上拉格朗日插值法实际情况下并不常用因为有一个很大的问题——龙格现象。
龙格现象
龙格现象是指在使用多项式插值逼近函数时当插值节点等距分布且插值多项式的次数较高时在插值区间的两端会出现剧烈振荡的现象导致插值结果在区间端点附近与原函数偏差较大不能很好地逼近原函数。 所以在不清楚曲线具体运动趋势的情况下不要轻易使用高次插值。
分段线性插值
为了避免龙格现象带来的影响故我们使用分段低次插值。 这里以分段二次插值为例 选取跟节点 x x x最近的三个节点 x i − 1 , x i , x i 1 x_{i - 1},x_i,x_{i 1} xi−1,xi,xi1进行二次插值。 即在每一个区间 [ x i − 1 , x i 1 ] [x_{i - 1},x_{i 1}] [xi−1,xi1]上取 f ( x ) ≈ L 2 ( x ) ∑ k i − 1 i 1 [ y k ∏ j i − 1 j ≠ k i 1 ( x − x j ) ( x k − x j ) ] f(x)\approx L_2(x)\sum_{k i - 1}^{i 1}[y_k\prod_{\substack{j i - 1\\j\neq k}}^{i 1}\frac{(x - x_j)}{(x_k - x_j)}] f(x)≈L2(x)ki−1∑i1[ykji−1jk∏i1(xk−xj)(x−xj)]
这种分段的低次插值称为分段二次插值在几何上就是用分段抛物线代替 y f ( x ) y f(x) yf(x)故分段二次插值又称为分段抛物插值。
牛顿插值法
牛顿插值法每次插值只和前n项的值有关这样每次只要在原来的函数上添加新的项就能够产生新的函数。 f ( x ) f ( x 0 ) f [ x 0 , x 1 ] ( x − x 0 ) f(x)f(x_0)f[x_0,x_1](x - x_0) f(x)f(x0)f[x0,x1](x−x0) f [ x 0 , x 1 , x 2 ] ( x − x 0 ) ( x − x 1 ) ⋯ f[x_0,x_1,x_2](x - x_0)(x - x_1)\cdots f[x0,x1,x2](x−x0)(x−x1)⋯ f [ x 0 , x 1 , ⋯ , x n − 2 , x n − 1 ] ( x − x 0 ) ( x − x 1 ) ⋯ ( x − x n − 3 ) ( x − x n − 2 ) f[x_0,x_1,\cdots,x_{n - 2},x_{n - 1}](x - x_0)(x - x_1)\cdots(x - x_{n - 3})(x - x_{n - 2}) f[x0,x1,⋯,xn−2,xn−1](x−x0)(x−x1)⋯(x−xn−3)(x−xn−2) f [ x 0 , x 1 , ⋯ , x n − 1 , x n ] ( x − x 0 ) ( x − x 1 ) ⋯ ( x − x n − 2 ) ( x − x n − 1 ) f[x_0,x_1,\cdots,x_{n - 1},x_n](x - x_0)(x - x_1)\cdots(x - x_{n - 2})(x - x_{n - 1}) f[x0,x1,⋯,xn−1,xn](x−x0)(x−x1)⋯(x−xn−2)(x−xn−1) 差商定义: 称 f [ x 0 , x k ] f ( x k ) − f ( x 0 ) x k − x 0 f[x_0,x_k]\frac{f(x_k)-f(x_0)}{x_k - x_0} f[x0,xk]xk−x0f(xk)−f(x0)为函数 f ( x ) f(x) f(x)关于点 x 0 , x k x_0,x_k x0,xk的一阶差商亦称均差. 二阶差商: f [ x 0 , x 1 , x 2 ] f [ x 1 , x 2 ] − f [ x 0 , x 1 ] x 2 − x 0 f[x_0,x_1,x_2]\frac{f[x_1,x_2]-f[x_0,x_1]}{x_2 - x_0} f[x0,x1,x2]x2−x0f[x1,x2]−f[x0,x1] K阶差商: x k , x k − 1 x_k,x_{k - 1} xk,xk−1可以不相邻 f [ x 0 , x 1 , ⋯ , x k ] f [ x 1 , ⋯ , x k − 1 , x k ] − f [ x 0 , x 1 , ⋯ , x k − 1 ] x k − x 0 f[x_0,x_1,\cdots,x_k]\frac{f[x_1,\cdots,x_{k - 1},x_k]-f[x_0,x_1,\cdots,x_{k - 1}]}{x_k - x_0} f[x0,x1,⋯,xk]xk−x0f[x1,⋯,xk−1,xk]−f[x0,x1,⋯,xk−1]
与拉格朗日插值法相比牛顿插值法的计算过程具有继承性。但是牛顿插值也存在龙格现象的问题。同时两种插值仅仅要求插值多项式在插值节点处与被插函数有相等的函数值而这种插值多项式却不能全面反映被插值函数的性态。然而在许多实际问题中不仅要求插值函数与被插值函数在所有节点处有相同的函数值它也需要在一个或全部节点上插值多项式与被插函数有相同的低阶甚至高阶的导数值。
Hermite埃尔米特插值
原理
设函数 f ( x ) f(x) f(x)在区间 [ a , b ] [a,b] [a,b]上有 n 1 n 1 n1个互异节点 a x 0 x 1 x 2 ⋯ x n b a x_0 x_1 x_2 \cdots x_n b ax0x1x2⋯xnb定义在 [ a , b ] [a,b] [a,b]上函数 f ( x ) f(x) f(x)在节点上满足 f ( x i ) y i , f ′ ( x i ) y i ′ ( i 0 , 1 , 2 , ⋯ , n ) f(x_i)y_i,f(x_i)y_i\ (i 0,1,2,\cdots,n) f(xi)yi,f′(xi)yi′ (i0,1,2,⋯,n) 2 n 2 2n 2 2n2个条件可唯一确定一个次数不超过 2 n 1 2n 1 2n1的多项式 H 2 n 1 ( x ) H ( x ) H_{2n 1}(x)H(x) H2n1(x)H(x)满足 H ( x j ) y j , H ′ ( x j ) m j ( j 0 , 1 , ⋯ , n ) H(x_j)y_j,\ H(x_j)m_j\ (j 0,1,\cdots,n) H(xj)yj, H′(xj)mj (j0,1,⋯,n).其余项为 R ( x ) f ( x ) − H ( x ) f ( 2 n 2 ) ( ξ ) ( 2 n 2 ) ! ω 2 n 2 ( x ) R(x)f(x)-H(x)\frac{f^{(2n 2)}(\xi)}{(2n 2)!}\omega_{2n 2}(x) R(x)f(x)−H(x)(2n2)!f(2n2)(ξ)ω2n2(x)
可以看到直接使用Hermite插值得到的次数较高也存在龙格现象因此在实际应用中往往使用分段三次Hermite插值多项式(pchip).
分段三次埃尔米特插值
构造 已知给定 n 1 n 1 n1个节点 x 0 , x 1 , ⋯ , x n x_0,x_1,\cdots,x_n x0,x1,⋯,xn以及对应的函数值 y 0 , y 1 , ⋯ , y n y_0,y_1,\cdots,y_n y0,y1,⋯,yn和导数值 m 0 , m 1 , ⋯ , m n m_0,m_1,\cdots,m_n m0,m1,⋯,mn。 函数值基函数 h i ( x ) h_i(x) hi(x) 对于节点 x i x_i xi构造函数值基函数 h i ( x ) h_i(x) hi(x)使其满足 h i ( x j ) δ i j h_i(x_j)\delta_{ij} hi(xj)δij其中 δ i j \delta_{ij} δij是克罗内克符号当 i j i j ij时 δ i j 1 \delta_{ij}1 δij1当 i ≠ j i\neq j ij时 δ i j 0 \delta_{ij}0 δij0。一种常见的构造方式是 h i ( x ) ( 1 2 x − x i x i 1 − x i ) ( x − x i 1 x i − x i 1 ) 2 h_i(x)\left(1 2\frac{x - x_i}{x_{i 1}-x_i}\right)\left(\frac{x - x_{i 1}}{x_i - x_{i 1}}\right)^2 hi(x)(12xi1−xix−xi)(xi−xi1x−xi1)2当 i 0 , 1 , ⋯ , n − 1 i 0,1,\cdots,n - 1 i0,1,⋯,n−1 h n ( x ) ( 1 2 x − x n x n − 1 − x n ) ( x − x n − 1 x n − x n − 1 ) 2 h_n(x)\left(1 2\frac{x - x_n}{x_{n - 1}-x_n}\right)\left(\frac{x - x_{n - 1}}{x_n - x_{n - 1}}\right)^2 hn(x)(12xn−1−xnx−xn)(xn−xn−1x−xn−1)2。 导数值基函数 h ^ i ( x ) \hat{h}_i(x) h^i(x) 构造导数值基函数 h ^ i ( x ) \hat{h}_i(x) h^i(x)使其满足 h ^ i ′ ( x j ) δ i j \hat{h}_i(x_j)\delta_{ij} h^i′(xj)δij。例如 h ^ i ( x ) ( x − x i ) ( x − x i 1 x i − x i 1 ) 2 \hat{h}_i(x)(x - x_i)\left(\frac{x - x_{i 1}}{x_i - x_{i 1}}\right)^2 h^i(x)(x−xi)(xi−xi1x−xi1)2当 i 0 , 1 , ⋯ , n − 1 i 0,1,\cdots,n - 1 i0,1,⋯,n−1 h ^ n ( x ) ( x − x n ) ( x − x n − 1 x n − x n − 1 ) 2 \hat{h}_n(x)(x - x_n)\left(\frac{x - x_{n - 1}}{x_n - x_{n - 1}}\right)^2 h^n(x)(x−xn)(xn−xn−1x−xn−1)2。 构造三次埃尔米特插值多项式 H ( x ) H(x) H(x) H ( x ) ∑ i 0 n [ y i h i ( x ) m i h ^ i ( x ) ] H(x)\sum_{i 0}^{n}[y_ih_i(x)m_i\hat{h}_i(x)] H(x)∑i0n[yihi(x)mih^i(x)]
应用
Matlab有内置的函数直接调用就行了 p pchip(x,y,new_x) x,y是已知的样本点的坐标new_x是插值横坐标
x -pi:pi; y sin(x);
new_x -pi:0.1:pi;
p pchip(x,y,new_x);
plot(x, y, o, new_x, p, r-)三次样条插值
定义
设 y f ( x ) y f(x) yf(x)在点 x 0 , x 1 , x 2 , ⋯ x n x_0,x_1,x_2,\cdots x_n x0,x1,x2,⋯xn的值为 y 0 , y 1 , y 2 , ⋯ y n y_0,y_1,y_2,\cdots y_n y0,y1,y2,⋯yn若函数 S ( x ) S(x) S(x)满足下列条件 S ( x i ) f ( x i ) y i S(x_i)f(x_i)y_i S(xi)f(xi)yi i 0 , 1 , 2 , ⋯ , n i 0,1,2,\cdots,n i0,1,2,⋯,n在每个子区间 [ x i , x i 1 ] ( i 0 , 1 , 2 , ⋯ , n − 1 ) [x_i,x_{i 1}](i 0,1,2,\cdots,n - 1) [xi,xi1](i0,1,2,⋯,n−1)上 S ( x ) S(x) S(x)是三次多项式 S ( x ) S(x) S(x)在 [ a , b ] [a,b] [a,b]上二阶连续可微。
则称 S ( x ) S(x) S(x)为函数 f ( x ) f(x) f(x)的三次样条插值函数。 S ( x ) S(x) S(x)除了满足基本插值条件 s ( x i ) f i s(x_i)f_i s(xi)fi外还应具有如下形式 S ( x ) { S 0 ( x ) , x ∈ [ x 0 , x 1 ] , S 1 ( x ) , x ∈ [ x 1 , x 2 ] , ⋮ S n − 1 ( x ) , x ∈ [ x n − 1 , x n ] ; S(x)\begin{cases} S_0(x), x\in[x_0,x_1], \\ S_1(x), x\in[x_1,x_2], \\ \vdots \\ S_{n - 1}(x), x\in[x_{n - 1},x_n]; \end{cases} S(x)⎩ ⎨ ⎧S0(x),S1(x),⋮Sn−1(x),x∈[x0,x1],x∈[x1,x2],x∈[xn−1,xn]; S i ( x ) ∈ C 3 ( [ x i , x i 1 ] ) S_i(x)\in C^3([x_i,x_{i 1}]) Si(x)∈C3([xi,xi1]). 并且满足条件 { S i − 1 ( x i ) S i ( x i ) , S i − 1 ′ ( x i ) S i ′ ( x i ) , S i − 1 ′ ′ ( x i ) S i ′ ′ ( x i ) , \begin{cases} S_{i - 1}(x_i)S_i(x_i), \\ S_{i - 1}(x_i)S_i(x_i), \\ S_{i - 1}(x_i)S_i(x_i), \end{cases} ⎩ ⎨ ⎧Si−1(xi)Si(xi),Si−1′(xi)Si′(xi),Si−1′′(xi)Si′′(xi),
应用
Matlab有内置的函数 p spline (x,y, new_x) x,y是已知的样本点的坐标new_x是要插入处对应的横坐标
x -pi:pi;
y sin(x);
new_x -pi:0.1:pi;
p1 pchip(x,y,new_x); %分段三次埃尔米特插值
p2 spline(x,y,new_x); %三次样条插值
plot(x,y,o,new_x,p1,r-,new_x,p2,b-)
Legend(样本点,三次埃尔米特插值,三次样条插值,Location,SouthEast) %标注显示在东南方向三次Hermite插值和三次样条插值的对比 三次样条比Hermite更光滑实际应用中两种都可以使用。
n维数据的插值
这个了解即可。 p interpn(x1,x2,...,xn, y, new_x1,new_x2,...,new_xn, method) x1,x2,...,xn是已知的样本点的横坐标 y是已知的样本点的纵坐标 new_x1,new_x2,...,new_xn是要插入点的横坐标 method是要插值的方法 ‘linear’: 线性插值 (默认算法); ‘cubic’: 三次插值; ‘spline’: 三次样条插值法; (最为精准) ‘nearest’: 最邻近插值算法。
x -pi:pi; y sin(x);
new_x -pi:0.1:pi;
p interpn (x, y, new_x,spline);
%等价于p spline(x, y, new_x);
plot(x, y, o, new_x, p, r-)省流本栏最有用的两行 p pchip(x,y,new_x) p spline (x,y, new_x) 拟合算法
插值算法中得到的多项式 f ( x ) f(x) f(x)要经过所有样本点。但是如果样本点太多那么这个多项式次数过高会造成龙格现象。 比如说这个例子 尽管我们可以选择分段的方法避免龙格现象但是更多时候我们更倾向于得到一个确定的曲线尽管这条曲线不能经过每一个样本点但只要保证误差足够小即可这就是拟合的思想。拟合的结果是得到一个确定的较为简单的曲线
最小二乘法
高中学的 懒得解释。 解释一下 k , b arg min k , b ( f ( x ) ) k,b\argmin_{\substack{k,b}}(f(x)) k,bk,bargmin(f(x))指的是求使得 f ( x ) f(x) f(x)取最小值时的参数 k , b k,b k,b的值
推导过程 结论 k ^ n ∑ i 1 n x i y i − ∑ i 1 n y i ∑ i 1 n x i n ∑ i 1 n x i 2 − ∑ i 1 n x i ∑ i 1 n x i \hat{k} \frac{n\sum_{i 1}^{n}x_iy_i - \sum_{i 1}^{n}y_i\sum_{i 1}^{n}x_i}{n\sum_{i 1}^{n}x_i^2 - \sum_{i 1}^{n}x_i\sum_{i 1}^{n}x_i} k^n∑i1nxi2−∑i1nxi∑i1nxin∑i1nxiyi−∑i1nyi∑i1nxi b ^ ∑ i 1 n x i 2 ∑ i 1 n y i − ∑ i 1 n x i ∑ i 1 n x i y i n ∑ i 1 n x i 2 − ∑ i 1 n x i ∑ i 1 n x i \hat{b} \frac{\sum_{i 1}^{n}x_i^2\sum_{i 1}^{n}y_i - \sum_{i 1}^{n}x_i\sum_{i 1}^{n}x_iy_i}{n\sum_{i 1}^{n}x_i^2 - \sum_{i 1}^{n}x_i\sum_{i 1}^{n}x_i} b^n∑i1nxi2−∑i1nxi∑i1nxi∑i1nxi2∑i1nyi−∑i1nxi∑i1nxiyi
拟合评价 三个基本概念 总体平方和 S S T SST SST (Total sum of squares) S S T ∑ i 1 n ( y i − y ˉ ) 2 SST \sum_{i 1}^{n}(y_i - \bar{y})^2 SST∑i1n(yi−yˉ)2 误差平方和 S S E SSE SSE (The sum of squares due to error) S S E ∑ i 1 n ( y i − y ^ i ) 2 SSE \sum_{i 1}^{n}(y_i - \hat{y}_i)^2 SSE∑i1n(yi−y^i)2 回归平方和 S S R SSR SSR (Sum of squares of the regression) S S R ∑ i 1 n ( y ^ i − y ˉ ) 2 SSR \sum_{i 1}^{n}(\hat{y}_i - \bar{y})^2 SSR∑i1n(y^i−yˉ)2 可以证明 S S T S S E S S R SST SSE SSR SSTSSESSR 证明 拟合优度 0 ≤ R 2 S S R S S T S S T − S S E S S T 1 − S S E S S T ≤ 1 0 \leq R^2 \frac{SSR}{SST} \frac{SST - SSE}{SST} 1 - \frac{SSE}{SST} \leq 1 0≤R2SSTSSRSSTSST−SSE1−SSTSSE≤1 R 2 R^2 R2越接近 1 1 1说明误差平方和 S S E SSE SSE 越接近 0 0 0说明拟合的越好。 注意 R 2 R^2 R2只能用于拟合函数是线性函数时拟合结果的评价. 线性函数和其他函数如复杂的指数函数比较拟合的好坏直接看 S S E SSE SSE即可 线性函数的界定这里我们说的线性函数指的是对参数为线性即在函数中参数仅以一次方出现并且不能乘以或除以其他任何参数也不能出现参数的复合函数形式。 形如 y a b x 2 yabx^2 yabx2 是线性函数这里的参数就是 a , b a,b a,b. 形如 y a b 3 x y a b^{3}x yab3x、 y a b x b c x 2 y a bx bcx^{2} yabxbcx2、 y a ( x − b ) 2 y a(x - b)^{2} ya(x−b)2、 y a sin ( b c x ) y a\sin(b cx) yasin(bcx)都不是线性函数不能用 R 2 R^{2} R2.
实际应用中常使用多个不同的拟合函数去试然后比较这几个拟合函数的拟合优度取最优的一个拟合函数。 曲线复杂度与拟合误差之间的取舍次数越高或者指数函数越复杂误差越小但是产生龙格现象的概率越大且函数曲线越复杂。一个理解在插值算法中曲线经过每一个点那么曲线的 S S E SSE SSE即为0如果放在拟合算法中来讲这种曲线的拟合优度极高但是曲线本身十分复杂这就违背了拟合的初衷 因此实际上如果两个拟合函数的 S S E SSE SSE相差不大但是曲线复杂度相差较大那么就取明显更为简单的那个函数哪怕它的 S S E SSE SSE相较另一个更大一点。
y\_hat k*xb; % y的拟合值
SSR sum((y\_hat-mean(y)).^2) % 回归平方和 注:mean()是求均值的函数。
SSE sum((y\_hat-y).^2) % 误差平方和
SST sum((y-mean(y)).^2) % 总体平方和
SST-SSE-SSR
R\_2 SSR / SSTMatlab自带拟合工具箱cftool 拟合方法 Custom Equations用户自定义的函数类型 注意修改范围Fit optionsExponential指数逼近有2种类型 a ∗ e x p ( b ∗ x ) a*exp(b*x) a∗exp(b∗x) 、 a ∗ e x p ( b ∗ x ) c ∗ e x p ( d ∗ x ) a*exp(b*x) c*exp(d*x) a∗exp(b∗x)c∗exp(d∗x)Fourier傅立叶逼近有7种类型基础型是 a 0 a 1 ∗ c o s ( x ∗ w ) b 1 ∗ s i n ( x ∗ w ) a0 a1*cos(x*w) b1*sin(x*w) a0a1∗cos(x∗w)b1∗sin(x∗w)Gaussian高斯逼近有8种类型基础型是 a 1 ∗ e x p ( − ( ( x − b 1 ) / c 1 ) 2 ) a1*exp(-((x-b1)/c1)^2) a1∗exp(−((x−b1)/c1)2)Interpolant插值逼近有4种类型linear、nearest neighbor、cubic spline、shape-preservingPolynomial多形式逼近有9种类型linear ~、quadratic ~、cubic ~、4-9th degree ~Power幂逼近有2种类型 a ∗ x b 、 a ∗ x b c a*x^b 、a*x^b c a∗xb、a∗xbcRational有理数逼近分子、分母共有的类型是linear ~、quadratic ~、cubic ~、4-5th degree ~此外分子还包括constant型Smoothing Spline平滑逼近Sum of Sin Functions正弦曲线逼近有8种类型基础型是 a 1 ∗ s i n ( b 1 ∗ x c 1 ) a1*sin(b1*x c1) a1∗sin(b1∗xc1)Weibull只有一种 a ∗ b ∗ x ( b − 1 ) ∗ e x p ( − a ∗ x b ) a*b*x^{(b-1)}*exp(-a*x^b) a∗b∗x(b−1)∗exp(−a∗xb) 拟合工具还可以直接生成代码Generate code可放在附录里。 小tips导出图像可以在文件-导出设置-渲染-修改分辨率使图更清晰
