数学建模:插值与拟合简介

插值

多项式插值

给定 n + 1 n+1 n+1 个互异点 a ⩽ x 0 < x 1 . . . < x n ⩽ b a \leqslant {x_0} < {x_1}... < {x_n} \leqslant b ax0<x1...<xnb 的值 y 0 , y 1 , ⋯   , y n y_0, y_1 , \cdots , y_n y0,y1,,yn, 求一个次数不超过 n n n 次的多项式 p ( x ) = a 0 + a 1 x + . . . + a n x n p(x) = {a_0} + {a_1}x + ... + {a_n}{x^n} p(x)=a0+a1x+...+anxn, s . t . s.t. s.t.

p ( x i ) = y i ( i = 0 , … n ) p(x_i) = y_i (i = 0, … n) p(xi)=yi(i=0,n)

p ( x ) p(x) p(x) f ( x ) f(x) f(x) n n n 次插值多项式

Lagrange 插值

l k ( x ) = ( x − x 0 ) ( x − x 1 ) ⋯ ( x − x k − 1 ) ( x − x k + 1 ) ⋯ ( x − x n ) ( x k − x 0 ) ( x k − x 1 ) ⋯ ( x k − x k − 1 ) ( x k − x k + 1 ) ⋯ ( x k − x n ) l_{k}(x)=\frac{\left(x-x_{0}\right)\left(x-x_{1}\right) \cdots\left(x-x_{k-1}\right)\left(x-x_{k+1}\right) \cdots\left(x-x_{n}\right)}{\left(x_{k}-x_{0}\right)\left(x_{k}-x_{1}\right) \cdots\left(x_{k}-x_{k-1}\right)\left(x_{k}-x_{k+1}\right) \cdots\left(x_{k}-x_{n}\right)} lk(x)=(xkx0)(xkx1)(xkxk1)(xkxk+1)(xkxn)(xx0)(xx1)(xxk1)(xxk+1)(xxn)

∏ i ≠ k x − x i x k − x i \prod_{{ i \neq k}} \frac{x-x_{i}}{x_{k}-x_{i}} i=kxkxixxi

k = 0 , 1 , ⋯   , n k=0, 1, \cdots , n k=0,1,,n

满足

