《数值分析》课程教学资源（PPT讲稿）第七章 曲线拟合与函数逼近（2/3）

§2正交多项式与最小二乘拟合 /*Orthogonal Polynomials Least-Squares Approximation * 已知x1…,xmn;y1…,ym,求一个简单易算的近 似函数P(x)≈fx)使得∑Px)-P最小。 已知[a,b上定义的fx),求一个简单易算的 近似函数P(x)使得Px)-f(x)最小 定义线性无关 nearly independent函数族{a(x) q(x),…,gp(x),…}满足条件:其中任意函数的线性组合 ag(x)+a1g1(x)+…+angn(x)=0对任意x∈Ia,b成立 当且仅当a0=a1=.=an=0

§2 正交多项式与最小二乘拟合 /* Orthogonal Polynomials & Least-Squares Approximation */ 已知 x1 … xm ; y1 … ym, 求一个简单易算的近 似函数 P(x)  f(x) 使得  最小。 = − m i i i P x y 1 2 | ( ) | 已知 [a, b]上定义的 f(x)，求一个简单易算的 近似函数 P(x) 使得  − 最小。 b a P x f x dx 2 [ ( ) ( )] 定义 线性 无 关/* linearly independent */ 函数 族 { 0 (x), 1 (x), … , n (x), … } 满足条件：其中任意函数的线性组合 a00 (x)+a11 (x)+…+ann (x)=0 对任意 x[a, b]成立 当且仅当 a0= a1=… =an =0

8 2 Orthogonal Polynomials L-S Approximation 定义考虑一般的线性无关函数族φ={(ax),q1(x),…, gn(x),…},其有限项的线性组合P(x)=∑a(x)称为广义 j=0 多项式/ generalized polynomial+l 纔常见多项式: >{q(x)=x}对应代数多项式/ algebraic polynomial >{q(x)= cos x}、{v(x)=sinx}→{q(x),v(x)对应三 角多项式/ trigonometric polynomial* >{q(x)=e,k≠k}对应指数多项式/ exponential polynomial */

§2 Orthogonal Polynomials & L-S Approximation 定义 考虑一般的线性无关函数族={ 0 (x), 1 (x), … , n (x), … }，其有限项的线性组合 称为广义 多项式 /* generalized polynomial*/. = = n j P x j j x 0 ( )   ( ) 常见多项式： ➢ { j (x) = x j } 对应代数多项式 /* algebraic polynomial */ ➢ { j (x) = cos jx }、{ j (x) = sin jx }  { j (x), j (x)}对应三 角多项式 /* trigonometric polynomial */ ➢ { j (x) = e kj x , ki  kj } 对应指数多项式 /* exponential polynomial */

8 2 Orthogonal Polynomials L-S Approximation 定义权函数: 离散型/ discrete type 根据一系列离散点(x;,y)(=1,…,n)拟合时,在每一误 差前乘一正数w,即误差函数=2wPx)-y2,这个 就称作权 weight,反映该点的重要程度。 ②连续型/ continuous type 在[a,b上用广义多项式P(x)拟合连续函数f(x)时,定义权 函数(x)∈Ca,b,即误差函数p=Jp(x)P(x)-y(x)d 权函数必须px)满足:非负、可积,且在a,b的任何子区 间上p(x)≠0

§2 Orthogonal Polynomials & L-S Approximation 定义 权函数： ① 离散型 /*discrete type */ 根据一系列离散点 拟合时，在每一误 差前乘一正数wi ，即 误差函数 ，这个wi 就称作权/* weight*/，反映该点的重要程度。 (x , y ) (i 1, ... , n) i i = = = − n i wi P xi yi 1 2 [ ( ) ] ② 连续型 /*continuous type */ 在[a, b]上用广义多项式 P(x) 拟合连续函数 f(x) 时，定义权 函数 (x) C[a, b]，即误差函数 = 。 权函数必须(x)满足：非负、可积，且在[a, b]的任何子区 间上(x)  0。 x P x y x dx b a 2 ( )[ ( ) − ( )]  

