計算方法曲線擬合

計算方法曲線擬合第1頁，共16頁，2023年，2月20日，星期四1.4OrthogonalPolynomialsL-SApproximation則稱0(x),1(x),…,n(x)是關(guān)于點集{xi}的正交函數(shù)族/*orthogonalfunctions

*/。如果k(x)都是多項式時，就稱之為關(guān)于點集{xi}的正交多項式/*orthogonal

polynomials

*/。對于點集{xi}(i=1,2,…,m)，若一組函數(shù)0(x),1(x),…,n(x)(n<m)滿足條件定義第2頁，共16頁，2023年，2月20日，星期四構(gòu)造正交多項式1.4OrthogonalPolynomialsL-SApproximation定理對于給定的點集{xi}(i=1,2,…,m)，遞推公式構(gòu)造的函數(shù)族0(x),1(x),…,n(x)是關(guān)于{xi}的一組正交多項式，且k(x)(k=0,1,…,n)是k次多項式，其最高項xn的系數(shù)為1。注：構(gòu)造的次序是：

證明略…第3頁，共16頁，2023年，2月20日，星期四例：用來擬合。解：通過正交多項式0(x),1(x),2(x)求解設(shè))()()(221100xaxaxayjjj++=1)(0=xj229),(),(0000==jjjya25),(),(00001==jjjjax25)()()(011-=-=xxxxjaj537),(),(1111==jjjya25),(),(11112==jjjjax45),(),(00111==jjjjb55)(45)()25()(2012+-=--=xxxxxxjjj21),(),(2222==jjjya1.4OrthogonalPolynomialsL-SApproximation第4頁，共16頁，2023年，2月20日，星期四程序設(shè)計uk+1,vk

,pi

,

qi

,

s1,

s0

1.4OrthogonalPolynomialsL-SApproximation第5頁，共16頁，2023年，2月20日，星期四程序設(shè)計(續(xù))1.4OrthogonalPolynomialsL-SApproximation第6頁，共16頁，2023年，2月20日，星期四§2最佳一致逼近

/*Optimal

UniformApproximation*/

對于一組數(shù)據(jù)(xi,yi),i=1,2,…,m,要求在某個函數(shù)類

(n<<m)中構(gòu)造一個函數(shù)，使得最大偏差最小，即ak在I的極值點引入變量D，則最大偏差最小化問題可歸結(jié)為以下問題：

besubjectedto以…為條件線性規(guī)劃問題

linearprogramming第7頁，共16頁，2023年，2月20日，星期四§2Optimal

UniformApproximation

線性規(guī)劃規(guī)劃變量

programming

variables目標函數(shù)

objectivfunction(線性)約束方程

constraintequations單純形法、有效集法

第8頁，共16頁，2023年，2月20日，星期四§2Optimal

UniformApproximation

解仍然采用線性擬合φ(x)=a0+a1x。相應(yīng)的線性規(guī)劃問題為：例某種鋁合金的含鋁量為x(%)，熔解溫度為y(oC)，由實驗測得的數(shù)據(jù)如下表。試用最佳一致逼近建立x與y之間的經(jīng)驗公式。

xi36.946.763.777.884.087.5yi181197235270283292規(guī)劃變量為a0,a1,D用單純法解得：a0=96.4886,a1=2.2132,D=2.8450,故有就

φ(x)=96.4886+2.2132x最小二乘法xi36.946.763.777.884.087.5φ(xi)178.16199.85237.47268.68282.40290.14δi-2.852.852.47-1.32-0.60-1.86擬合方法均方誤差最大偏差最小二乘法5.163.22最佳一致逼近6.282.85第9頁，共16頁，2023年，2月20日，星期四§2Optimal

UniformApproximation

例某圓形工件圓周上的測量數(shù)據(jù)為(xi,yi),i=1,2,…,m，求該工件的圓心坐標和半徑。

為了確定待定量x0,y0

和r，讓最大偏差最小化：

設(shè)(x0’,y0’)是理想圓心坐標的近似值，把偏差表達式在x0’,y0’附近展成泰勒級數(shù)，并取線性項得到

解由于加工誤差存在，工件輪廓非標準圓。設(shè)理想圓的圓心坐標為(x0,y0),半徑為r，則測點到理想圓的偏差為It’seasythatthecirclecanbedeterminedthroughanythreepoints這是個非線性優(yōu)化問題，求解較困難，需線性化。

第10頁，共16頁，2023年，2月20日，星期四§2Optimal

UniformApproximation

用單純形法或有效集法解出規(guī)劃變量Δx、Δy、r和D，便有x0=x0’+Δx,y0=y0’+Δy和r

，問題得解。設(shè)D=于是有如下線性規(guī)劃問題：第11頁，共16頁，2023年，2月20日，星期四§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.注：第12頁，共16頁，2023年，2月20日，星期四

AnothervonNeumannquote: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.第13頁，共16頁，2023年，2月20日，星期四§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.第14頁，共16頁，2023年，2月20日，星期四§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=