数值分析与计算方法matlab例程讲解

1、计算方法学号： 专业班级：姓名：ET6V目录1 拉格朗日插值 32 牛顿插值多项式 43 不动点迭代法 64 斯特芬森迭代法 75 牛顿迭代法 86 欧拉法 97 改进欧拉法 118 曲线拟合的最小二乘法 12拉格朗日插值n 次插值基数为：lk x (x x0)(x x1).(x xk 1)(x xk 1).( x xn) ，k=(xk x0)(xk xk 1)(xk xk 1).( xk xn)n0,1,2 n 则我们把形如Ln( x)yk lk ( x) 的插值多项式叫做拉格朗日插值k0多项式当 n=1 时 L1(x) yklk (x) yk 1lk 1(x)当n=2 时L2(x) yk1

2、lk 1(x) yklk(x) yk1lk 1(x)Lagrange.m 文件内容如下function yi= Lagrange(x,y,xi)m = length(x); n=length(y);p=length ( xi);if m=n error( 向量x 与 y 的长度必须一致); ends=0;for k = 1:nt= ones(1,p);for j =1:nif j =kt=t.*(x(i)-x(j)/(x(k)-x(j);ends=s+t*y(k);endyi= s;例题：计算 115在 matlab 的命令窗口执行x=100,121;y=10,11;y1=Lagrangex,

3、y,115得到： y1=10.7143;牛顿插值多项式牛顿插值多项式的表达式如下：Nn (x) c0 c1(x x0)cn (x x0)(x x1) (x xn 1)其中每一项的系数 ci 的表达式如下：xix0ci f(x0,x1, ,xi) f(x1,x2, ,xi) f(x0,x1, ,xi 1)根据ci 以上公式，计算的步骤如下：1、计算 f(x0), f(x1), , f (xn)2、计算 f(x0,x1), , f(xn 1,xn)n、计算 f(x0,x1, ,xn 1), f (x1, ,xn 1,xn) n 1、计算 f(x0,x1, ,xn 1,xn)Newtint.m 文件

4、内容如下：function yi=newtonint(x,y,xi)m=length(x);n=length(y);endif m=n error( 向量x 与 y 的长度必须一致)k=2; f(1)=y(1)while k=n+1f(1)=y(k);k,x(k)for i= 1:k-1if i=k-1f(i+1)=(f(i)-y(i)/(x(k)-x(i);endendcs(i)=f(i+1);y(k)=f(k);k=k+1;endcfwh=0;for i=0:n-2w=1;for j=1:iw=w* (xi-x(j) )；endcfwh=cfwh+cs(i)*w;endyi=y(1)+cf

5、wh;例题：已知函数 f (x) shx 的函数值如下表，构造 4 次 Newton 插值多项式并计算 f (0.596) sh0.596 的值。k012345xk0.400.550.650.800.901.05f(xk )0.410750.578150.696750.888111.026521.25386解：在 matlab 的命令窗口执行 ： x=0.40,0.55,0.65,0.80,0.90,1.05;y=0.41075,0.57815,0.69675,0.88811,1.02652,1.25382;xi=0.596;ni=newtonint(x,y,xhat);得:kxkfxkfx0

6、k fx01kf x012k f x0123k f x01234k00.40000.410810.55000.57821.116020.65000.69671.14400.280030.80000.88811.19340.30960.197340.90001.02651.23150.33010.20050.031251.05001.25391.29710.36220.20550.03250.0085则：N4(x) 0.4108 1.116(x 0.4) 0.2800(x 0.4)(x 0.55)0.1973(x 0.4)(x 0.55)(x 0.65) 0.0312(x 0.4)(x 0.55

7、)(x 0.65)(x 0.8)7sh0.596 N4(0.596) 0.6319三不动点迭代法对于 x (x) 若要求 x* 满足 f(x*) 0,则x* f (x* ) ；反之亦然，称x* 为 f(x )的一个不动点。求 f(x )的零点就等价于求 (x)的不动点，选择一个初 始近似值 x0 ，将它带入 x (x)，即可求得 x1 (x0) ，如此就可反复计算 xk 1 (xk ) ,k=0,1, , (x) 称为迭代函数。如果对于 x0 a,b,由 xk 1 (xk ) ,k=0,1, ,得到序列 xk 有极限 lim xk x* ，则称迭代方程 xk 1 (xk ) ,k=0,1, 收

8、敛，且 x* (x*) 为 ( x)的不动点，此种方法称为不动点 迭代法iterate.m 文件如下functionx_star,k=iterate(fun,x0,ep,Nmax)if nargin4 Nmax=500; endif narginep&k fun =inline( exp(-x) );x_star,k=iterate(fun,0.5);迭代 18 次后得： x_star=0.5671四斯特芬森迭代法将埃特金加速技巧与不动点迭代结合，则可得到如下的迭代法 ：yn 1 yn f(xn ,yn )xn 1 xnk=0 ，1 称之为特芬森迭代法yn 1 yn hf(x n,yn )St

9、effensen.m 文件functionx_star,k=steffensen(fun,x0,ep,Nmax);if nargin4 Nmax=500; endif narginep&k fun =inline( x3-1 );x_star,k=iterate(fun,0.5);得:kxkkxk01.5000000000000041.3289513241715011.5000200000000051.324804592489501021.4163095909551961.3247179940566631.35566023877734五牛顿迭代法设已知方程 f(x) 0有近似根 xk(假定 f

10、(xk) 0 )，将函数 f(x)在点 xk展 开，有 f(x) f(xk) f(xk)(x xk)于是方程 f (x) 0 可近似地表示为f(xk) f(xk)(x xk) 0 。这是个线性方程，记其根为 xk 1，则xk 1的计算公式为 xk 1 xk f(xk) ，k=0,1 这就是牛顿法。f (xk )Newton.m 文件functionx_star,k=Newton(fname,dfname,x0,ep,Nmax)if nargin5 Nmax=500; endif narginep&kfname=inline( x3+2*x2+10*x-20 );dfname=inline( 3

11、*x2+4*x+10 );x_star,k=Newton(fname,dfname,1);迭代 4 次后的计算结果如下表：k1234xk1.411764711.369336741.368808191.36880810六欧拉法一般的在一曲线中，做折线，设折线的顶点 Pn，过 Pn(xn,yn)依方向场的方向yn 1 yn f(xn, yn) 在推进到 Pn1(xn1,yn 1) ，显然两个顶点 Pn ，Pn 1的坐标有关系 xn 1 xnn nyn 1 yn hf (xn,yn)这就是欧拉法。实际上，这是对常微分方程中的导数用均差近似，即12n y(xn) f (xn, y(xn) 直接得到的，

12、若初值 y0已知，则可由 hyn 1 yn f(xn,yn)y1 y0 hf(x 0,y 0)xn 1 xn逐次算出 y2 y1 hf(x 1,y 1)yn 1 yn hf ( xn , yn ) eular1.m 文件如下 functionx,y=eular1(dyfun,xspan,y0,h) x=xspan(1) :h :xspan(2) ; y(1)=y0 ;for n=1 :length(x)-1y(n+1)=y(n)+h*feval(dyfun,x(n),y(n); endx=x;y=y;例题：dy2xy y ,x 0,1 求解初值问题 dyyy(0) 1解：在 matlab 的命

13、令窗口执行输入：dyfun=inline( y-2*x/y );x,y=eular1(dyfun,0,1,1,0.2) ；x,y得到ans= 0 1.0000130.2000 1.20000.4000 1.37330.6000 1.53150.8000 1.6811 1.0000 1.8269七改进欧拉法先用欧拉公式求的一个初步的近似值 yn 1，称之为预测值，预测值 yn 1的 精度可能很差，在用梯形公式 yn 1 yn h(f(xn,yn) f (xn 1,yn 1)将它校正一次， 得到的结果 yn 1称为校正值， 则这样的预测校正系统通常称之为该进的欧拉公式预测yn 1 yn hf (x

14、n,yn)h 或表示为以下平均化形式校正yn 1 yn 2(f(xn,yn) f(xn 1,yn 1)yp yn hf (xn,yn)ye yn hf (xn 1,yp)1yn 1 (yp ye)2eular1.m 文件如下functionx,y=eular2(dyfun,xspan,y0,h)x=xspan(1) :h :xspan(2) ; y(1)=y0 ;for n=1 :length(x)-114k1=feval(dyfun,x(n),y(n);y(n+1)=y(n)+h*k1;k2=feval(dyfun,x(n+1),y(n+1);y(n+1)=y(n)+h*(k1+k2)/2;

15、endx=x;y=y;例题：用改进欧拉法求解初值问题ddy ydy2x,x 0,1 yy(0) 1解：在 matlab 的命令窗口执行输入： dyfun = inline( y-2*x/y );x,y=eularg2(dyfun,0,1,1,0.2);x,y得 ans = 01.00000.2000 1.18670.4000 1.34830.6000 1.49370.8000 1.6279151.00001.7542八曲线拟合的最小二乘法设 0 1, n是 n个线性无关的连续函数， n是由 i(i 0,1,2., n)的所有线性组合构成的函数集合，记 n span 0 1, n，任取 x n(x) ，则n(x)ak k(x) ，对于已知点 (xi,yi) (i 0,1,2., m) ，在 n中求一函数k0m m n(x) ， 使得 | 22|(xi( )yi 2 )(ak i xi( )y i2达到)最小，这i 0 i 0 k 0就是最小二乘拟合 t 问题。ployfit.m 文件function p= ployfit(x,y,m)A=zeros(m+1,m+1) ;for i=0:mfor j=0:mA(i+1,j+1)=sum(x.(i+j);endb(i+1)=sum(x.i.*y);end16a=Ab ;p=fliplar(a );例题：