§ 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 iii
yxPw
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 iii
yxPw
1
2])([
② 连续型 /*continuous type */
已知 y(x)?C[a,b] 以及权函数?(x),求广义多项式 P(x) 使得误差函数?= 最小 。dxxyxPxb
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
nkya k
n
j
jjk,...,0,),(),(
0
设则完全类似地有：
)(...)()()( 1100 xaxaxaxP nn +++?
0
ka
法方程组
/*normal equations */
定理 Ba = c 存在唯一解0(x),?1(x),…,?n(x) 线性无关。
即：
),(
),(
),(
00
y
y
a
a
b
nn
jiij
= c
证明,若存在一组系数 {?i } 使得 0...1100?+++ nn
则等式两边分别与?0,?1,…,?n作内积，得到：
+++
+++
+++
0),(.,,),(),(
.
.
.
0),(.,,),(),(
0),(.,,),(),(
1100
1111100
0011000
nnnnn
nn
nn
即,B? = 0
… …
§ 2 Orthogonal Polynomials & L-S Approximation
例,用 来拟合,w? 12
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)(463||||484,|||| 1 BcondBB
§ 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? 12
210 xcxccy ++?
x 1 2 3 4 y 4 10 18 26
解,通过正交多项式?0(x),?1(x),?2(x) 求解设 )()()( 221100 xaxaxay ++? ),( ),( kkkk ya
1)(0?x? 229),( ),(
00
0
0
ya
2
5
),(
),(
00
00
1
x
2
5)()()(
011 xxxx 5
37
),(
),(
11
1
1
ya
2
5
),(
),(
11
11
2
x
4
5
),(
),(
00
11
1
55)(45)()25()( 2012 + xxxxxx 21),( ),(
22
2
2
ya
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) = (x1)?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) = (xk)?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+0007.2279e?0021.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 iii
yxPw
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 iii
yxPw
1
2])([
② 连续型 /*continuous type */
已知 y(x)?C[a,b] 以及权函数?(x),求广义多项式 P(x) 使得误差函数?= 最小 。dxxyxPxb
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
nkya k
n
j
jjk,...,0,),(),(
0
设则完全类似地有：
)(...)()()( 1100 xaxaxaxP nn +++?
0
ka
法方程组
/*normal equations */
定理 Ba = c 存在唯一解0(x),?1(x),…,?n(x) 线性无关。
即：
),(
),(
),(
00
y
y
a
a
b
nn
jiij
= c
证明,若存在一组系数 {?i } 使得 0...1100?+++ nn
则等式两边分别与?0,?1,…,?n作内积，得到：
+++
+++
+++
0),(.,,),(),(
.
.
.
0),(.,,),(),(
0),(.,,),(),(
1100
1111100
0011000
nnnnn
nn
nn
即,B? = 0
… …
§ 2 Orthogonal Polynomials & L-S Approximation
例,用 来拟合,w? 12
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)(463||||484,|||| 1 BcondBB
§ 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? 12
210 xcxccy ++?
x 1 2 3 4 y 4 10 18 26
解,通过正交多项式?0(x),?1(x),?2(x) 求解设 )()()( 221100 xaxaxay ++? ),( ),( kkkk ya
1)(0?x? 229),( ),(
00
0
0
ya
2
5
),(
),(
00
00
1
x
2
5)()()(
011 xxxx 5
37
),(
),(
11
1
1
ya
2
5
),(
),(
11
11
2
x
4
5
),(
),(
00
11
1
55)(45)()25()( 2012 + xxxxxx 21),( ),(
22
2
2
ya
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) = (x1)?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) = (xk)?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+0007.2279e?0021.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
