1、第六章第六章 曲线拟合与函数逼近曲线拟合与函数逼近/*Approximation Theory*/仍然是已知仍然是已知 x1 xm;y1 ym,求一个简单易求一个简单易算的近似函数算的近似函数 P(x)f(x)。但是但是 m 很大；很大；yi 本身是测量值，不准确，即本身是测量值，不准确，即 yi f(xi)这时没必要取这时没必要取 P(xi)=yi,而要使而要使 P(xi)yi 总体上总体上尽可能小。尽可能小。常见做法：常见做法：使使 最小最小/*minimax problem*/|)(|max1iimiyxP 太复杂太复杂 使使 最小最小 miiiyxP1|)(|不可导，求解困难不可导，求
2、解困难 使使 最小最小 /*Least-Squares method*/miiiyxP12|)(|1 最小二乘拟合最小二乘拟合多项式多项式 /*L-S approximating polynomials*/确定多项式确定多项式 ，对于一组数，对于一组数据据(xi,yi)(i=1,2,n)使得使得 达到达到极小极小，这里这里 n m。nnxaxaaxP .)(10 miiiyxP12)(naaa10 实际上是实际上是 a0,a1,an 的多元函数，即的多元函数，即 miinininyxaxaaaaa121010.),.,(在在 的极值点应有的极值点应有nkak,.,0,0 kimiiikaxPy
3、xPa )()(201 kiminjijijxyxa 102 njmikiimikjijxyxa0112记记 mikiikmikikxycxb11,nnnnnnccaabbbb.000000法方程组法方程组(或或正规方程组正规方程组)/*normal equations*/回归系数回归系数/*regression coefficients*/1 L-S Approximating Polynomials定理定理L-S 拟合多项式拟合多项式存在唯一存在唯一(n 0,b 0)线性化：线性化：由由 可做变换可做变换xbay lnlnbBaAxXyY ,ln,1,lnBXAY 就是个就是个线性问题线性
4、问题将将 化为化为 后易解后易解 A 和和B),(iiYX),(iiyxxbAeaxPBbea/)(,HW:p.233#7，#9,#10,#11 例例 用用 来拟合来拟合 。xbaexp/)(2 正交多项式与最小二乘拟合正交多项式与最小二乘拟合/*Orthogonal Polynomials&Least-Squares Approximation */已知已知 x1 xm;y1 ym,求一个简单易算的近求一个简单易算的近似函数似函数 P(x)f(x)使得使得 最小。最小。miiiyxP12|)(|已知已知 a,b上定义的上定义的 f(x)，求一个简单易算的，求一个简单易算的近似函数近似函数 P
5、(x)使得使得 最小。最小。badxxfxP2)()(定义定义线性无关线性无关/*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),，其有限项的线性组合，其有限项的线性组合 称为称为广义广义多项式多项
6、式/*generalized polynomial*/.njjjxxP0)()(常见多项式：常见多项式：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定义定义权函数：权函数：离散型离散型/*discre
7、te type*/根据一系列离散点根据一系列离散点 拟合时，在每一误拟合时，在每一误差前乘一正数差前乘一正数wi，即，即 误差函数误差函数 ，这个，这个wi 就称作就称作权权/*weight*/，反映该点的重要程度。，反映该点的重要程度。),.,1(),(niyxii niiiiyxPw12)(连续型连续型/*continuous type*/在在a,b上用广义多项式上用广义多项式 P(x)拟合连续函数拟合连续函数 f(x)时，定义时，定义权权函数函数 (x)Ca,b，即误差函数，即误差函数 =。权函数必须权函数必须(x)满足：非负、可积，且在满足：非负、可积，且在a,b的任何子区的任何子区间
8、上间上(x)0。dxxyxPxba2)()()(2 Orthogonal Polynomials&L-S Approximation定义定义广义广义 L-S 拟合：拟合：离散型离散型/*discrete type*/在点集在点集 x1 xm 上测得上测得 y1 ym，在一组权系数，在一组权系数 w1 wm 下求广义多项式下求广义多项式 P(x)使得使得误差函数误差函数 最小。最小。niiiiyxPw12)(连续型连续型/*continuous type*/已知已知 y(x)Ca,b 以及权函数以及权函数 (x)，求广义多项式求广义多项式 P(x)使使得误差函数得误差函数 =最小最小。dxxyx
9、Pxba2)()()(内积内积与与范数范数 bamiiiidxxgxfxxgxfwgf)()()()()(),(1 离散型离散型连续型连续型则易证则易证(f,g)是是内积内积，而而 是是范数范数。),(|fff(f,g)=0 表示表示 f 与与 g 带权正交带权正交。广义广义 L-S 问题可叙述为：求广义多项式问题可叙述为：求广义多项式P(x)使得使得 最小。最小。2|),(yPyPyP 2 Orthogonal Polynomials&L-S Approximationnkyaknjjjk,.,0,),(),(0 设设则完全类似地有：则完全类似地有：)(.)()()(1100 xaxaxax
10、Pnn 0 ka 法方程组法方程组/*normal equations*/定理定理 Ba=c 存在唯一解存在唯一解 0(x),1(x),n(x)线性无关。线性无关。即：即：),(),(),(00yyaabnnjiij =c证明：证明：若存在一组系数若存在一组系数 i 使得使得0.1100 nn 则等式两边分别与则等式两边分别与 0,1,n作内积，得到：作内积，得到：0),(.),(),(.0),(.),(),(0),(.),(),(110011111000011000nnnnnnnnn 即：即：B =0 2 Orthogonal Polynomials&L-S Approximation例：例
11、：用用 来拟合来拟合 ，w 12210 xaxaay x 1 2 3 4 y 4 10 18 26 解：解：0(x)=1，1(x)=x，2(x)=x2622),(182),(581),(354),(301),(30),(101),(100),(411),(241104144122220412411110412412100 yyyyxxxxxxiiiiiiiiiiiiii 6221825835410030100301030104210aaa21,1049,23210 aaa23104921)(2 xxxPyIt is soooo simple!What can possibly go wrong
12、?7623)(463|484,|1 BcondBB2 Orthogonal Polynomials&L-S Approximation例：例：连续型拟合中，取连续型拟合中，取1,0)(,1)(,)(Cxyxxxjj 则则 1011),(jidxxxjiji Hilbert阵！阵！改进：改进：若能取函数族若能取函数族=0(x),1(x),n(x),，使得任意一对使得任意一对 i(x)和和 j(x)两两两两（带权）正交（带权）正交，则则 B 就化为就化为对角阵对角阵！这时直接可算出这时直接可算出ak=),(),(kkky Well,no free lunch anyway 正交正交多项式多项式的构
13、造：的构造：将正交函数族中的将正交函数族中的 k 取为取为k 阶阶多项式多项式，为简单起见，可取，为简单起见，可取 k 的的首项系数为首项系数为 1。有递推有递推关系式：关系式：)()()(,1)(0110 xxxx )()()()(111xxxxkkkkk 其中其中),(),(,),(),(111 kkkkkkkkkkx 证明略证明略2 Orthogonal Polynomials&L-S Approximation例：例：用用 来拟合来拟合 ，w 12210 xcxccy x 1 2 3 4 y 4 10 18 26 解：解：通过正交多项式通过正交多项式 0(x),1(x),2(x)求解求
14、解设设)()()(221100 xaxaxay ),(),(kkkkya 1)(0 x 229),(),(0000 ya25),(),(00001 x25)()()(011 xxxx 537),(),(1111 ya25),(),(11112 x45),(),(00111 55)(45)()25()(2012 xxxxxx 21),(),(2222 ya23104921)55(21)25(537122922 xxxxxy与前例结果一致。与前例结果一致。注：注：手算时也可手算时也可用待定系数法确用待定系数法确定函数族。定函数族。2 Orthogonal Polynomials&L-S Appro
15、ximation Algorithm:Orthogonal Polynomials Approximation To approximate a given function by a polynomial with error boundedby a given tolerance.Input:number of data m;xm;ym;weight wm;tolerance TOL;maximum degree of polynomial Max_n.Output:coefficients of the approximating polynomial.Step 1 Set 0(x)1;
16、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 0 distinct points .You are supposed to write a functionvoid OPA(double(*f)(),double x,double w,int m,double tol,FILE*outfile)to approximate t
17、he 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 total error must be no larger than tol.mxxx .21 miiniixPxfxwerr12)()()(2 Orthogonal Polynomials&L-S ApproximationInputTher
18、e 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
19、the degree of the polynomial in the format:fprintf(outfile,%dn,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(o
20、utfile,error=%12.8en,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.,10nnnxaxaaxP .)(102 Orthogonal Polynomials&L-S ApproximationSample Judge ProgramSample Judge Program#i
21、nclude#include#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 xMAX_m,wMAX_m,tol;m=90;for(i=0;im;i+)xi=3.1415926535897932;xi=xi*(double)(i+1)/180.0;wi=1.0;tol=0.001;OPA(f1,
22、x,w,m,tol,outfile);m=200;for(i=0;im;i+)xi=0.01*(double)i;wi=1.0;tol=0.001;OPA(f2,x,w,m,tol,outfile);fclose(outfile);2 Orthogonal Polynomials&L-S ApproximationSample OutputSample Output(represents a space)represents a space)3 2.5301e 0031.0287e+000 7.2279e 002 1.1287e 001error=6.33097847e 005 41.0025
23、e+0009.6180e 0016.2900e 0017.0907e 0031.1792e 001error=1.61711536e 0042 函数的最佳逼近函数的最佳逼近 /*Optimal Approximation*/最佳最佳平方平方逼近：即连续型逼近：即连续型L-S逼近，在逼近，在 意义意义 下，使得下，使得 最小。最小。),(|2fff 2|yP 最佳最佳一致一致逼近逼近/*uniform approximation*/|)(|max|,xffbax 在在 意义下，使得意义下，使得 最小。也称最小。也称为为minimax problem。|yP偏差偏差/*deviation*/若若
24、，则称，则称 x0 为为 偏差点偏差点。|)()(00yPxyxP Didnt you say its a very difficult problem?Take it easy.Its not so difficult if we consider polynomials only.3 Optimal Approximationv 1.0 最佳一致逼近多项式最佳一致逼近多项式/*optimal uniform approximating polynomial*/的构造：求的构造：求 n 阶多项式阶多项式 Pn(x)使得使得|Pn y|最最小。小。直接构造直接构造 OUAP 的确比较困难，不妨
25、换个角度，先的确比较困难，不妨换个角度，先考察它应该具备的考察它应该具备的性质性质。有如下结论：。有如下结论：OUAP 存在，且必同时有存在，且必同时有 偏差点。偏差点。证明：证明：存在性证明略。后者用反证法，设只有正偏差点。存在性证明略。后者用反证法，设只有正偏差点。设设nnbaxnExyxPyP|)()(|max|,而对于所有的而对于所有的 x a,b 都有都有nnExyxP )()(nnnExyxPE )()(2/|)(2/)(|nnExyxP是是n阶多项式阶多项式是误差更小的多项式是误差更小的多项式3 Optimal Approximation（Chebyshev定理）定理）Pn 是是
26、 y 的的OUAP Pn 关于关于 y 在定义域在定义域上至少有上至少有n+2个个交错交错的的 偏差点。偏差点。即存在点集即存在点集 a t1 tn+2 b 使得使得 tk 称为称为切比雪夫交错组切比雪夫交错组/*Chebyshev alternating sequence*/|)1()()(yPtytPnkkkn 若若 且且 y 不是不是 n 次多项式，则次多项式，则 n 次次OUAP 唯一唯一。,baCy 证明：证明：反证，设有反证，设有2个个OUAPs，分别是，分别是Pn 和和 Qn 。则它们的平均函数则它们的平均函数 也是一个也是一个OUAP。2)()()(xQxPxRnnn 对于对于
27、Rn 有有Chebyshev交错组交错组 t1,tn+2 使得使得nkknkknkknnEtytQtytPtytRE|)()(|21|)()(|21|)()(|nkknkknEtytQtytP|)()(|)()(|则至少在一个点上必须有则至少在一个点上必须有)()()()(knkkkntQtytytP 0)()(kkntytR0 nE3 Optimal Approximation 由由Chebyshev定理可推出：定理可推出：Pn(x)y(x)在定义域上至少变号在定义域上至少变号 次，故至少有次，故至少有 个个根根。xy0yy x()yy xEn ()yy xEn ()yP xn()n+1n+
28、1可见可见Pn(x)是是 y(x)的的某一个某一个插值多项式插值多项式 如何确定如何确定插值节点插值节点 x0,xn 的位置，使得的位置，使得Pn(x)刚好是刚好是 y 的的OUAP？即，使插值余项即，使插值余项v 2.0达到极小？达到极小？niinnxxnyxR0)1()()!1()(|)(|3 Optimal Approximationv 2.1 在在 1,1上求上求 x1,xn 使得使得 的的|wn|最小。最小。niinxxxw1)()(注意到注意到 ，要使，要使|wn|最小就意味着最小就意味着)()(1xPxxwnnn v 3.0 在在 1,1上求函数上求函数 xn 的的n 1阶阶 O
29、UAP。由由Chebyshev定理可推出：定理可推出：Pn 1(x)关于关于xn 有有n+1个偏个偏差点，即差点，即wn(x)在在 n+1个点上交错取极大、极小值。个点上交错取极大、极小值。v 3.1 在在 1,1上求上求切比雪夫交错组切比雪夫交错组 t1,tn+1 。切比雪夫多项式切比雪夫多项式/*Chebyshev polynomials*/3 Optimal Approximation考虑三角函数考虑三角函数 cos(n )在在 0,上的上的 个极值点。个极值点。n+1当当 时，时，cos(n )交错达到极大值交错达到极大值 1 和极和极小值小值 1，且存在系数，且存在系数 a0,an
30、使得使得),.,1,0(nknkk nkkkan0)(cos)cos(令令 x=cos()，则，则 x 1,1。)cos arccos()cos()(xnnxTn 称为称为Chebyshev多项式多项式 Tn 的重要性质：的重要性质：当当 时，时，交错取到极大值交错取到极大值 1 和极小值和极小值 1，即，即),.,1,0(cosnknktk )(kntT|)(|)1()(xTtTnkkn1 当当 时时 ，即，即 x1,xn 为为Tn(x)的的n个零点。个零点。),.,1(212cosnknkxk 0)(knxT3 Optimal Approximation Tn(x)满足递推关系：满足递推关
31、系：T0(x)=1，T1(x)=x，)()(2)(11xTxTxxTnnn Tn(x)为为 n 次多项式，首项系数为次多项式，首项系数为 。且。且T2n(x)只含只含 x 的的 次幂，次幂，T2n+1(x)只含只含x 的的 次幂。次幂。2n 1偶偶奇奇 T0(x),T1(x),是是 1,1 上关于权上关于权 正交正交的函数族。即在内积的函数族。即在内积 的意义下有的意义下有211)(xx 11)()()(),(dxxTxTxTTlklk 0200),(lklklkTTlk OKOK,I think its enough for us Whats our target again?v 3.1 在
32、在 1,1上求上求切比雪夫交错组切比雪夫交错组 t1,tn+1 。v 3.0 在在 1,1上求函数上求函数 xn 的的n 1阶阶 OUAP。Tn(x)的的n个个零点零点。3 Optimal Approximation 可见：若取可见：若取 ，则，则wn在在 1,1 上有上有 n+1 个个极极值点值点 tk，也即，也即Pn 1(x)=xn wn(x)关于关于xn在在 1,1 上有上有n+1个交错个交错偏差点偏差点 tk 。12)()(nnnxTxwv3.0 OKv 2.1 在在 1,1上求上求 x1,xn 使得使得 的的|wn|最小。最小。niinxxxw1)()(取最小值取最小值|)(|1nn
33、xxP )(21|min1xTwnnnwnn121 n n=首项系数为首项系数为1的的 n 阶多项式阶多项式/*monic polynomials of degree n*/x1,xn 即为即为 如何确定插值节点如何确定插值节点 x0,xn 的位置，使得的位置，使得Pn(x)刚刚好是好是 y 的的OUAP？即，使插值余项？即，使插值余项达到极小？达到极小？v 2.0 niinnxxnyxR0)1()()!1()(|)(|取取 x0,xn 为为Tn+1(x)的的n+1个个零点零点，做，做 y 的插值多项式的插值多项式Pn(x)，则插值余项的上界可达极小，则插值余项的上界可达极小 。)!1(2 n
34、Mn3 Optimal Approximation注：注：上界上界最小不表示最小不表示|Rn(x)|最小，故最小，故Pn(x)严格意义上只是严格意义上只是y(x)的的近似近似最佳逼近多项式；最佳逼近多项式；对于一般区间对于一般区间 x a,b，可作变量替换，可作变量替换 ，则，则 t 1,1 ，这时，这时tabbax22 ).()()(220222211nabbaabbaabbannxtxttwxw ).(012nnabtttt )()(12)(12112121tTtTnabnnabnnn 即以即以 为插值节点为插值节点(k=0,n)，得得Pn(x)，余项，余项 有最小上界。有最小上界。221
35、2cos22nkabbaxk)(2)()!1()()(1121)1(tTabnyxRnnnnn 3 Optimal Approximation例：例：求求 f(x)=ex 在在0,1上的近似最佳逼近多项式，使其误差上的近似最佳逼近多项式，使其误差不超过不超过 0.5 10 4。解：解：根据误差上界确定根据误差上界确定 n:4102121)!1(|12 nnneRn=4 计算计算 T5(t)的根：的根：109cos,107cos,105cos,103cos,10cos43210 ttttt)1(2122 ttabbax02.0)1109(cos2121.0)1107(cos21,50.0)110
36、5(cos2179.0)1103(cos21,98.0)110(cos2143210 xxxxx 以以 x0,x4 为节点作为节点作L4(x)3 Optimal Approximation Chebyshev 多项式的其它应用多项式的其它应用 多项式降次多项式降次/*reduce the degree of polynomial with a minimal loss of accuracy*/设设 f(x)Pn(x)。在降低在降低 Pn(x)次数的同时次数的同时,使使因此增加的误差尽可能小因此增加的误差尽可能小,也叫也叫 economiza-tion of power series。从从 P
37、n中去掉一个含有其最高次项的中去掉一个含有其最高次项的 ,结果降结果降次为次为 ,则：则：PnPn 1|)(|max|)()(|max|)()(|max1,11,111,1xPxPxfxPxfnnn 因降次而增的误差因降次而增的误差设设 Pn 的首项系数为的首项系数为an，则取，则取 可使可使精度尽可能少损失。精度尽可能少损失。12)()(nnnnxTaxP3 Optimal Approximation例：例：f(x)=ex 在在 1,1上的上的4 阶阶 Taylor 展开为展开为246214324xxxxP ，此时误差，此时误差023.0|!5|)(|54 xexR请将其请将其降为降为2阶多
38、项式阶多项式。解：解：取取)81(241)(2124124434 xxxTP188244 xxT（查表知（查表知 ）)81(24162123244 xxxxPP32612413192191xxx 取取)43(61)(21613323xxxTP xxT3433 （查表知（查表知 ）192191892413233 xxPP057.0|)(|2 xPex若简单取若简单取 ，则误差，则误差21)(22xxxP 45.0!3 e另类解法可查另类解法可查p.163表表7-2，将将x3 和和x4 中的中的T3 和和T4 删除。删除。注：注：对一般区间对一般区间a,b，先将，先将 x 换为换为 t，考虑，考虑 f(t)在在 1,1上上的逼近的逼近Pn(t)，再将，再将 t 换回换回x，最后得到，最后得到Pn(x)。HW:p.164#3
