南大数值分析课件第六章-曲线拟合与函数逼近

第六章曲线拟合与函数逼近/*ApproximationTheory*/仍然是已知x1…xm

;y1…ym,求一个简单易算的近似函数P(x)

f(x)。但是①

m很大；②

yi本身是测量值，不准确，即yi

f(xi)这时没必要取P(xi)=yi,而要使P(xi)yi总体上尽可能小。常见做法：

使最小/*minimaxproblem*/

太复杂使最小不可导，求解困难使最小/*Least-Squaresmethod*/第六章曲线拟合与函数逼近仍然是已知x1…xm;1§1最小二乘拟合多项式

/*L-Sapproximatingpolynomials*/确定多项式，对于一组数据(xi,yi)(i=1,2,…,n)使得达到极小，这里n

<<

m。naaa10实际上是a0,a1,…,an的多元函数，即[]=-+++=miinininyxaxaaaaa121010...),...,,(j在的极值点应有kiminjijijxyxa==-=10][2-====+njmikiimikjijxyxa0112记====mikiikmikikxycxb11,法方程组(或正规方程组)/*normalequations*/回归系数/*regressioncoefficients*/§1最小二乘拟合多项式/*L-Sapproxim2§1L-SApproximatingPolynomials定理L-S拟合多项式存在唯一

(n<m)。证明：记法方程组为Ba=c.则有其中对任意，必有。若不然，则存在一个使得…即是n

阶多项式的根则B为正定阵，则非奇异，所以法方程组存在唯一解。Waitasecond!Youonlygavemeacriticalpoint,butit’snotnecessarilyaminimumpoint!§1L-SApproximatingPolynomi3§1L-SApproximatingPolynomials定理

Ba=c的解确是的极小点。即：设a

为解，则任意b=(b0

b1…bn)T

对应的多项式必有==njjjxbxF0)(===--=mimiiiiibyxFyxPa1122)(])([])([)(jj证明：==---=-miiimiiiyxPyxFab1212])([])([)()(jj==---+-=miiimiiiiiyxPyxPxPxF1212])([])()()([==--+-=miiiiimiiiyxPxPxFxPxF112])()][()([2)]()([0注：L-Smethod首先要求设定P(x)的形式。若设n=m1，则可取P(x)为过m个点的m1阶插值多项式，这时=0。P(x)不一定是多项式，通常根据经验确定。§1L-SApproximatingPolynomi4例用来拟合。例用5§1L-SApproximatingPolynomials例：xy(xi,yi),i=1,2,…,m方案一：设baxxxPy+=)(求a和b使得最小。=-+=miiiiybaxxba12)(),(jButhey,thesystemofequationsforaandbisnonlinear!Takeiteasy!Wejusthavetolinearizeit…线性化

/*linearization*/：令，则bXaY+就是个线性问题将化为后易解a和b。),(iiYX),(iiyx§1L-SApproximatingPolynomi6例用来拟合。例用7§1L-SApproximatingPolynomials方案二：设xbeaxPy/)(-=(a>0,b>0)线性化：由可做变换xbay-lnlnbBaAxXyY-====,ln,1,lnBXAY+就是个线性问题将化为后易解A和B),(iiYX),(iiyxHW:p.233#7，#9,#10,#11§1L-SApproximatingPolynomi8例用来拟合。例用9§2正交多项式与最小二乘拟合

/*OrthogonalPolynomials&Least-SquaresApproximation*/已知x1…xm

;y1…ym,求一个简单易算的近似函数P(x)

f(x)使得最小。已知[a,b]上定义的f(x)，求一个简单易算的近似函数P(x)使得最小。定义

线性无关/*linearlyindependent*/函数族{0(x),1(x),…,n(x),…}满足条件：其中任意函数的线性组合

a00(x)+a11(x)+…+ann(x)=0对任意x[a,b]成立当且仅当a0=a1=…=an=0。§2正交多项式与最小二乘拟合已知x1…xm;y10§2OrthogonalPolynomials&L-SApproximation定义考虑一般的线性无关函数族={0(x),1(x),…,n(x),…}，其有限项的线性组合称为广义多项式

/*generalizedpolynomial*/.常见多项式：

{j(x)=xj}对应代数多项式/*algebraicpolynomial*/

{j(x)=cosjx}、{j(x)=sinjx}{j(x),j(x)

}对应三角多项式/*trigonometricpolynomial*/

{j(x)=ekjx,ki

kj

}对应指数多项式/*exponentialpolynomial*/§2OrthogonalPolynomials&L11§2OrthogonalPolynomials&L-SApproximation定义权函数：①

离散型/*discretetype*/根据一系列离散点拟合时，在每一误差前乘一正数wi

，即误差函数

，这个wi

就称作权/*weight*/，反映该点的重要程度。=-=niiiiyxPw12])([②

连续型

/*continuoustype*/在[a,b]上用广义多项式P(x)拟合连续函数f(x)时，定义权函数(x)C[a,b]，即误差函数=。权函数必须(x)满足：非负、可积，且在[a,b]的任何子区间上(x)0。§2OrthogonalPolynomials&L12§2OrthogonalPolynomials&L-SApproximation定义广义L-S拟合：①

离散型/*discretetype*/在点集{x1…xm}

上测得{y1…ym}，在一组权系数{w1…wm}下求广义多项式P(x)使得误差函数最小。

=-=niiiiyxPw12])([②

连续型

/*continuoustype*/已知y(x)

C[a,b]以及权函数(x)，求广义多项式P(x)使得误差函数=最小。dxxyxPxba2)]()([)(-r内积与范数离散型连续型则易证(f,g)是内积，而是范数。(f,g)=0表示f与g

带权正交。广义L-S问题可叙述为：求广义多项式P(x)使得最小。§2OrthogonalPolynomials&L13§2OrthogonalPolynomials&L-SApproximationnkyaknjjjk,...,0,),(),(0===jjj设则完全类似地有：)(...)()()(1100xaxaxaxPnnjjj+++=法方程组

/*normalequations*/定理

Ba=c存在唯一解

0(x),1(x),…,n(x)线性无关。即：),(),(),(00yyaabnnjiijjjjj===c证明：若存在一组系数{i

}使得0...1100=+++nnjajaja则等式两边分别与0,1,…,n作内积，得到：即：B=0……§2OrthogonalPolynomials&L14§2OrthogonalPolynomials&L-SApproximation例：用来拟合，w1解：0(x)=1，1(x)=x，2(x)=x2Itissoooosimple!Whatcanpossiblygowrong?7623)(463||||484,||||1==-=BcondBB§2OrthogonalPolynomials&L15§2OrthogonalPolynomials&L-SApproximation例：连续型拟合中，取则Hilbert阵！改进：若能取函数族={0(x),1(x),…,n(x),…}，使得任意一对i(x)和j(x)两两（带权）正交，则B就化为对角阵！这时直接可算出ak=Well,nofreelunchanyway…

正交多项式的构造：将正交函数族中的k取为k阶多项式，为简单起见，可取k的首项系数为1。有递推关系式：其中证明略§2OrthogonalPolynomials&L16§2OrthogonalPolynomials&L-SApproximation例：用来拟合，w1解：通过正交多项式0(x),1(x),2(x)求解设)()()(221100xaxaxayjjj++=1)(0=xj229),(),(0000==jjjya25),(),(00001==jjjjax25)()()(011-=-=xxxxjaj537),(),(1111==jjjya25),(),(11112==jjjjax45),(),(00111==jjjjb55)(45)()25()(2012+-=--=xxxxxxjjj21),(),(2222==jjjya与前例结果一致。注：手算时也可用待定系数法确定函数族。§2OrthogonalPolynomials&L17§2OrthogonalPolynomials&L-SApproximation

Algorithm:OrthogonalPolynomialsApproximation

Toapproximateagivenfunctionbyapolynomialwitherrorboundedbyagiventolerance.Input:numberofdatam;x[m];y[m];weightw[m];toleranceTOL;maximumdegreeofpolynomialMax_n.Output:coefficientsoftheapproximatingpolynomial.Step1Set0(x)

1;a0=(0,y)/(0,0);P(x)=a00(x);err=(y,y)a0(0,y);Step2Set1=

(x0,0)/(0,0);1(x)

=(x1)0(x);

a1=(1,y)/(1,1);P(x)+=a11(x);err

=a1(1,y);Step3Setk=1;Step4While((k<Max_n)&&(|err|TOL))dosteps5-7

Step5k++;

Step6k=

(x1,1)/(1,1);k1=(1,1)/(0,0);

2(x)

=(xk)1(x)k10(x);ak

=(2,y)/(2,2);

P(x)+=ak

2(x);err

=ak

(2,y);

Step7Set0(x)=1(x);1(x)=2(x);Step8Output();STOP.注：§2OrthogonalPolynomials&L18AnothervonNeumannquote:Youngman,inmathematicsyoudon'tunderstandthings,youjustgetusedtothem.HW:p.152#1§2OrthogonalPolynomials&L-SApproximationLab12.OrthogonalPolynomialsApproximationGivenafunctionfandasetof200

m>0distinctpoints.YouaresupposedtowriteafunctionvoidOPA(double(*f)(),doublex[],doublew[],intm,doubletol,FILE*outfile)toapproximatethefunctionfbyanorthogonalpolynomialusingtheexactfunctionvaluesatthegivenmpointsx[].Thearrayw[m]containsthevaluesofaweightfunctionatthegivenpointsx[].Thetotalerrormustbenolargerthantol.AnothervonNeumannqu19§2OrthogonalPolynomials&L-SApproximationInputThereisnoinputfile.Instead,youmusthandinyourfunctionina*.hfile.Theruleofnamingthe*.hfileisthesameasthatofnamingthe*.cor*.cppfiles.Output(representsaspace)Foreachtestcase,youaresupposedtooutputthefollowinginformation:

The1stlinecontainstheinteger6

n>0whichisthedegreeofthepolynomialintheformat:

fprintf(outfile,"%d\n",n);

The2ndlinecontainsthen+1coefficientsoftheapproximationpolynomialwhere.EachofthecoefficientistobeprintedasinCprintf:fprintf(outfile,"%8.4e",coefficient);

The3rdlinecontainsthetotalerrorintheformat:fprintf(outfile,"error=%12.8e\n",err);Note:Ifthetotalerrorisstillnotsmallenoughwhenn=6,simplyoutputtheresultobtainedwhenn=6.Theoutputsoftwotestcasesmustbeseperatedbyablankline.§2OrthogonalPolynomials&L20§2OrthogonalPolynomials&L-SApproximationSampleJudgeProgram#include<stdio.h>#include<math.h>#defineMAX_m200#defineMAX_n6#include"98115001_12.h"

doublef1(doublex){returnsin(x);}

doublef2(doublex){returnexp(x);}

voidmain(){FILE*outfile=fopen("out.txt","w");

intm,i;doublex[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);}§2OrthogonalPolynomials&L21§2OrthogonalPolynomials&L-SApproximationSampleOutput(representsaspace)32.5301e0031.0287e+0007.2279e0021.1287e001error=6.33097847e005

41.0025e+0009.6180e0016.2900e0017.0907e0031.1792e001error=1.61711536e004

§2OrthogonalPolynomials&L22§2函数的最佳逼近/*OptimalApproximation*/

最佳平方逼近：即连续型L-S逼近，在意义下，使得最小。最佳一致逼近/*uniformapproximation*/在意义下，使得最小。也称为minimaxproblem。偏差/*deviation*/若，则称x0为偏差点。Didn’tyousayit’saverydifficultproblem?Takeiteasy.It’snotsodifficultifweconsiderpolynomialsonly.§2函数的最佳逼近/*OptimalApprox23§3OptimalApproximationv1.0最佳一致逼近多项式

/*optimaluniformapproximatingpolynomial*/的构造：求n

阶多项式Pn(x)使得||Pn

y

||最小。直接构造OUAP

的确比较困难，不妨换个角度，先考察它应该具备的性质。有如下结论：

OUAP存在，且必同时有偏差点。证明：存在性证明略。后者用反证法，设只有正偏差点。设而对于所有的x[a,b]都有是n阶多项式是误差更小的多项式§3OptimalApproximationv1.024§3OptimalApproximation（Chebyshev定理）Pn是y的OUAP

Pn关于y在定义域上至少有n+2个交错的偏差点。即存在点集at1<…<tn+2b使得{tk}称为切比雪夫交错组

/*Chebyshevalternatingsequence*/若且y不是n

次多项式，则n次OUAP

唯一。证明：反证，设有2个OUAP’s，分别是Pn

和Qn。则它们的平均函数也是一个OUAP。2)()()(xQxPxRnnn+=对于Rn

有Chebyshev交错组{t1,…,tn+2}使得nkknkknkknnEtytQtytPtytRE-+--=|)()(|21|)()(|21|)()(|nkknkknEtytQtytP=-=-|)()(||)()(|则至少在一个点上必须有)()()()(knkkkntQtytytP-=-0)()(=-kkntytR0=nE§3OptimalApproximation（Ch25§3OptimalApproximation由Chebyshev定理可推出：Pn(x)

y(x)在定义域上至少变号

次，故至少有个根。xy0yyx=()yyxEn=+()yyxEn=-()yPxn=()n+1n+1可见Pn(x)是y(x)的某一个插值多项式

如何确定插值节点{x0,…,xn

}的位置，使得Pn(x)刚好是

y

的OUAP？即，使插值余项v2.0达到极小？§3OptimalApproximation由Ch26§3OptimalApproximationv2.1

在[1,1]上求{x1,…,xn}使得的||wn||最小。=-=niinxxxw1)()(注意到，要使||wn||最小就意味着)()(1xPxxwnnn--=v3.0

在[1,1]上求函数xn的n1阶

OUAP。由Chebyshev定理可推出：Pn1(x)关于xn有n+1个偏差点，即wn(x)在n+1个点上交错取极大、极小值。v3.1

在[1,1]上求切比雪夫交错组{t1,…,tn+1

}。§3OptimalApproximationv2.127切比雪夫多项式/*Chebyshevpolynomials*/§3OptimalApproximation考虑三角函数cos(n)在[0,]上的个极值点。n+1当时，cos(n)交错达到极大值1和极小值1，且存在系数a0,…,an使得

令x=cos()，则x[1,1

]。)cos

arccos()cos()(xn·nxTn==q称为Chebyshev多项式Tn的重要性质：当时，交错取到极大值1和极小值1，即1当时，即{x1,…,xn}为Tn(x)的n个零点。切比雪夫多项式/*Chebyshevpolynom28§3OptimalApproximationTn(x)满足递推关系：T0(x)=1，T1(x)=x，Tn(x)为n

次多项式，首项系数为。且T2n(x)只含x

的次幂，T2n+1(x)只含x

的次幂。2n1偶奇{T0(x),T1(x),…}是[1,1

]上关于权正交的函数族。即在内积的意义下有

OKOK,Ithinkit’senoughforus…What’sourtargetagain?v3.1

在[1,1]上求切比雪夫交错组{t1,…,tn+1

}。v3.0

在[1,1]上求函数xn的n1阶

OUAP。§3OptimalApproximationTn(29Tn(x)的n个零点。§3OptimalApproximation可见：若取，则wn在[1,1

]上有n+1

个极值点{tk}，也即Pn1(x)=xn

wn(x)关于xn在[1,1

]上有n+1个交错偏差点{tk}

。v3.0OKv2.1

在[1,1]上求{x1,…,xn}使得的||wn||最小。=-=niinxxxw1)()(取最小值n={首项系数为1的n

阶多项式/*monicpolynomialsofdegreen*/}{x1,…,xn}即为

如何确定插值节点{x0,…,xn

}的位置，使得Pn(x)刚好是

y

的OUAP？即，使插值余项达到极小？v2.0取{x0,…,xn}为Tn+1(x)的n+1个零点，做y

的插值多项式Pn(x)，则插值余项的上界可达极小。Tn(x)的n个零点。§3OptimalApproxi30§3OptimalApproximation注：上界最小不表示|Rn(x)|最小，故Pn(x)严格意义上只是y(x)的近似最佳逼近多项式；对于一般区间x[a,b]，可作变量替换，则t[1,1

]，这时即以为插值节点(k=0,…,n)，得Pn(x)，余项有最小上界。§3OptimalApproximation注：对31§3OptimalApproximation例：求f(x)=ex在[0,1]上的近似最佳逼近多项式，使其误差不超过0.5104。解：根据误差上界确定n:n=4计算T5(t)的根：以x0,…,x4为节点作L4(x)§3OptimalApproximation例：求f32§3OptimalApproximationChebyshev多项式的其它应用——多项式降次

/*reducethedegreeofpolynomialwithaminimallossofaccuracy*/设f(x)Pn(x)。在降低Pn(x)次数的同时,使因此增加的误差尽可能小,也叫economiza-tionofpowerseries。从Pn中去掉一个含有其最高次项的,结果降次为,则：Pn~Pn1|)(|max|)()(|max|)()(|max]1,1[]1,1[1]1,1[xPxPxfxPxfnnn----+--~因降次而增的误差设Pn的首项系数为an，则取可使精度尽可能少损失。12)()(-=nnnnxTaxP§3OptimalApproximationCh33§3OptimalApproximation例：

f(x)=ex在[1,1]上的4阶Taylor展开为，此时误差请将其降为2阶多项式。解：取（查表知）取（查表知）若简单取，则误差另类解法可查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§3OptimalApproximation例：f34第六章曲线拟合与函数逼近/*ApproximationTheory*/仍然是已知x1…xm

;y1…ym,求一个简单易算的近似函数P(x)

f(x)。但是①

m很大；②

yi本身是测量值，不准确，即yi

f(xi)这时没必要取P(xi)=yi,而要使P(xi)yi总体上尽可能小。常见做法：

使最小/*minimaxproblem*/

太复杂使最小不可导，求解困难使最小/*Least-Squaresmethod*/第六章曲线拟合与函数逼近仍然是已知x1…xm;35§1最小二乘拟合多项式

/*L-Sapproximatingpolynomials*/确定多项式，对于一组数据(xi,yi)(i=1,2,…,n)使得达到极小，这里n

<<

m。naaa10实际上是a0,a1,…,an的多元函数，即[]=-+++=miinininyxaxaaaaa121010...),...,,(j在的极值点应有kiminjijijxyxa==-=10][2-====+njmikiimikjijxyxa0112记====mikiikmikikxycxb11,法方程组(或正规方程组)/*normalequations*/回归系数/*regressioncoefficients*/§1最小二乘拟合多项式/*L-Sapproxim36§1L-SApproximatingPolynomials定理L-S拟合多项式存在唯一

(n<m)。证明：记法方程组为Ba=c.则有其中对任意，必有。若不然，则存在一个使得…即是n

阶多项式的根则B为正定阵，则非奇异，所以法方程组存在唯一解。Waitasecond!Youonlygavemeacriticalpoint,butit’snotnecessarilyaminimumpoint!§1L-SApproximatingPolynomi37§1L-SApproximatingPolynomials定理

Ba=c的解确是的极小点。即：设a

为解，则任意b=(b0

b1…bn)T

对应的多项式必有==njjjxbxF0)(===--=mimiiiiibyxFyxPa1122)(])([])([)(jj证明：==---=-miiimiiiyxPyxFab1212])([])([)()(jj==---+-=miiimiiiiiyxPyxPxPxF1212])([])()()([==--+-=miiiiimiiiyxPxPxFxPxF112])()][()([2)]()([0注：L-Smethod首先要求设定P(x)的形式。若设n=m1，则可取P(x)为过m个点的m1阶插值多项式，这时=0。P(x)不一定是多项式，通常根据经验确定。§1L-SApproximatingPolynomi38例用来拟合。例用39§1L-SApproximatingPolynomials例：xy(xi,yi),i=1,2,…,m方案一：设baxxxPy+=)(求a和b使得最小。=-+=miiiiybaxxba12)(),(jButhey,thesystemofequationsforaandbisnonlinear!Takeiteasy!Wejusthavetolinearizeit…线性化

/*linearization*/：令，则bXaY+就是个线性问题将化为后易解a和b。),(iiYX),(iiyx§1L-SApproximatingPolynomi40例用来拟合。例用41§1L-SApproximatingPolynomials方案二：设xbeaxPy/)(-=(a>0,b>0)线性化：由可做变换xbay-lnlnbBaAxXyY-====,ln,1,lnBXAY+就是个线性问题将化为后易解A和B),(iiYX),(iiyxHW:p.233#7，#9,#10,#11§1L-SApproximatingPolynomi42例用来拟合。例用43§2正交多项式与最小二乘拟合

/*OrthogonalPolynomials&Least-SquaresApproximation*/已知x1…xm

;y1…ym,求一个简单易算的近似函数P(x)

f(x)使得最小。已知[a,b]上定义的f(x)，求一个简单易算的近似函数P(x)使得最小。定义

线性无关/*linearlyindependent*/函数族{0(x),1(x),…,n(x),…}满足条件：其中任意函数的线性组合

a00(x)+a11(x)+…+ann(x)=0对任意x[a,b]成立当且仅当a0=a1=…=an=0。§2正交多项式与最小二乘拟合已知x1…xm;y44§2OrthogonalPolynomials&L-SApproximation定义考虑一般的线性无关函数族={0(x),1(x),…,n(x),…}，其有限项的线性组合称为广义多项式

/*generalizedpolynomial*/.常见多项式：

{j(x)=xj}对应代数多项式/*algebraicpolynomial*/

{j(x)=cosjx}、{j(x)=sinjx}{j(x),j(x)

}对应三角多项式/*trigonometricpolynomial*/

{j(x)=ekjx,ki

kj

}对应指数多项式/*exponentialpolynomial*/§2OrthogonalPolynomials&L45§2OrthogonalPolynomials&L-SApproximation定义权函数：①

离散型/*discretetype*/根据一系列离散点拟合时，在每一误差前乘一正数wi

，即误差函数

，这个wi

就称作权/*weight*/，反映该点的重要程度。=-=niiiiyxPw12])([②

连续型

/*continuoustype*/在[a,b]上用广义多项式P(x)拟合连续函数f(x)时，定义权函数(x)C[a,b]，即误差函数=。权函数必须(x)满足：非负、可积，且在[a,b]的任何子区间上(x)0。§2OrthogonalPolynomials&L46§2OrthogonalPolynomials&L-SApproximation定义广义L-S拟合：①

离散型/*discretetype*/在点集{x1…xm}

上测得{y1…ym}，在一组权系数{w1…wm}下求广义多项式P(x)使得误差函数最小。

=-=niiiiyxPw12])([②

连续型

/*continuoustype*/已知y(x)

C[a,b]以及权函数(x)，求广义多项式P(x)使得误差函数=最小。dxxyxPxba2)]()([)(-r内积与范数离散型连续型则易证(f,g)是内积，而是范数。(f,g)=0表示f与g

带权正交。广义L-S问题可叙述为：求广义多项式P(x)使得最小。§2OrthogonalPolynomials&L47§2OrthogonalPolynomials&L-SApproximationnkyaknjjjk,...,0,),(),(0===jjj设则完全类似地有：)(...)()()(1100xaxaxaxPnnjjj+++=法方程组

/*normalequations*/定理

Ba=c存在唯一解

0(x),1(x),…,n(x)线性无关。即：),(),(),(00yyaabnnjiijjjjj===c证明：若存在一组系数{i

}使得0...1100=+++nnjajaja则等式两边分别与0,1,…,n作内积，得到：即：B=0……§2OrthogonalPolynomials&L48§2OrthogonalPolynomials&L-SApproximation例：用来拟合，w1解：0(x)=1，1(x)=x，2(x)=x2Itissoooosimple!Whatcanpossiblygowrong?7623)(463||||484,||||1==-=BcondBB§2OrthogonalPolynomials&L49§2OrthogonalPolynomials&L-SApproximation例：连续型拟合中，取则Hilbert阵！改进：若能取函数族={0(x),1(x),…,n(x),…}，使得任意一对i(x)和j(x)两两（带权）正交，则B就化为对角阵！这时直接可算出ak=Well,nofreelunchanyway…

正交多项式的构造：将正交函数族中的k取为k阶多项式，为简单起见，可取k的首项系数为1。有递推关系式：其中证明略§2OrthogonalPolynomials&L50§2OrthogonalPolynomials&L-SApproximation例：用来拟合，w1解：通过正交多项式0(x),1(x),2(x)求解设)()()(221100xaxaxayjjj++=1)(0=xj229),(),(0000==jjjya25),(),(00001==jjjjax25)()()(011-=-=xxxxjaj537),(),(1111==jjjya25),(),(11112==jjjjax45),(),(00111==jjjjb55)(45)()25()(2012+-=--=xxxxxxjjj21),(),(2222==jjjya与前例结果一致。注：手算时也可用待定系数法确定函数族。§2OrthogonalPolynomials&L51§2OrthogonalPolynomials&L-SApproximation

Algorithm:OrthogonalPolynomialsApproximation

Toapproximateagivenfunctionbyapolynomialwitherrorboundedbyagiventolerance.Input:numberofdatam;x[m];y[m];weightw[m];toleranceTOL;maximumdegreeofpolynomialMax_n.Output:coefficientsoftheapproximatingpolynomial.Step1Set0(x)

1;a0=(0,y)/(0,0);P(x)=a00(x);err=(y,y)a0(0,y);Step2Set1=

(x0,0)/(0,0);1(x)

=(x1)0(x);

a1=(1,y)/(1,1);P(x)+=a11(x);err

=a1(1,y);Step3Setk=1;Step4While((k<Max_n)&&(|err|TOL))dosteps5-7

Step5k++;

Step6k=

(x1,1)/(1,1);k1=(1,1)/(0,0);

2(x)

=(xk)1(x)k10(x);ak

=(2,y)/(2,2);

P(x)+=ak

2(x);err

=ak

(2,y);

Step7Set0(x)=1(x);1(x)=2(x);Step8Output();STOP.注：§2OrthogonalPolynomials&L52AnothervonNeumannquote:Youngman,inmathematicsyoudon'tunderstandthings,youjustgetusedtothem.HW:p.152#1§2OrthogonalPolynomials&L-SApproximationLab12.OrthogonalPolynomialsApproximationGivenafunctionfandasetof200

m>0distinctpoints.YouaresupposedtowriteafunctionvoidOPA(double(*f)(),doublex[],doublew[],intm,doubletol,FILE*outfile)toapproximatethefunctionfbyanorthogonalpolynomialusingtheexactfunctionvaluesatthegivenmpointsx[].Thearrayw[m]containstheva