第七讲MATLAB在计算方法中的应用

MATLAB在计算方法中的应用MATLAB入门到精通插值与拟合插值与拟合来源于实际，又广泛应用于实际。随着计算机的不断发展及计算水平的不断提高，他们在国民生产和科学研究等方面扮演着越来越重要的角色。插值法主要有Lagrange插值、分段线性插值、Hermite插值及三次样条插值等拟合法主要有最小二乘法拟合和快速Fourier变换等Lagrange插值对给定的n个插值节点及相应的函数值，利用n次Lagrange插值多项式对插值区间内任意x的函数值y可通过下式求得：nkjkjnkjjkxxxxyxy11)()(MATLAB实现functiony=lagrange(x0,y0,x)%lagrangeinsertn=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;end例：f(x)=lnxx0.40.50.60.70.8ln(x)-0.916291-0.693147-0.510826-0.357765-0.223144x=[0.4:0.1:0.8];y=[-0.916291-0.693147-0.510826-0.356675-0.223144];lagrange(x,y,0.54)ans=-0.61614log(0.54)ans=-0.61619Runge现象19世纪Runge给出了一个等距节点插值多项式不收敛的例子：在区间[-5，5]上各阶导数存在，但在此区间上取n个节点构造的Lagrange插值多项式在全区间上不收敛。211)(xxfRunge现象在区间[-5，5]上，取n=10，用lagrange插值法进行插值计算:x=[-5:1:5];y=1./(1+x.^2);x0=[-5:0.1:5];y0=lagrange(x,y,x0);y1=1./(1+x0.^2);plot(x0,y0,'--r')holdonplot(x0,y1,'-b')-5-4-3-2-1012345-0.500.511.52分段线性插值分段线性插值就是通过插值点用折线段连接起来逼近原曲线。MATLAB实现interp1一维插值yi=interp1(x,y,xi)%对一组节点(x,y)进行插值，计算插值点xi的函数值。yi=interp1(y,xi)%默认x=1:nyi=interp1(x,y,xi,’method’)%method为指定插值算法nearest线性最近项插值linear线性插值spline三次样条插值cubic三次插值同类的函数还有inter1q,interpft,spline,interp2,interp3,interpn等例：正弦曲线插值x=0:0.1:10;y=sin(x);xi=0:.25:10;yi=interp1(x,y,xi);plot(x,y,'o',xi,yi)0246810-1-0.500.51例：Rouge现象的解决x=[-5:1:5];y=1./(1+x.^2);x0=[-5:0.1:5];y1=1./(1+x0.^2);y2=interp1(x,y,x0);plot(x0,y0,'--r')holdonplot(x0,y1,'-b')plot(x0,y2,'*m')-505-0.500.511.52Hermite插值不少问题不但要求在节点上函数值相等，而且要求导数值也相等，甚至要求高阶导数也相等，可用Hermite插值多项式。已知n个插值节点及其对应的函数值和一阶导数值，则计算插值区域内任意x的函数值y的Hermite插值公式为：nijijijijnijjiniiiiiiixxaxxxxhyyyaxxhxy12111;)(])2)([()(其中MATLAB实现functiony=hermite(x0,y0,y1,x)%hermiteinsertn=length(x0);m=length(x);fork=1:myy=0.0;fori=1:nh=1.0;a=0.0;forj=1:nifj~=ih=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2;a=1/(x0(i)-x0(j))+a;endendyy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i));endy(k)=yy;end例：根据如下给定数据，给出0.34处的值X0.300.320.35Sinx0.295520.314570.34290Cosx0.955340.949240.93937x0=[0.30.320.35];y0=[0.295520.314570.34290];y1=[0.955340.949240.93937];x=[0.3:0.005:0.35];y=hermite(x0,y0,y1,0.34)y=0.33349sin(0.34)ans=0.33349y=hermite(x0,y0,y1,x);plot(x,y)holdonplot(x,sin(x),'--r')0.30.3050.310.3150.320.3250.330.3350.340.3450.350.2950.30.3050.310.3150.320.3250.330.3350.340.345三次样条插值样条函数可以给出平滑的插值曲线和曲面。方法介绍：设在[a,b]上给定n+1个点a=x0x1x2…xn=b和相应的函数值yi=f(xi),i=0,1,…,n求三次样条函数s(x)满足插值条件:s(xi)=yi以及下列三类边界条件之一：Is′(x0)=y0′,s′(xn)=yn′Ⅱs″(x0)=y0″,s″(xn)=yn″IIIs′(x0)=s′(xn),s″(x0)=s″(xn)当f(x0)=f(xn)时此s(x)称为三次样条插值函数.MATLAB实现SPLINECubicsplinedatainterpolation.YY=SPLINE(X,Y,XX)usescubicsplineinterpolationtofindavectorYYcorrespondingtoXX.XandYarethegivendatavectorsandXXisthenewabscissavector.TheordinatesYmaybevector-valued,inwhichcaseY(:,j)isthej-thordinate.PP=SPLINE(X,Y)returnstheppformofthecubicsplineinterpolantinstead,forlaterusewithPPVAL,etc.Ordinarily,thenot-a-knotendconditionsareused.However,ifYcontainstwomoreordinatesthanXhasentries,thenthefirstandlastordinateinYareusedastheendslopesforthecubicspline.例：正弦函数x=0:10;y=sin(x);xx=0:.25:10;yy=spline(x,y,xx);plot(x,y,'o',xx,yy)0246810-1-0.500.51例：如下表所示数据x28.7282930f(x)4.14.34.13.0x=[2828.72930];y=[4.14.34.13.0];x0=[28:0.15:30];y1=spline(x,y,x0);plot(x,y,'--r')holdonplot(x0,y1)2828.228.428.628.82929.229.429.629.83033.544.5最小二乘法拟合在科学实验的统计方法中，往往要从一组实验数据(xi,yi)中寻找出自变量x和因变量y之间的函数关系y=f(x)。由于观测数据往往不够准确，因此并不要求y=f(x)经过所有的点(xi,yi)，而只要求在给定点xi上误差ei=f(xi)-yi按照某种标准达到最小，通常采用欧氏范数||e||2作为误差度量的标准，这就是最小二乘法。MATLAB实现利用polyfit函数进行多项式拟合POLYFITFitpolynomialtodata.POLYFIT(X,Y,N)findsthecoefficientsofapolynomialP(X)ofdegreeNthatfitsthedata,P(X(I))~=Y(I),inaleast-squaressense.[P,S]=POLYFIT(X,Y,N)returnsthepolynomialcoefficientsPandastructureSforusewithPOLYVALtoobtainerrorestimatesonpredictions.Iftheerrorsinthedata,Y,areindependentnormalwithconstantvariance,POLYVALwillproduceerrorboundswhichcontainatleast50%ofthepredictions.ThestructureScontainstheCholeskyfactoroftheVandermondematrix(R),thedegreesoffreedom(df),andthenormoftheresiduals(normr)asfields.利用矩阵除法解决复杂型函数的拟合例：拟合以下数据x0.51.01.52.02.53.0y1.752.453.814.808.008.60?x=[0.51.01.52.02.53.0];?y=[1.752.453.814.808.008.60];?a=polyfit(x,y,2)a=0.49001.25010.8560?x1=[0.5:0.05:3.0];?y1=a(3)+a(2)*x1+a(1)*x1.^2;?plot(x,y,'*')?holdon?plot(x1,y1,'-r')0.511.522.530246810例：根据经验公式y=a+bx2，拟合如下数据：xi1925313844yi19.032.349.073.398.8?x=[1925313844];?y=[19.032.349.073.398.8];?x1=x.^2;?x1=[ones(5,1),x1'];?ab=x1\y'ab=0.59370.0506?x0=[19:0.2:44];?y0=ab(1)+ab(2)*x0.^2;?plot(x,y,'o',x0,y0,'-r')15202530354045020406080100快速Fourier变换在函数逼近中，当函数为周期函数时，用三角多项式逼近比用代数多项式更合适，因此引入Fourier逼近和快速Fourier变换。MATLAB实现fft快速Fourier变换fft(x)fft(x,n)fft(x,[],DIM)或fft(x,n,DIM)fft2二维快速Fourier变换fft2(x)fft2(x,MRows,NCols)fftnn维快速Fourier变换fftn(x)fftn(x,size)ifft快速Fourier逆变换ifft2二维快速Fourier逆变换ifftnn维快速Fourier逆变换积分与微分Newton-Cotes系列数值求积将积分区间[a,b]划分为n等分，步长h=(b-a)/n，选取等距节点xk=a+kh，则求积分公式：Ck(n)由Cotes系数表给出。n=0时为矩形公式；n=1时为梯形公式;n=2时为Simpson公式