1、MATLAB在计算方法中的应用,MATLAB入门到精通,插值与拟合,插值与拟合来源于实际，又广泛应用于实际。随着计算机的不断发展及计算水平的不断提高，他们在国民生产和科学研究等方面扮演着越来越重要的角色。插值法主要有Lagrange插值、分段线性插值、Hermite插值及三次样条插值等拟合法主要有最小二乘法拟合和快速Fourier变换等,Lagrange插值,对给定的n个插值节点及相应的函数值，利用n次Lagrange插值多项式对插值区间内任意x的函数值y可通过下式求得：,MATLAB实现,function y=lagrange(x0,y0,x)%lagrange insertn=length
2、(x0);m=length(x);for i=1:m z=x(i);s=0.0;for k=1:n p=1.0;for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j);end end s=p*y0(k)+s;end y(i)=s;end,例：f(x)=lnx,x=0.4:0.1:0.8;y=-0.916291-0.693147-0.510826-0.356675-0.223144;lagrange(x,y,0.54)ans=-0.61614 log(0.54)ans=-0.61619,Runge现象,19世纪Runge给出了一个等距节点插值多项式不收敛的例子：在区
3、间-5，5上各阶导数存在，但在此区间上取n个节点构造的Lagrange插值多项式在全区间上不收敛。,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)hold on plot(x0,y1,-b),分段线性插值,分段线性插值就是通过插值点用折线段连接起来逼近原曲线。MATLAB实现 interp1 一维插值yi=interp1(x,y,xi)%对一组节点(x,y)进行插值，计算插值点xi的函数值。y
4、i=interp1(y,xi)%默认x=1:nyi=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),例：Rouge现象的解决,x=-5:1:5;y=1./(1+x.2);x0=-5:0.1:5;y1=1./(1
5、+x0.2);y2=interp1(x,y,x0);plot(x0,y0,-r)hold on plot(x0,y1,-b)plot(x0,y2,*m),Hermite插值,不少问题不但要求在节点上函数值相等，而且要求导数值也相等，甚至要求高阶导数也相等，可用Hermite插值多项式。已知n个插值节点及其对应的函数值和一阶导数值，则计算插值区域内任意x的函数值y的Hermite插值公式为：,MATLAB实现,function y=hermite(x0,y0,y1,x)%hermite insertn=length(x0);m=length(x);for k=1:m yy=0.0;for i=1
6、:n h=1.0;a=0.0;for j=1:n if j=i h=h*(x(k)-x0(j)/(x0(i)-x0(j)2;a=1/(x0(i)-x0(j)+a;end end yy=yy+h*(x0(i)-x(k)*(2*a*y0(i)-y1(i)+y0(i);end y(k)=yy;end,例：根据如下给定数据，给出0.34处的值,x0=0.3 0.32 0.35;y0=0.29552 0.31457 0.34290;y1=0.95534 0.94924 0.93937;x=0.3:0.005:0.35;y=hermite(x0,y0,y1,0.34)y=0.33349 sin(0.34)
7、ans=0.33349 y=hermite(x0,y0,y1,x);plot(x,y)hold on plot(x,sin(x),-r),三次样条插值,样条函数可以给出平滑的插值曲线和曲面。方法介绍：,设在a,b上给定n+1个点 a=x0 x1x2xn=b和相应的函数值yi=f(xi),i=0,1,n求三次样条函数s(x)满足插值条件:s(xi)=yi以及下列三类边界条件之一：I s(x0)=y0,s(xn)=yn s(x0)=y0,s(xn)=yn III s(x0)=s(xn),s(x0)=s(xn)当f(x0)=f(xn)时 此s(x)称为三次样条插值函数.,MATLAB实现,SPLIN
8、E Cubic spline data interpolation.YY=SPLINE(X,Y,XX)uses cubic spline interpolation to find a vector YY corresponding to XX.X and Y are the given data vectors and XX is the new abscissa vector.The ordinates Y may be vector-valued,in which case Y(:,j)is the j-th ordinate.PP=SPLINE(X,Y)returns the ppfo
9、rm of the cubic spline interpolant instead,for later use with PPVAL,etc.Ordinarily,the not-a-knot end conditions are used.However,if Y contains two more ordinates than X has entries,then the first and last ordinate in Y are used as the endslopes for the cubic spline.,例：正弦函数,x=0:10;y=sin(x);xx=0:.25:
10、10;yy=spline(x,y,xx);plot(x,y,o,xx,yy),例：如下表所示数据,x=28 28.7 29 30;y=4.1 4.3 4.1 3.0;x0=28:0.15:30;y1=spline(x,y,x0);plot(x,y,-r)hold onplot(x0,y1),最小二乘法拟合,在科学实验的统计方法中，往往要从一组实验数据(xi,yi)中寻找出自变量x和因变量y之间的函数关系y=f(x)。由于观测数据往往不够准确，因此并不要求y=f(x)经过所有的点(xi,yi)，而只要求在给定点xi上误差ei=f(xi)-yi按照某种标准达到最小，通常采用欧氏范数|e|2作为误差
11、度量的标准，这就是最小二乘法。,MATLAB实现,利用polyfit函数进行多项式拟合POLYFIT Fit polynomial to data.POLYFIT(X,Y,N)finds the coefficients of a polynomial P(X)of degree N that fits the data,P(X(I)=Y(I),in a least-squares sense.P,S=POLYFIT(X,Y,N)returns the polynomial coefficients P and a structure S for use with POLYVAL to obt
12、ain error estimates on predictions.If the errors in the data,Y,are independent normal with constant variance,POLYVAL will produce error bounds which contain at least 50%of the predictions.The structure S contains the Cholesky factor of the Vandermonde matrix(R),the degrees of freedom(df),and the nor
13、m of the residuals(normr)as fields.利用矩阵除法解决复杂型函数的拟合,例：拟合以下数据,?x=0.5 1.0 1.5 2.0 2.5 3.0;?y=1.75 2.45 3.81 4.80 8.00 8.60;?a=polyfit(x,y,2)a=0.4900 1.2501 0.8560?x1=0.5:0.05:3.0;?y1=a(3)+a(2)*x1+a(1)*x1.2;?plot(x,y,*)?hold on?plot(x1,y1,-r),例：根据经验公式y=a+bx2，拟合如下数据：,?x=19 25 31 38 44;?y=19.0 32.3 49.0
14、73.3 98.8;?x1=x.2;?x1=ones(5,1),x1;?ab=x1yab=0.5937 0.0506?x0=19:0.2:44;?y0=ab(1)+ab(2)*x0.2;?plot(x,y,o,x0,y0,-r),快速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)ff
15、tn n维快速Fourier变换 fftn(x)fftn(x,size)ifft 快速Fourier逆变换ifft2 二维快速Fourier逆变换ifftn n维快速Fourier逆变换,积分与微分,Newton-Cotes系列数值求积将积分区间a,b划分为n等分，步长h=(b-a)/n，选取等距节点xk=a+kh，则求积分公式：Ck(n)由Cotes系数表给出。n=0 时为矩形公式；n=1时为梯形公式;n=2时为Simpson公式；n=3时为Cotes公式；,矩形求积公式Cumsum 元素累积和函数Cumsum(x)对于向量x，函数返回一个向量，向量的第n个元素为向量x的前n个元素的和；对于
16、矩阵X，函数返回一个矩阵，矩阵的列对应X的列的累积和返回值Cumsum(x,DIM)参数DIM指明求和从第DIM维开始。,例：x=0 1 2;3 4 5;cumsum(x)ans=0 1 2 3 5 7,梯形求积分公式Z=trapz(Y)用梯形方法计算Y的积分近似值。对向量Y，函数返回Y的积分；对矩阵Y，函数返回一个向量，向量各元素为矩阵Y对应列向量的积分值。Z=trapz(X,Y)Z=trapz(X,Y,DIM)Z=trapz(Y,DIM),例：求积分,?d=pi/1000;?t=0:d:3*pi;?nt=length(t);?y=fun(t);?sc=cumsum(y)*d;?scf=sc
17、(nt)scf=0.901618619310013?z=trapz(y)*dz=0.900840276606884,建立被积函数：Fun.mfunction y=fun(t)y=exp(-0.5*t).*sin(t+pi/6);,自适应法-simpson法求积,QuadQ=quad(F,a,b)求函数F(x)从a到b的相对误差为1e-3的积分近似值；Q=quad(F,a,b,tol)向量tol给出相对误差和绝对误差；Q=quad(F,a,b,tol,TRACE)TRACE不为零时给出图形；,例：求积分,建立被积函数：fun1.mfunction f=fun1(x)f=x./(4+x.2);,在
18、MATLAB命令窗口中：?quad(fun1,0,1)ans=0.111572382538912,自适应法-cotes法求积,Quad8Q=quad8(F,a,b)求函数F(x)从a到b的相对误差为1e-3的积分近似值；Q=quad8(F,a,b,tol)向量tol给出相对误差和绝对误差；Q=quad8(F,a,b,tol,TRACE)TRACE不为零时给出图形；,例：求积分,建立被积函数：fun2.mfunction f=fun2(x)f=exp(-x/2);,在MATLAB命令窗口：?quad8(fun2,1,3,1e-10,1)ans=0.766800999128407,Gauss求积公
19、式,为了得到较高的代数精度，可以使用Gauss公式：式中Ak为Gauss-Legender系数对任意区间a,b的积分，可变换为如下公式：,Romberg求积公式,梯形公式积分的算法简单、编程容易，但精度较差、收敛速度较慢；Romberg求积公式可自动改变积分步长，使相邻的两次值的绝对误差或相对误差小于预先设定的允许误差。,符号积分,SymbolicToolBoxint(S)%计算符号表达式的不定积分；int(S,v)%计算符号表达式对自变量v的不定积分；int(S,a,b)%计算从a到b的定积分；int(S,v,a,b)%计算对变量v从a到b的定积分；,符号积分,?syms x t;?a=co
20、s(x*t),sin(x*t);-sin(x*t),cos(x*t);?int(a,t)ans=1/x*sin(x*t),-cos(x*t)/x cos(x*t)/x,1/x*sin(x*t),?syms x1?int(x1*log(1+x1),0,1)ans=1/4,微分和差分,数值微分和差分Diff(x)%对向量x，其值为 x(2)-x(1)x(3)-x(2)x(n)-x(n-1)Diff(X)%对矩阵X,其值为矩阵列的差分 x(2:n,:)-x(1:n-1,:)符号微分和差分Diff(S)%对符号表达式S求微分Diff(S,v)%对自变量v求微分Diff(S,n)%n为正整数，对S求n次
21、微分Diff(S,v,n),?x=sym(x)；?diff(sin(x2)ans=2*cos(x2)*x,?syms x?diff(sin(x2),2)ans=-4*sin(x2)*x2+2*cos(x2),多元函数导数Jacobian(f,v)%计算对向量v的Jacobi矩阵,例：,?f2=3*x-cos(x*y)-0.5;x2-81*(y+0.1)2+sin(z)+1.06;exp(-x*y)+20*z+(10*pi/3-1)；?jacobian(f2,x,y,z)ans=3+sin(x*y)*y,sin(x*y)*x,0 2*x,-162*y-16.2,cos(z)-y*exp(-x*y),-x*exp(-x*y),20,