8 2 Orthogonal Polynomials L-S Approximation 定义广义LS拟合: ④离散型/ discrete type" 在点集{x1…,xmn}上测得{y1…,ym},在一组权系数{v1 wm}下求广义多项式P(x)使得误差函数p=EwP(x)-y (,g)=0表示∫与g ②连续型P 带权正交。 已知y(x)∈Ca,b以及px),求广义多项式P(x)使 得误差函数-∫p最小 4p内积与范数 ∑ wif(ig() 离散型 则易证(g)是内积,(,8) 而‖∫=√,∫是范数。 p(x)f(x)g(x)dx-连续型 广义LS问题可叙述为:求广义多项式P(x)使得 =(P-y,P-y)=‖P-y2最小

§2 Orthogonal Polynomials & L-S Approximation 定义 广义 L-S 拟合： ① 离散型 /*discrete type */ 在点集{ x1 … xm } 上测得{ y1 … ym }，在一组权系数{ w1 … wm }下求广义多项式 P(x) 使得误差函数  最小。 = = − n i wi P xi yi 1 2 [ ( ) ] ② 连续型 /*continuous type */ 已知 y(x) C[a, b] 以及权函数 (x)，求广义多项式 P(x) 使 得误差函数 = x P x y x dx 最小。 b a 2  ( )[ ( )− ( )] 内积与范数      =  = b a m i i i i x f x g x dx w f x g x f g ( ) ( ) ( ) ( ) ( ) ( , ) 1  离散型 连续型 则易证( f, g ) 是内积， 而 || f || = ( f , f ) 是范数。 ( f, g )=0 表示 f 与 g 带权正交。 广义 L-S 问题可叙述为：求广义多项式P(x)使得  = (P − y, P − y) = || P − y ||2 最小

8 2 Orthogonal Polynomials L-S approximation 设P(x)=a0p0(x)+a1g1(x)+…+an9n(x) 则完全类似地有:如=0→∑(,g)=(91,y),k=0,m,n 即: b=(9,甲) 法方程组 /normal equations"/ n2 定理|m=c在胜一解9(09(,.9(线性无关 证明:若存在一组系数{cx}使得a090+a191+…+ann=0 则等式两边分别与q,q1,…,q作内积,得到: a(9n,9)+a(,%)+…+a(,9)=0即:Ba=0 a0(q0,1)+a1(q1,91)+…+an(qn,q)=0 a0(q0,n)+a1(q1,n)+…+an(qn,n)=0

§2 Orthogonal Polynomials & L-S Approximation a k y k n n j k j j ( , ) ( , ) , 0, ... , 0  = = =    设 则完全类似地有： P( x) = a00 ( x) + a11 ( x) + ... + a n n ( x) = 0   ak  法方程组 /*normal equations */ 定理 Ba = c 存在唯一解  0 (x), 1 (x), … , n (x) 线性无关。 即： ( , ) ( , ) ( , ) 0 0 y y a a b n n ij i j   =   = = c 证明：若存在一组系数 {i } 使得  0  0 +1 1 + ...+ n  n = 0 则等式两边分别与0 , 1 , … , n作内积，得到：        + + + = + + + = + + + = ( , ) ( , ) ... ( , ) 0 . . . ( , ) ( , ) ... ( , ) 0 ( , ) ( , ) ... ( , ) 0 0 0 1 1 0 0 1 1 1 1 1 0 0 0 1 1 0 0 n n n n n n n n n                            即：B  = 0 … …

8 2 Orthogonal Polynomials L-S Approximation 2|3 例:用y=a0+a1x+a2x2来拟合 y4101826,w≡1 解:qx)=1 It is soooo simple! What (Po, po can possibly go wrong? (9,g)=∑1x=10(,9)= =30 (9,2)=∑1x2=30(q2,q2)=∑x=354 ∑1y1=58(q,y)=182(92,y)=622 41030)a 02 1 1030100a1=182 30100354八a 622 y=P(x) x+—x ‖B‖=484,‖B-‖ 63 4cond(B)=7623

§2 Orthogonal Polynomials & L-S Approximation 例：用 来拟合 ，w  1 2 0 1 2 y = a + a x + a x x 1 2 3 4 y 4 10 18 26 解： 0 (x) = 1， 1 (x) = x， 2 (x) = x 2 ( , ) 1 58 ( , ) 182 ( , ) 622 ( , ) 1 30 ( , ) 354 ( , ) 1 10 ( , ) 30 ( , ) 1 1 4 ( , ) 100 2 4 1 0 1 4 1 4 4 1 2 2 2 0 2 4 1 2 4 1 0 1 1 1 4 1 2 4 1 0 0 1 2 =  = = = =  = = = =  = = = =  = =  =        = = = = = = = y y y y x x x x x x i i i i i i i i i i i i i i                          =                     622 182 58 30 100 354 10 30 100 4 10 30 2 1 0 a a a 2 1 , 10 49 , 2 3 a0 = − a1 = a2 = 2 3 10 49 2 1 ( ) 2 y = P x = x + x − It is soooo simple! What can possibly go wrong? ( ) 7623 4 63 || || 484, || || 1  = = − B  = B cond B

