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

MATLAB在计算方法中的应用MATLAB入门到精通MATLAB在计算方法中的应用MATLAB入门到精通1插值与拟合l插值与拟合来源于实际，又广泛应用于实际。随着计算机的不断发展及计算水平的不断提高，他们在国民生产和科学研究等方面扮演着越来越重要的角色。l插值法主要有Lagrange插值、分段线性插值、Hermite插值及三次样条插值等l拟合法主要有最小二乘法拟合和快速Fourier变换等插值与拟合插值与拟合来源于实际，又广泛应用于实际。随着计算机2Lagrange插值l对给定的n个插值节点及相应的函数值，利用n次Lagrange插值多项式对插值区间内任意x的函数值y可通过下式求得：Lagrange插值对给定的n个插值节点及相应的函数值，利用3MATLAB实现function y=lagrange(x0,y0,x)%lagrange insertn=length(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;endMATLAB实现function y=lagrange(x04例：f(x)=lnxx0.40.50.60.70.8ln(x)-0.916291-0.693147-0.510826-0.357765-0.223144 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例：f(x)=lnxx0.40.50.60.70.8ln(x5Runge现象l19世纪Runge给出了一个等距节点插值多项式不收敛的例子：在区间-5，5上各阶导数存在，但在此区间上取n个节点构造的Lagrange插值多项式在全区间上不收敛。Runge现象19世纪Runge给出了一个等距节点插值多项式6Runge现象l在区间-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)Runge现象在区间-5，5上，取n=10，用lagra7分段线性插值l分段线性插值就是通过插值点用折线段连接起来逼近原曲线。lMATLAB实现 interp1 一维插值yi=interp1(x,y,xi)%对一组节点(x,y)进行插值，计算插值点xi的函数值。yi=interp1(y,xi)%默认x=1:nyi=interp1(x,y,xi,method)%method为指定插值算法lnearest 线性最近项插值llinear 线性插值lspline 三次样条插值lcubic 三次插值同类的函数还有inter1q,interpft,spline,interp2,interp3,interpn等分段线性插值分段线性插值就是通过插值点用折线段连接起来逼近原8例：正弦曲线插值 x=0:0.1:10;y=sin(x);xi=0:.25:10;yi=interp1(x,y,xi);plot(x,y,o,xi,yi)例：正弦曲线插值 x=0:0.1:10;9例：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)hold on plot(x0,y1,-b)plot(x0,y2,*m)例：Rouge现象的解决 x=-5:1:5;10Hermite插值l不少问题不但要求在节点上函数值相等，而且要求导数值也相等，甚至要求高阶导数也相等，可用Hermite插值多项式。l已知n个插值节点及其对应的函数值和一阶导数值，则计算插值区域内任意x的函数值y的Hermite插值公式为：Hermite插值不少问题不但要求在节点上函数值相等，而且要11MATLAB实现function y=hermite(x0,y0,y1,x)%hermite insertn=length(x0);m=length(x);for k=1:m yy=0.0;for i=1: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;endMATLAB实现function y=hermite(x0,12例：根据如下给定数据，给出0.34处的值X0.300.320.35Sinx0.295520.314570.34290Cosx0.955340.949240.93937 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)ans=0.33349 y=hermite(x0,y0,y1,x);plot(x,y)hold on plot(x,sin(x),-r)例：根据如下给定数据，给出0.34处的值X0.300.32013三次样条插值l样条函数可以给出平滑的插值曲线和曲面。l方法介绍：设在a,b上给定n+1个点 a=x0 x1x2x=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)例：如下表所示数据x28.7282930f(x)4.14.317最小二乘法拟合l在科学实验的统计方法中，往往要从一组实验数据(xi,yi)中寻找出自变量x和因变量y之间的函数关系y=f(x)。l由于观测数据往往不够准确，因此并不要求y=f(x)经过所有的点(xi,yi)，而只要求在给定点xi上误差ei=f(xi)-yi按照某种标准达到最小，通常采用欧氏范数|e|2作为误差度量的标准，这就是最小二乘法。最小二乘法拟合在科学实验的统计方法中，往往要从一组实验数据(18MATLAB实现l利用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 obtain 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 norm of the residuals(normr)as fields.l利用矩阵除法解决复杂型函数的拟合MATLAB实现利用polyfit函数进行多项式拟合19例：拟合以下数据x0.51.01.52.02.53.0y1.752.453.814.808.008.60?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)例：拟合以下数据x0.51.01.52.02.53.0y1.20例：根据经验公式y=a+bx2，拟合如下数据：xi1925313844yi19.032.349.073.398.8?x=19 25 31 38 44;?y=19.0 32.3 49.0 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)例：根据经验公式y=a+bx2，拟合如下数据：xi1925321快速Fourier变换l在函数逼近中，当函数为周期函数时，用三角多项式逼近比用代数多项式更合适，因此引入Fourier逼近和快速Fourier变换。lMATLAB实现 fft 快速Fourier变换l fft(x)l fft(x,n)l fft(x,DIM)或 fft(x,n,DIM)快速Fourier变换在函数逼近中，当函数为周期函数时，用三22 fft2 二维快速Fourier变换l fft2(x)l fft2(x,MRows,NCols)fftn n维快速Fourier变换l fftn(x)l fftn(x,size)ifft 快速Fourier逆变换ifft2 二维快速Fourier逆变换ifftn n维快速Fourier逆变换 fft2 二维快速Fourier变换23积分与微分lNewton-Cotes系列数值求积将积分区间a,b划分为n等分，步长h=(b-a)/n，选取等距节点xk=a+kh，则求积分公式：Ck(n)由Cotes系数表给出。ln=0 时为矩形公式；ln=1时为梯形公式;ln=2时为Simpson公式；ln=3时为Cotes公式；l积分与微分Newton-Cotes系列数值求积24l矩形求积公式lCumsum 元素累积和函数Cumsum(x)对于向量x，函数返回一个向量，向量的第n个元素为向量x的前n个元素的和；对于矩阵X，函数返回一个矩阵，矩阵的列对应X的列的累积和返回值Cumsum(x,DIM)参数DIM指明求和从第DIM维开始。例：x=0 1 2;3 4 5;cumsum(x)ans=0 1 2 3 5 7矩形求积公式例：x=0 1 2;3 4 5;25l梯形求积分公式lZ=trapz(Y)用梯形方法计算Y的积分近似值。对向量Y，函数返回Y的积分；对矩阵Y，函数返回一个向量，向量各元素为矩阵Y对应列向量的积分值。lZ=trapz(X,Y)lZ=trapz(X,Y,DIM)lZ=trapz(Y,DIM)梯形求积分公式26例：求积分?d=pi/1000;?t=0:d:3*pi;?nt=length(t);?y=fun(t);?sc=cumsum(y)*d;?scf=sc(nt)scf=0.901618619310013?z=trapz(y)*dz=0.900840276606884建立被积函数：Fun.mfunction y=fun(t)y=exp(-0.5*t).*sin(t+pi/6);例：求积分?d=pi/1000;建立被积函数：27自适应法-simpson法求积lQuadQ=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不为零时给出图形；自适应法-simpson法求积Quad28例：求积分建立被积函数：fun1.mfunction f=fun1(x)f=x./(4+x.2);在MATLAB命令窗口中：?quad(fun1,0,1)ans=0.111572382538912例：求积分建立被积函数：fun1.m在MATLAB命令窗口29自适应法-cotes法求积lQuad8Q=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不为零时给出图形；自适应法-cotes法求积Quad830例：求积分建立被积函数：fun2.mfunction f=fun2(x)f=exp(-x/2);在MATLAB命令窗口：?quad8(fun2,1,3,1e-10,1)ans=0.766800999128407例：求积分建立被积函数：fun2.m在MATLAB命令窗口：31符号积分lSymbolicToolBoxint(S)%计算符号表达式的不定积分；int(S,v)%计算符号表达式对自变量v的不定积分；int(S,a,b)%计算从a到b的定积分；int(S,v,a,b)%计算对变量v从a到b的定积分；符号积分SymbolicToolBox34符号积分?syms x t;?a=cos(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符号积分?syms x t;?syms x135微分和差分l数值微分和差分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,:)l符号微分和差分Diff(S)%对符号表达式S求微分Diff(S,v)%对自变量v求微分Diff(S,n)%n为正整数，对S求n次微分Diff(S,v,n)微分和差分数值微分和差分36?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)?x=sym(x)；?syms x37l多元函数导数lJacobian(f,v)%计算对向量v的Jacobi矩阵多元函数导数38例：?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例：?f2=3*x-cos(x*y)-0.5;x2-839