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

1、MATLAB在计算方法中的应用MATLAB入门到精通插值与拟合l插值与拟合来源于实际，又广泛应用于实际。随着计算机的不断发展及计算水平的不断提高，他们在国民生产和科学研究等方面扮演着越来越重要的角色。l插值法主要有Lagrange插值、分段线性插值、Hermite插值及三次样条插值等l拟合法主要有最小二乘法拟合和快速Fourier变换等Lagrange插值l 对给定的n个插值节点及相应的函数值，利用n次Lagrange插值多项式对插值区间内任意x的函数值y可通过下式求得：nkjkjnkjjkxxxxyxy11)()(MATLAB实现function y=lagrange(x0,y0,x)% l

2、agrange 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;end例：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.223

3、144; lagrange(x,y,0.54)ans = -0.61614 log(0.54)ans = -0.61619Runge现象l 19世纪Runge给出了一个等距节点插值多项式不收敛的例子： 在区间-5，5上各阶导数存在，但在此区间上取n个节点构造的Lagrange插值多项式在全区间上不收敛。211)(xxfRunge现象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

4、 plot(x0,y1,-b)-5-4-3-2-1012345-0.500.511.52分段线性插值l 分段线性插值就是通过插值点用折线段连接起来逼近原曲线。l MATLAB实现 interp1 一维插值 yi=interp1(x,y,xi) %对一组节点(x,y)进行插值，计算插值点xi的函数值。 yi=interp1(y,xi) %默认x=1:n yi=interp1(x,y,xi,method) %method为指定插值算法l nearest 线性最近项插值l linear 线性插值l spline 三次样条插值l cubic 三次插值 同类的函数还有inter1q,interpft,s

5、pline,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) hold on plot(x0,y1,-b) plot(x0,y2,*m)-505-0.500.511.52Hermite插值l 不少问题不

6、但要求在节点上函数值相等，而且要求导数值也相等，甚至要求高阶导数也相等，可用Hermite插值多项式。l 已知n个插值节点及其对应的函数值和一阶导数值，则计算插值区域内任意x的函数值y的Hermite插值公式为：nijijijijnijjiniiiiiiixxaxxxxhyyyaxxhxy12111;)( )2)()(其中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:n h=1.0;a=0.0; for j=1:n if j=i h=

7、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.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=hermit

8、e(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.30.3050.310.3150.320.3250.330.3350.340.3450.350.2950.30.3050.310.3150.320.3250.330.3350.340.345三次样条插值l 样条函数可以给出平滑的插值曲线和曲面。l 方法介绍：设在a,b上给定n+1个点 a=x0 x1x2x=28 28.7 29 30;y=4.1 4.3 4.1 3.0;x0=

9、28:0.15:30;y1=spline(x,y,x0);plot(x,y,-r)hold onplot(x0,y1)2828.228.428.628.82929.229.429.629.83033.544.5最小二乘法拟合l在科学实验的统计方法中，往往要从一组实验数据(xi,yi)中寻找出自变量x和因变量y之间的函数关系y=f(x)。l由于观测数据往往不够准确，因此并不要求y=f(x)经过所有的点(xi,yi)，而只要求在给定点xi上误差ei=f(xi)-yi按照某种标准达到最小，通常采用欧氏范数|e|2作为误差度量的标准，这就是最小二乘法。MATLAB实现l 利用polyfit函数进行多项

10、式拟合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 predicti

11、ons. 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 (n

12、ormr) as fields. l 利用矩阵除法解决复杂型函数的拟合例：拟合以下数据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)0.511.522.530246810例：根据经验公式y=a+bx

13、2，拟合如下数据：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)15202530354045020406080100快速Fourier变换l在函数逼近中，当函数为周期函数时，用三角多项式逼近比用代数多项式更合适，因此引入Fourier逼近和快速Fourier变换。lMATL

14、AB实现 fft 快速Fourier变换l fft(x)l fft(x,n)l fft(x,DIM) 或 fft(x,n,DIM) 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逆变换积分与微分l Newton-Cotes系列数值求积将积分区间a,b划分为n等分，步长h=(b-a)/n，选取等距节点xk=a+kh，则求积分公式：Ck(n)由Cotes系数表

15、给出。ln=0 时为矩形公式；ln=1时为梯形公式;ln=2时为Simpson公式；ln=3时为Cotes公式；lsxfCabInkknkn0)()()(l矩形求积公式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 7l梯形求积分公式lZ=trapz(Y) 用梯形方法计算Y的积分近似值。对向量Y，函数返回Y的积分；

16、对矩阵Y，函数返回一个向量，向量各元素为矩阵Y对应列向量的积分值。lZ=trapz(X,Y) lZ=trapz(X,Y,DIM)lZ=trapz(Y,DIM) 例：求积分305 . 0)6/sin(dttet?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);自适应法-simp

17、son法求积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不为零时给出图形；例：求积分dxxx1024建立被积函数： fun1.mfunction f=fun1(x)f=x./(4+x.2);在MATLAB命令窗口中：?quad(fun1,0,1)ans = 0.111572382538912自适应法-cotes法求积lQuad8Q=quad8(F,a,b) 求函数F(x)从a到b的相对误差为1e-3的积分近似值；Q=q

18、uad8(F,a,b,tol) 向量tol给出相对误差和绝对误差；Q=quad8(F,a,b,tol,TRACE) TRACE不为零时给出图形；例：求积分dxex312建立被积函数：fun2.mfunction f=fun2(x)f=exp(-x/2);在MATLAB命令窗口：?quad8(fun2,1,3,1e-10,1)ans = 0.76680099912840711.522.530.20.250.30.350.40.450.50.550.60.65Gauss求积公式l 为了得到较高的代数精度，可以使用Gauss公式： 式中Ak为Gauss-Legender系数l 对任意区间a,b的积分

19、，可变换为如下公式：110)()(nkkkxfAdxxf11)22(2)(dtbatabfabdxxfbaRomberg求积公式l梯形公式积分的算法简单、编程容易，但精度较差、收敛速度较慢；lRomberg求积公式可自动改变积分步长，使相邻的两次值的绝对误差或相对误差小于预先设定的允许误差。符号积分lSymbolicToolBoxint(S) %计算符号表达式的不定积分；int(S,v) %计算符号表达式对自变量v的不定积分；int(S,a,b) %计算从a到b的定积分；int(S,v,a,b) %计算对变量v从a到b的定积分；符号积分?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?