8 2 Orthogonal Polynomials L-S approximation 例:连续型拟合中,取q,(x)=x,p(x)≡1,y(x)∈C0,1 则(q,q)=「mxx!de Hilbert阵! i+j+1 若能取函数族={q(x),q1(x) 改进:使得任意一对qx)和q(x)两两(带权)正交 则B就化为对角阵! 这时直接可算出ak= (k, y) (Pk, PPK >正交多项式的构造: 将正交函数族中的列聊为阶多顺式,为简单起见,可取 q的首项系数为1。 anyway 有递推 q0(x)=1,q1(x)=(x-a1)q0(x) Pu()=(x-ak+ok (x)-Bk u(r) 证明略 关系式: 其中ak+1= (xφk,9) B (x, Pu p148-149 (x, P) (q4-1,9k-1)

§2 Orthogonal Polynomials & L-S Approximation 例：连续型拟合中，取 (x) x , (x) 1, y(x) C[0, 1] j  j =    则  + + = = 1 0 1 1 ( , ) i j x x dx i j i  j Hilbert阵！ 改进： 若能取函数族={ 0 (x), 1 (x), … , n (x), … }， 使得任意一对i (x)和j (x)两两（带权）正交， 则 B 就化为对角阵！ 这时直接可算出ak = ( , ) ( , ) k k k y    Well, no free lunch anyway… ➢ 正交多项式的构造： 将正交函数族中的k 取为k 阶多项式，为简单起见，可取 k 的首项系数为 1 。 有递推 关系式： ( ) 1, ( ) ( ) ( ) 0 x = 1 x = x −1 0 x ( ) ( ) ( ) ( ) k+1 x = x −k+1 k x −  k k−1 x 其中 ( , ) ( , ) , ( , ) ( , ) 1 1 1 − − + = = k k k k k k k k k k x           证明略 p.148-149

8 2 Orthogonal Polynomials L-S Approximation 例:用y=C+c1x+c2x2来拟合P410i826,=1 解:通过正交多项式(x),g(x),q2(x)求解 k 设y=a090(x)+a191(x)+a292(x) (or, P) X)= (90,y)_29 0 (90,90)2 (x9,90)_5 5 (90,90)2 91(x)=(x-a1)90(x)=x (91,y)37 2 (91,91)5 (x91,91) (91,91)_5 (91,91)2 (90,90)4 92(x)=(x-=)01(x)-90(x)=x2-5+5 (92,y)_1 (92,92)2 37 y=-×1+(x-。)+(x2-5x+5) 22 注:手算时也可 210x2与前例结果一致。用待定系数法确 x十 定函数族

§2 Orthogonal Polynomials & L-S Approximation 例：用 来拟合 ，w  1 2 0 1 2 y = c + c x + c x x 1 2 3 4 y 4 10 18 26 解：通过正交多项式 0 (x), 1 (x), 2 (x) 求解 设 ( ) ( ) ( ) y = a0 0 x + a1 1 x + a2 2 x ( , ) ( , ) k k k k y a    =  0 ( x) = 1 2 29 ( , ) ( , ) 0 0 0 0 = =    y a 2 5 ( , ) ( , ) 0 0 0 0 1 = =      x 2 5 ( ) ( ) ( ) 1 1 0  x = x −  x = x − 5 37 ( , ) ( , ) 1 1 1 1 = =    y a 2 5 ( , ) ( , ) 1 1 1 1 2 = =      x 4 5 ( , ) ( , ) 0 0 1 1 1 = =      ( ) 5 5 4 5 ) ( ) 2 5 ( ) ( 2 2 1 0  x = x −  x −  x = x − x + 2 1 ( , ) ( , ) 2 2 2 2 = =    y a 2 3 10 49 2 1 ( 5 5) 2 1 ) 2 5 ( 5 37 1 2 29 2 2 = + − =  + − + − + x x y x x x 与前例结果一致。 注：手算时也可 用待定系数法确 定函数族

8 2 Orthogonal Polynomials L-S approximation Algorithm: Orthogonal Polynomials Approximation To approximate a given function by a polynomial with error bounded by a given tolerance. Input: number of data m; x m: y m weight wm; tolerance Tol; maximum degree of polynomial Mar n. Output: coefficients of the approximating polynomial. Step 1 set po(x)=l; 0=(o,y) /(o, o); P(x)=0 Po(r); err=v,y)-ao (o,y); Step 2 Set a(xo, o)/(o, o);,()=(x-a1polx) a1=(q1,y)/(q1,g1);P(x)+=a1q1(x);err-=a1(q1,y); Step 3 Set k=1 Step 4 While((k< Max n)&&errrOr)) do steps 5-7 Step 5 k++; Se6a=(x1,1(q1,91);Bk-1=(91,g1/(qo,90); g2(x)=(x-ap1(x)-Bk-19(x);ak=(q2,y)(q292); P(x)+=ak P2(x); err -=k(2, v); Step 7 Set po(x)=p(x); p,x)=2() Step 8 Output ( STOP. 注:c=‖P-y=(P-y,P-y)=∑a9-y,∑a1-y) k=0 ∑以9-2(9,)+(,y=(,y)-∑a1(9,y k=0 k=0

§2 Orthogonal Polynomials & L-S Approximation Algorithm: Orthogonal Polynomials Approximation To approximate a given function by a polynomial with error bounded by a given tolerance. Input: number of data m; x[m]; y[m]; weight w[m]; tolerance TOL; maximum degree of polynomial Max_n. Output: coefficients of the approximating polynomial. Step 1 Set 0 (x)  1; a0 = (0 , y)/(0 , 0 ); P(x) = a0 0 (x); err = (y, y) − a0 (0 , y); Step 2 Set 1= (x0 , 0 )/(0 , 0 ); 1 (x) = (x − 1 ) 0 (x); a1 = (1 , y)/(1 , 1 ); P(x) += a1 1 (x); err −= a1 (1 , y); Step 3 Set k = 1; Step 4 While (( k < Max_n)&&(|err|TOL)) do steps 5-7 Step 5 k ++; Step 6 k= (x1 , 1 )/(1 , 1 );  k−1 = (1 , 1 )/(0 , 0 ); 2 (x) = (x − k ) 1 (x) −  k−1 0 (x); ak= (2 , y)/(2 , 2 ); P(x) += ak 2 (x); err −= ak (2 , y); Step 7 Set 0 (x) = 1 (x); 1 (x) = 2 (x); Step 8 Output ( ); STOP. 注： 2 err = || P − y ||   = = = − − = − − n k n i k k i i P y P y a y a y 0 0 ( , ) (  ,  )   = = = − + n k k k n k k k k a a y y y 0 0 2 ( , ) 2 ( , ) ( , ) = = − n k k k y y a y 0 ( , ) ( , )

82 Orthogonal Polynomials& L-S Approximation Another von Neumann quote Young man. in mathematics you dont HW understand things, you just get used to p.152#1 them Lab 12 Orthogonal Polynomials Approximation Given a functionf and a set of 2002m>0 distinct points x, <x,<.<x You are supposed to write a function void opa( double( fo, double xl l, double wl l, int m, double tol, FILE*outfile to approximate the function f by an orthogonal polynomial using the exact function values at the given m points x. The array wm contains the values of a weight function at the given points x. The totalerror er=∑w(xf(x)-P(x)2 must be no larger than tol

Another von Neumann quote : Young man, in mathematics you don't understand things, you just get used to them. HW: p.152 #1 §2 Orthogonal Polynomials & L-S Approximation Lab 12. Orthogonal Polynomials Approximation Given a function f and a set of 200  m > 0 distinct points . You are supposed to write a function void OPA ( double (*f)( ), double x[ ], double w[ ], int m, double tol, FILE *outfile ) to approximate the function f by an orthogonal polynomial using the exact function values at the given m points x[ ]. The array w[m] contains the values of a weight function at the given points x[ ]. The total error must be no larger than tol. x  x   xm ... 1 2  = = − m i i xi Pn xi err w x f 1 2 ( )[ ( ) ( )]