§ 2 正交多项式与最小二乘拟合
/* Orthogonal Polynomials & Least-Squares Approximation */
已知 x1 … xm ; y1 … ym,求一个简单易算的近
似函数 P(x) ? f(x) 使得 最小。 ?
?
?m
i
ii yxP
1
2|)(|
已知 [a,b]上定义的 f(x)，求一个简单易算的
近似函数 P(x) 使得 最小。 ? ?b
a dxxfxP
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。
§ 2 Orthogonal Polynomials & L-S Approximation
定义 考虑一般的线性无关函数族 ?={ ?0(x),?1(x),…,
?n(x),… },其有限项的线性组合 称为 广义
多项式 /* generalized polynomial */,
?
?
?
n
j
jj xxP
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 */
§ 2 Orthogonal Polynomials & L-S Approximation
定义 权函数,
① 离散型 /*discrete type */
根据一系列离散点 拟合时，在每一误
差前乘一正数 wi,即 误差函数 ?,这个 wi
就称作 权 /* weight*/，反映该点的重要程度。
),.,,,1(),( niyx ii ?
?
?
? ?
n
i i i i
y x P w
1
2 ] ) ( [
② 连续型 /*continuous type */
在 [a,b]上用广义多项式 P(x) 拟合连续函数 f(x) 时，定义 权
函数 ?(x) ?C[a,b]，即误差函数 ? = 。
权函数必须 ?(x)满足：非负、可积，且在 [a,b]的任何子区
间上 ?(x) ? 0。
dxxyxPxba 2)]()([)( ?? ?
§ 2 Orthogonal Polynomials & L-S Approximation
定义 广义 L-S 拟合,
① 离散型 /*discrete type */
在点集 { x1 … xm } 上测得 { y1 … ym }，在一组权系数 { w1 …
wm }下求广义多项式 P(x) 使得 误差函数 ?
最小。
?
?
? ?
n
i i i i
y x P w
1
2 ] ) ( [
② 连续型 /*continuous type */
已知 y(x) ?C[a,b] 以及权函数 ?(x),求广义多项式 P(x) 使
得误差函数 ? = 最小 。 dx x y x P x b
a
2 )] ( ) ( [ ) ( ? ? ?
内积 与 范数
?
?
?
?
?
?
?
?
?
b
a
m
i
iii
dxxgxfx
xgxfw
gf
)()()(
)()(
),( 1
?
离散型
连续型 则易证 ( f,g ) 是 内积,而 是 范数 。
),(|||| fff ?
( f,g )=0 表示 f 与 g
带权正交 。
广义 L-S 问题可叙述为：求广义多项式 P(x)使得
最小。 2||||),( yPyPyP ??????
§ 2 Orthogonal Polynomials & L-S Approximation
n k y a k
n
j
j j k,...,0,),( ),(
0
? ? ?
?
? ? ?
设
则完全类似地有,
) (,.,) ( ) ( ) ( 1 1 0 0 x a x a x a x P n n ? ? ? + + + ?
0???
ka
?
法方程组
/*normal equations */
定理 Ba = c 存在唯一解 ? ?0(x),?1(x),…,?n(x) 线性无关。
即,
),(
),(
),(
0 0
y
y
a
a
b
n n
j i ij
?
?
? ? ? ? = c
证明,若存在一组系数 {?i } 使得 0,.,1 1 0 0 ? + + + n n ? ? ? ? ? ?
则等式两边分别与 ?0,?1,…,?n作内积，得到,
?
?
?
?
?
?
?
?+++
?+++
?+++
0),(.,,),(),(
.
.
.
0),(.,,),(),(
0),(.,,),(),(
1100
1111100
0011000
nnnnn
nn
nn
?????????
?????????
?????????
即,B ? = 0
… …
§ 2 Orthogonal Polynomials & L-S Approximation
例,用 来拟合, w ? 1 2
210 xaxaay ++?
x 1 2 3 4 y 4 10 18 26
解,?0(x) = 1,?1(x) = x,?2(x) = x2
6 2 2),(1 8 2),(581),(
3 5 4),(301),(
30),(101),(
1 0 0),(411),(
2
4
1
10
4
1
4
4
1
22
2
20
4
1
2
4
1
1110
4
1
2
4
1
2100
?????
?????
?????
??????
?
??
??
??
?
??
??
??
yyyy
xx
xx
xx
i
i
i
i
i
i
i
i
i
i
i
ii
i
???
????
????
????
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
622
182
58
35410030
1003010
30104
2
1
0
a
a
a 2
1,
10
49,
2
3
210 ???? aaa
2
3
10
49
2
1)( 2 ?+?? xxxPy
It is soooo simple! What
can possibly go wrong?
7623 ) ( 4 63 || || 484,|| || 1 ? ? ? ? ? ? B cond B B
§ 2 Orthogonal Polynomials & L-S Approximation
例,连续型拟合中，取 ]1,0[)(,1)(,)( Cxyxxx j
j ??? ??
则 ?
++??
1
0 1
1),(
jidxxx
ji
ji ??
Hilbert阵！
改进,
若能取函数族 ?={ ?0(x),?1(x),…,?n(x),… }，
使得任意一对 ?i(x)和 ?j(x)两两 （带权）正交,
则 B 就化为 对角阵 ！
这时直接可算出 ak =
),(
),(
kk
k y
??
?
Well,no free lunch
anyway…
? 正交 多项式 的构造,
将正交函数族中的 ?k 取为 k 阶 多项式,为简单起见，可取
?k 的 首项系数为 1 。
有递推
关系式,
)()()(,1)( 0110 xxxx ???? ???
)()()()( 111 xxxx kkkkk ?++ ??? ?????
其中
),(
),(,
),(
),(
11
1
??
+ ??
kk
kkk
kk
kkk x
??
???
??
???
证明略
p.148-149
§ 2 Orthogonal Polynomials & L-S Approximation
例,用 来拟合, w ? 1 2
210 xcxccy ++?
x 1 2 3 4 y 4 10 18 26
解,通过正交多项式 ?0(x),?1(x),?2(x) 求解
设 ) ( ) ( ) ( 2 2 1 1 0 0 x a x a x a y ? ? ? + + ? ),( ),( kkkk ya ????
1 ) ( 0 ? x ? 2 29 ),( ),(
0 0
0
0 ? ? ? ?
? y a
2
5
),(
),(
0 0
0 0
1 ? ? ? ?
? ? ? x
2
5 ) ( ) ( ) (
0 1 1 ? ? ? ? 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 0 1 2 + ? ? ? ? ? x x x x x x ? ? ? 2 1 ),( ),(
2 2
2
2 ? ? ? ?
? y a
2
3
10
49
2
1
)55(
2
1
)
2
5
(
5
37
1
2
29
2
2
?+?
+?+?+??
xx
xxxy
与前例结果一致。
注,手算时也可
用待定系数法确
定函数族。
§ 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|||| yPer r ?? ? ?
? ?
?????? n
k
n
i
iikk yayayPyP
0 0
),(),( ??
??
??
+?? n
k
kk
n
k
kkk yyyaa
00
2 ),(),(2),( ??? ?
?
?? n
k
kk yayy
0
),(),( ?
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,
mxxx ???,..21
? ? ?? mi inii xPxfxwe r r 1 2)]()([)(
§ 2 Orthogonal Polynomials & L-S Approximation
Input
There is no input file,Instead,you must hand in your function in a *.h file,
The rule of naming the *.h file is the same as that of naming the *.c or
*.cpp files,
Output (? represents a space)
For each test case,you are supposed to output the following information,
? The 1st line contains the integer 6 ? n >0 which is the degree of the
polynomial in the format,fprintf(outfile,"%d\n",n );
? The 2nd line contains the n+1 coefficients of the approximation
polynomial where, Each of the
coefficient is to be printed as in C printf,
fprintf(outfile,"%8.4e?",coefficient );
? The 3rd line contains the total error in the format,
fprintf(outfile,"error?=?%12.8e\n",err );
Note,If the total error is still not small enough when n = 6,simply output
the result obtained when n = 6,
The outputs of two test cases must be seperated by a blank line,
naaa,.,,,,10 nnn xaxaaxP +++?,..)( 10
§ 2 Orthogonal Polynomials & L-S Approximation
Sample Judge Program
#include <stdio.h>
#include <math.h>
#define MAX_m 200
#define MAX_n 6
#include "98115001_12.h"
double f1 ( double x )
{ return sin(x); }
double f2 ( double x )
{ return exp(x); }
void main( )
{
FILE *outfile = fopen("out.txt","w");
int m,i;
double x[MAX_m],w[MAX_m],tol;
m = 90;
for (i=0; i<m; i++) {
x[i] = 3.1415926535897932;
x[i] = x[i]* (double)(i+1)/180.0;
w[i] = 1.0;
}
tol = 0.001;
OPA(f1,x,w,m,tol,outfile);
m = 200;
for (i=0; i<m; i++) {
x[i] = 0.01*(double)i;
w[i] = 1.0;
}
tol = 0.001;
OPA(f2,x,w,m,tol,outfile);
fclose(outfile);
}
§ 2 Orthogonal Polynomials & L-S Approximation
Sample Output (? represents a space)
3
?2.5301e?003?1.0287e+000??7.2279e?002??1.1287e?001?
error?=?6.33097847e?005
4
1.0025e+000?9.6180e?001?6.2900e?001?7.0907e?003?1.1792e?001?
error?=?1.61711536e?004
/* Orthogonal Polynomials & Least-Squares Approximation */
已知 x1 … xm ; y1 … ym,求一个简单易算的近
似函数 P(x) ? f(x) 使得 最小。 ?
?
?m
i
ii yxP
1
2|)(|
已知 [a,b]上定义的 f(x)，求一个简单易算的
近似函数 P(x) 使得 最小。 ? ?b
a dxxfxP
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。
§ 2 Orthogonal Polynomials & L-S Approximation
定义 考虑一般的线性无关函数族 ?={ ?0(x),?1(x),…,
?n(x),… },其有限项的线性组合 称为 广义
多项式 /* generalized polynomial */,
?
?
?
n
j
jj xxP
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 */
§ 2 Orthogonal Polynomials & L-S Approximation
定义 权函数,
① 离散型 /*discrete type */
根据一系列离散点 拟合时，在每一误
差前乘一正数 wi,即 误差函数 ?,这个 wi
就称作 权 /* weight*/，反映该点的重要程度。
),.,,,1(),( niyx ii ?
?
?
? ?
n
i i i i
y x P w
1
2 ] ) ( [
② 连续型 /*continuous type */
在 [a,b]上用广义多项式 P(x) 拟合连续函数 f(x) 时，定义 权
函数 ?(x) ?C[a,b]，即误差函数 ? = 。
权函数必须 ?(x)满足：非负、可积，且在 [a,b]的任何子区
间上 ?(x) ? 0。
dxxyxPxba 2)]()([)( ?? ?
§ 2 Orthogonal Polynomials & L-S Approximation
定义 广义 L-S 拟合,
① 离散型 /*discrete type */
在点集 { x1 … xm } 上测得 { y1 … ym }，在一组权系数 { w1 …
wm }下求广义多项式 P(x) 使得 误差函数 ?
最小。
?
?
? ?
n
i i i i
y x P w
1
2 ] ) ( [
② 连续型 /*continuous type */
已知 y(x) ?C[a,b] 以及权函数 ?(x),求广义多项式 P(x) 使
得误差函数 ? = 最小 。 dx x y x P x b
a
2 )] ( ) ( [ ) ( ? ? ?
内积 与 范数
?
?
?
?
?
?
?
?
?
b
a
m
i
iii
dxxgxfx
xgxfw
gf
)()()(
)()(
),( 1
?
离散型
连续型 则易证 ( f,g ) 是 内积,而 是 范数 。
),(|||| fff ?
( f,g )=0 表示 f 与 g
带权正交 。
广义 L-S 问题可叙述为：求广义多项式 P(x)使得
最小。 2||||),( yPyPyP ??????
§ 2 Orthogonal Polynomials & L-S Approximation
n k y a k
n
j
j j k,...,0,),( ),(
0
? ? ?
?
? ? ?
设
则完全类似地有,
) (,.,) ( ) ( ) ( 1 1 0 0 x a x a x a x P n n ? ? ? + + + ?
0???
ka
?
法方程组
/*normal equations */
定理 Ba = c 存在唯一解 ? ?0(x),?1(x),…,?n(x) 线性无关。
即,
),(
),(
),(
0 0
y
y
a
a
b
n n
j i ij
?
?
? ? ? ? = c
证明,若存在一组系数 {?i } 使得 0,.,1 1 0 0 ? + + + n n ? ? ? ? ? ?
则等式两边分别与 ?0,?1,…,?n作内积，得到,
?
?
?
?
?
?
?
?+++
?+++
?+++
0),(.,,),(),(
.
.
.
0),(.,,),(),(
0),(.,,),(),(
1100
1111100
0011000
nnnnn
nn
nn
?????????
?????????
?????????
即,B ? = 0
… …
§ 2 Orthogonal Polynomials & L-S Approximation
例,用 来拟合, w ? 1 2
210 xaxaay ++?
x 1 2 3 4 y 4 10 18 26
解,?0(x) = 1,?1(x) = x,?2(x) = x2
6 2 2),(1 8 2),(581),(
3 5 4),(301),(
30),(101),(
1 0 0),(411),(
2
4
1
10
4
1
4
4
1
22
2
20
4
1
2
4
1
1110
4
1
2
4
1
2100
?????
?????
?????
??????
?
??
??
??
?
??
??
??
yyyy
xx
xx
xx
i
i
i
i
i
i
i
i
i
i
i
ii
i
???
????
????
????
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
622
182
58
35410030
1003010
30104
2
1
0
a
a
a 2
1,
10
49,
2
3
210 ???? aaa
2
3
10
49
2
1)( 2 ?+?? xxxPy
It is soooo simple! What
can possibly go wrong?
7623 ) ( 4 63 || || 484,|| || 1 ? ? ? ? ? ? B cond B B
§ 2 Orthogonal Polynomials & L-S Approximation
例,连续型拟合中，取 ]1,0[)(,1)(,)( Cxyxxx j
j ??? ??
则 ?
++??
1
0 1
1),(
jidxxx
ji
ji ??
Hilbert阵！
改进,
若能取函数族 ?={ ?0(x),?1(x),…,?n(x),… }，
使得任意一对 ?i(x)和 ?j(x)两两 （带权）正交,
则 B 就化为 对角阵 ！
这时直接可算出 ak =
),(
),(
kk
k y
??
?
Well,no free lunch
anyway…
? 正交 多项式 的构造,
将正交函数族中的 ?k 取为 k 阶 多项式,为简单起见，可取
?k 的 首项系数为 1 。
有递推
关系式,
)()()(,1)( 0110 xxxx ???? ???
)()()()( 111 xxxx kkkkk ?++ ??? ?????
其中
),(
),(,
),(
),(
11
1
??
+ ??
kk
kkk
kk
kkk x
??
???
??
???
证明略
p.148-149
§ 2 Orthogonal Polynomials & L-S Approximation
例,用 来拟合, w ? 1 2
210 xcxccy ++?
x 1 2 3 4 y 4 10 18 26
解,通过正交多项式 ?0(x),?1(x),?2(x) 求解
设 ) ( ) ( ) ( 2 2 1 1 0 0 x a x a x a y ? ? ? + + ? ),( ),( kkkk ya ????
1 ) ( 0 ? x ? 2 29 ),( ),(
0 0
0
0 ? ? ? ?
? y a
2
5
),(
),(
0 0
0 0
1 ? ? ? ?
? ? ? x
2
5 ) ( ) ( ) (
0 1 1 ? ? ? ? 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 0 1 2 + ? ? ? ? ? x x x x x x ? ? ? 2 1 ),( ),(
2 2
2
2 ? ? ? ?
? y a
2
3
10
49
2
1
)55(
2
1
)
2
5
(
5
37
1
2
29
2
2
?+?
+?+?+??
xx
xxxy
与前例结果一致。
注,手算时也可
用待定系数法确
定函数族。
§ 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|||| yPer r ?? ? ?
? ?
?????? n
k
n
i
iikk yayayPyP
0 0
),(),( ??
??
??
+?? n
k
kk
n
k
kkk yyyaa
00
2 ),(),(2),( ??? ?
?
?? n
k
kk yayy
0
),(),( ?
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,
mxxx ???,..21
? ? ?? mi inii xPxfxwe r r 1 2)]()([)(
§ 2 Orthogonal Polynomials & L-S Approximation
Input
There is no input file,Instead,you must hand in your function in a *.h file,
The rule of naming the *.h file is the same as that of naming the *.c or
*.cpp files,
Output (? represents a space)
For each test case,you are supposed to output the following information,
? The 1st line contains the integer 6 ? n >0 which is the degree of the
polynomial in the format,fprintf(outfile,"%d\n",n );
? The 2nd line contains the n+1 coefficients of the approximation
polynomial where, Each of the
coefficient is to be printed as in C printf,
fprintf(outfile,"%8.4e?",coefficient );
? The 3rd line contains the total error in the format,
fprintf(outfile,"error?=?%12.8e\n",err );
Note,If the total error is still not small enough when n = 6,simply output
the result obtained when n = 6,
The outputs of two test cases must be seperated by a blank line,
naaa,.,,,,10 nnn xaxaaxP +++?,..)( 10
§ 2 Orthogonal Polynomials & L-S Approximation
Sample Judge Program
#include <stdio.h>
#include <math.h>
#define MAX_m 200
#define MAX_n 6
#include "98115001_12.h"
double f1 ( double x )
{ return sin(x); }
double f2 ( double x )
{ return exp(x); }
void main( )
{
FILE *outfile = fopen("out.txt","w");
int m,i;
double x[MAX_m],w[MAX_m],tol;
m = 90;
for (i=0; i<m; i++) {
x[i] = 3.1415926535897932;
x[i] = x[i]* (double)(i+1)/180.0;
w[i] = 1.0;
}
tol = 0.001;
OPA(f1,x,w,m,tol,outfile);
m = 200;
for (i=0; i<m; i++) {
x[i] = 0.01*(double)i;
w[i] = 1.0;
}
tol = 0.001;
OPA(f2,x,w,m,tol,outfile);
fclose(outfile);
}
§ 2 Orthogonal Polynomials & L-S Approximation
Sample Output (? represents a space)
3
?2.5301e?003?1.0287e+000??7.2279e?002??1.1287e?001?
error?=?6.33097847e?005
4
1.0025e+000?9.6180e?001?6.2900e?001?7.0907e?003?1.1792e?001?
error?=?1.61711536e?004