l k ( x i ) = δ k i = { 1 ,   i = k 0 ,   i ≠ k {l_k}({x_i}) = {\delta _{ki}} = \left\{ \begin{gathered} 1,{\text{ }}i = k \\ 0,{\text{ }}i \ne k \\ \end{gathered} \right. lk(xi)=δki={1, i=k0, i=k

l 0 ( x ) , l 1 ( x ) , ⋯   , l n ( x ) l_0(x), l_1(x), \cdots , l_n(x) l0(x),l1(x),,ln(x) 线性无关,即它们构成线性空间 P n ( x ) P_n(x) Pn(x) 的一组基

p ( x ) = y 0 l 0 ( x ) + y 1 l 1 ( x ) + ⋯ + y n l n ( x ) p(x)=y_{0} l_{0}(x)+y_{1} l_{1}(x)+\cdots+y_{n} l_{n}(x) p(x)=y0l0(x)+y1l1(x)++ynln(x)

满足插值条件 p ( x i ) = y i   ( i = 0 , … n ) p(x_i) = y_i \ (i = 0, … n) p(xi)=yi (i=0,n), 称为 拉格朗日插值多项式,记作 L n ( x ) L_n(x) Ln(x),即

L n ( x ) = ∑ j = 0 n y j l j ( x ) = ∑ j = 0 n y j ∏ i = 0 i ≠ j n x − x i x j − x i L_{n}(x)=\sum_{j=0}^{n} y_{j} l_{j}(x)=\sum_{j=0}^{n} y_{j} \prod_{{i=0 \\ i \neq j}}^{n} \frac{x-x_{i}}{x_{j}-x_{i}} Ln(x)=j=0nyjlj(x)=j=0nyji=0i=jnxjxixxi

三次样条插值

样条函数: 由一些按照某种光滑条件分段拼接起来的多项式组成的函数, 最常用的样条函数为三次样条函数,即由三次多项式组成,满足处处有二阶连续导数.

节点: x 0 < x 1 < ⋯ < x n − 1 < x n x_0 < x_1 < \cdots < x_n-1 < x_n x0<x1<<xn1<xn,
函数值: y i ,   i = 0 , 1 , 2 , ⋯   , n y_i, \ i = 0, 1, 2, \cdots , n yi, i=0,1,2,,n

确定三次样条函数 s ( x ) s(x) s(x) 满足 s ( x i ) = y i   ( i = 0 , 1 , 2 , … , n ) s\left(x_{i}\right)=y_{i} \ (i=0,1,2, \ldots, n) s(xi)=yi (i=0,1,2,,n)

由定义可设

s ( x ) = { s 1 ( x ) , x ∈ [ x 0 , x 1 ] s 2 ( x ) , x ∈ [ x 1 , x 2 ] ⋮ s n ( x ) , x ∈ [ x n − 1 , x n ] s(x)=\left\{\begin{array}{cl} s_{1}(x), & x \in\left[x_{0}, x_{1}\right] \\ s_{2}(x), & x \in\left[x_{1}, x_{2}\right] \\ \vdots & \\ s_{n}(x), & x \in\left[x_{n-1}, x_{n}\right] \end{array}\right. s(x)=s1(x),s2(x),sn(x),x[x0,x1]x[x1,x2]x[xn1,xn]

其中 s k ( x ) s_{k}(x) sk(x) [ x k − 1 , x k ] \left[x_{k-1}, x_{k}\right] [xk1,xk] 上的三次多项式, 且满足
s k ( x k − 1 ) = y k − 1 , s k ( x k ) = y k ( k = 1 , 2 , … , n ) s_{k}\left(x_{k-1}\right)=y_{k-1}, s_{k}\left(x_{k}\right)=y_{k} \quad(k=1,2, \ldots, n) sk(xk1)=yk1,sk(xk)=yk(k=1,2,,n)

三次样条函数的连接条件:

s ( x ) ∈ C 2 [ a , b ] s(x) \in C^{2}[a, b] s(x)C2[a,b]

⟹ s ′ ( x k − ) = s ′ ( x k + ) , s ′ ′ ( x k − ) = s ′ ′ ( x k + ) \Longrightarrow s^{\prime}\left(x_{k}^{-}\right)=s^{\prime}\left(x_{k}^{+}\right), s^{\prime \prime}\left(x_{k}^{-}\right)=s^{\prime \prime}\left(x_{k}^{+}\right) s(xk)=s(xk+),s(xk)=s(xk+)

⟹ s k ′ ( x k − ) = s k + 1 ′ ( x k + ) , s k ′ ′ ( x k − ) = s k + 1 ′ ′ ( x k + ) \Longrightarrow s_{k}^{\prime}\left(x_{k}^{-}\right)=s_{k+1}^{\prime}\left(x_{k}^{+}\right), s_{k}^{\prime \prime}\left(x_{k}^{-}\right)=s_{k+1}^{\prime \prime}\left(x_{k}^{+}\right) sk(xk)=sk+1(xk+),sk(xk)=sk+1(xk+)

k = 1 , 2 , … , n − 1 k=1,2, \ldots, n-1 k=1,2,,n1

每个 s k ( x ) s_k(x) sk(x) 均为三次多项式,有 4 4 4 个待定系数,所以共有 4 n 4_n 4n 个待定系数,需 4 n 4_n 4n 个方程才能确定。前面已经得到 2 n + 2 ( n − 1 ) = 4 n – 2 2n +2(n -1) = 4n –2 2n+2(n1)=4n2 个方程,还缺 2 2 2 个方程,实际问题通常对样条函数在端点处的状态有要求,即所谓的边界条件

三次样条函数的边界条件:

  • 第一类边界条件:给定函数在端点处的一阶导数,即

s ′ ( x 0 ) = f 0 ′ , s ′ ( x n ) = f n ′ s^{\prime}\left(x_{0}\right)=f_{0}^{\prime}, s^{\prime}\left(x_{n}\right)=f_{n}^{\prime} s(x0)=f0,s(xn)=fn

  • 第二类边界条件:给定函数在端点处的二阶导数,即

s ′ ′ ( x 0 ) = f 0 ′ ′ , s ′ ′ ( x n ) = f n ′ ′ s^{\prime \prime}\left(x_{0}\right)=f_{0}^{\prime \prime}, s^{\prime \prime}\left(x_{n}\right)=f_{n}^{\prime \prime} s(x0)=f0,s(xn)=fn

s ′ ′ ( x 0 ) = s ′ ′ ( x n ) = 0 s''({x_0}) = s''({x_n}) = 0 s(x0)=s(xn)=0 时, 称为自然边界条件,此时的样条函数称为自然样条函数

  • 第三类边界条件:设 f ( x ) f (x) f(x) 是周期函数,并设 x n – x 0 x_n – x_0 xnx0 是一个周期,于是 s ( x ) s(x) s(x) 满足

s ′ ( x 0 ) = s ′ ( x n ) , s ′ ′ ( x 0 ) = s ′ ′ ( x n ) s^{\prime}\left(x_{0}\right)=s^{\prime}\left(x_{n}\right), s^{\prime \prime}\left(x_{0}\right)=s^{\prime \prime}\left(x_{n}\right) s(x0)=s(xn),s(x0)=s(xn)

拟合

  • 插值: 过点 (适合精确数据)
  • 拟合: 不过点, 整体近似;(适合有经验公式或有误差的数据)

直线最小二乘拟合

对于给定的数据点 ( x i , y i ) \left(x_{i}, y_{i}\right) (xi,yi), i = 1 , 2 , ⋯   , N i=1,2, \cdots, N i=1,2,,N y = a + b x y=a+b x y=a+bx, 使总误差平方和为最小, 即

Q = ∑ i = 1 N [ y i − ( a + b x i ) ] 2 Q=\sum_{i=1}^{N}\left[y_{i}-\left(a+b x_{i}\right)\right]^{2} Q=i=1N[yi(a+bxi)]2

最小. Q Q Q 是关于未知数 a a a b b b 的二元函数, 这一问题就是要确定 a a a b b b 取何值时,二元函数 Q ( a , b ) Q(a,b) Q(a,b) 的值最小.

∂ Q ∂ a = 0 , ∂ Q ∂ b = 0 \frac{\partial \boldsymbol{Q}}{\partial a}=\mathbf{0}, \quad \frac{\partial \boldsymbol{Q}}{\partial b}=\mathbf{0} aQ=0,bQ=0

∂ Q ∂ a = ∑ i = 1 N [ y i − ( a + b x i ) ] ⋅ ( − 1 ) = 0 ⟹ N a + b ∑ i = 1 N x i = ∑ i = 1 N y i \frac{\partial Q}{\partial a}=\sum_{i=1}^{N}\left[y_{i}-\left(a+b x_{i}\right)\right] \cdot(-1)=0 \Longrightarrow N a+b \sum_{i=1}^{N} x_{i}=\sum_{i=1}^{N} y_{i} aQ=i=1N[yi(a+bxi)](1)=0Na+bi=1Nxi=i=1Nyi

∂ Q ∂ b = ∑ i = 1 N [ y i − ( a + b x i ) ] ⋅ ( − x i ) = 0 ⟹ a ∑ i = 1 N x i + b ∑ i = 1 N x i 2 = ∑ i = 1 N x i y i \frac{\partial Q}{\partial b}=\sum_{i=1}^{N}\left[y_{i}-\left(a+b x_{i}\right)\right] \cdot\left(-x_{i}\right)=0 \Longrightarrow a \sum_{i=1}^{N} x_{i}+b \sum_{i=1}^{N} x_{i}^{2}=\sum_{i=1}^{N} x_{i} y_{i} bQ=i=1N[yi(a+bxi)](xi)=0ai=1Nxi+bi=1Nxi2=i=1Nxiyi

从而,未知数 a , b a, b a,b 满足如下方程组 (正则方程组)

{ N a + b ∑ i = 1 N x i = ∑ i = 1 N y i a ∑ i = 1 N x i + b ∑ i = 1 N x i 2 = ∑ i = 1 N x i y i \begin{cases} N a+b \sum_{i=1}^{N} x_{i}=\sum_{i=1}^{N} y_{i} \\ a \sum_{i=1}^{N} x_{i}+b \sum_{i=1}^{N} x_{i}^{2}=\sum_{i=1}^{N} x_{i} y_{i} \end{cases} {Na+bi=1Nxi=i=1Nyiai=1Nxi+bi=1Nxi2=i=1Nxiyi

解出 a , b a, b a,b, 就可得最小二乘拟合直线 y = a + b x y=a+bx y=a+bx

一般线性最小二乘

给定一个线性无关的函数系 ${ {\varphi _k}(x) \ | \ k = 1,2, \cdots ,m} $ ,如果拟合函数以其线性组合的形式

f ( x ) = ∑ k = 1 m a k φ k ( x ) f(x) = \sum\limits_{k = 1}^m {{a_k}{\varphi _k}(x)} f(x)=k=1makφk(x)

出现,则 f ( x ) = f ( x , a 1 , a 2 , ⋯   , a m ) f(x) = f(x,{a_1},{a_2}, \cdots ,{a_m}) f(x)=f(x,a1,a2,,am) 就是关于参数 a 1 , a 2 , ⋯   , a m {a_1},{a_2}, \cdots ,{a_m} a1,a2,,am 的线性函数,代入

J = ∑ i = 1 n ( f ( x i ) − y i ) 2 J = \sum\limits_{i = 1}^n {{{(f({x_i}) - {y_i})}^2}} J=i=1n(f(xi)yi)2

则函数 J = J ( a 1 , a 2 , ⋯   , a k ) J = J({a_1},{a_2}, \cdots ,{a_k}) J=J(a1,a2,,ak) 是关于参数 a 1 , a 2 , ⋯   , a m {a_1},{a_2}, \cdots ,{a_m} a1,a2,,am 的多元函数。由

∂ J ∂ a k = 0 , k = 1 , 2 , ⋯   , m \frac{{\partial J}}{{\partial {a_k}}} = 0,\quad k = 1,2, \cdots ,m akJ=0,k=1,2,,m

亦即

∑ i = 1 n [ ( f ( x i ) − y i ) φ k ( x i ) ] = 0 , k = 1 , 2 , ⋯   , m \sum\limits_{i = 1}^n {[(f({x_i}) - {y_i}){\varphi _k}({x_i})]} = 0,\quad k = 1,2, \cdots ,m i=1n[(f(xi)yi)φk(xi)]=0,k=1,2,,m

可得

∑ j = 1 m [ ∑ i = 1 n φ j ( x i ) φ k ( x i ) ] a j = ∑ i = 1 n y i φ k ( x i ) , k = 1 , 2 , ⋯   , m \sum\limits_{j = 1}^m {\left[ {\sum\limits_{i = 1}^n {{\varphi _j}({x_i}){\varphi _k}({x_i})} } \right]{a_j}} = \sum\limits_{i = 1}^n {{y_i}{\varphi _k}({x_i})} ,\quad k = 1,2, \cdots ,m j=1m[i=1nφj(xi)φk(xi)]aj=i=1nyiφk(xi),k=1,2,,m

R = ( φ 1 ( x 1 ) φ 2 ( x 1 ) ⋯ φ m ( x 1 ) φ 1 ( x 2 ) φ 2 ( x 2 ) ⋯ φ m ( x 2 ) ⋮ ⋮ ⋮ φ 1 ( x n ) φ 2 ( x n ) ⋯ φ m ( x n ) ) {\mathbf{R}} = \left( {\begin{array}{l} {{\varphi _1}({x_1})}&{{\varphi _2}({x_1})}& \cdots &{{\varphi _m}({x_1})} \\ {{\varphi _1}({x_2})}&{{\varphi _2}({x_2})}& \cdots &{{\varphi _m}({x_2})} \\ \vdots & \vdots &{}& \vdots \\ {{\varphi _1}({x_n})}&{{\varphi _2}({x_n})}& \cdots &{{\varphi _m}({x_n})} \end{array}} \right) R=φ1(x1)φ1(x2)φ1(xn)φ2(x1)φ2(x2)φ2(xn)φm(x1)φm(x2)φm(xn)

A = ( a 1 a 2 ⋮ a m ) {\mathbf{A}} = \left( {\begin{array}{c} a_1 \\ {{a_2}} \\ \vdots \\ {{a_m}} \end{array}} \right) A=a1a2am

Y = ( y 1 y 2 ⋮ y n ) {\mathbf{Y}} = \left( {\begin{array}{c} {{y_1}} \\ {{y_2}} \\ \vdots \\ {{y_n}} \end{array}} \right) Y=y1y2yn

则正规方程组可表示为
R T R A = R T Y {{\mathbf{R}}^T}{\mathbf{RA}} = {{\mathbf{R}}^T}{\mathbf{Y}} RTRA=RTY

由代数知识可知,当矩阵 R {\mathbf{R}} R 是列满秩时 R T R {{\mathbf{R}}^T}{\mathbf{R}} RTR 是可逆的。于是正规方程组有唯一解,即
A = ( R T R ) − 1 R T Y {\mathbf{A}} = {({{\mathbf{R}}^T}{\mathbf{R}})^{ - 1}}{{\mathbf{R}}^T}{\mathbf{Y}} A=(RTR)1RTY

为所求的拟合函数的系数,就可得到最小二乘拟合函数 f ( x ) f(x) f(x)

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值