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

计算方法学号：专业班级ET6V姓名：目录1 拉格朗日插值 32 牛顿插值多项式 43 不动点迭代法 64 斯特芬森迭代法 75 牛顿迭代法 86 欧拉法 97 改进欧拉法 11128 曲线拟合的最小二乘法一. 拉格朗日插值n 次插值基数为：1 (x )=(x x 0)(x 心x(x X+( x xn)，k= k(x 一 x ).( x 一 x )(x 一 x ).( x 一 x )k0kk -1 kk +1kn0,l,2.n则我们把形如 l (x) = y 1(x)的插值多项式叫做拉格朗日插值 nk kk 二 0多项式当 n=l 时 l (x) = y 1 (x) + y 1 (x)1 k kk + 1 k +1当 n=2 时 L ( x) = y 1 ( x) + y 1 ( x) + y 1 ( x)2 k 一 1 k 一 1k kk + 1 k + 1Lagrange.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;例题：计算胃肓在 matlab 的命令窗口执行x=100,121;y=10,11;y1=Lagrangex,y,115得到： y1=10.7143;二牛顿插值多项式牛顿插值多项式的表达式如下：N (x) = c + c (x - x ) +b c (x - x )(x - x )n 0 1 0 n 0 1其中每一项的系数Cj的表达式如下：(x - x )n - 1ci = f (x0x 1,xi) =f (x , x，x ) f (x , x，x ) 12i01t1-x - xt0根据c以上公式，计算的步骤如下:t1、计算 f (x ), f (x ), f (x )0 1 n2、计算 f (x , x ), , f (x , x )0 1n - 1 nn、计算 f (x , x , , x ), f (x , , x , x ) 0 1n-11n -1 nn + 1、计算 f (x , x , , x , x )0 1n -1 nNewtint.m 文件内容如下：funCtion yi=newtonint(x,y,xi) m=length(x);n=length(y);if m =n error( 向量x与y的长度必须一致)；end 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) )； end cfwh=cfwh+cs(i)*w; endyi=y(1)+cfwh;例题：已知函数f（x）= shx的函数值如下表，构造4次Newton插值多项式并计算f (0.596) = sh0.596 的值。k012345xk0.400.550.650.800.901.05f (x )k0.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);得:kxkf x kf x 0 kf x 01 kf x 012 kf x 0123 kf x01234 k00.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则：N 4( 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)( x 0.65)( x 0.8)sh 0.596 = N (0.596) 0.63194三不动点迭代法对于x = p( x ) 若要求x*满足 f (x*) = 0 ,则x* = f (x*)；反之亦然，称x* 为f (x )的一个不动点。求f (x )的零点就等价于求9(x)的不动点，选择一个初 始近似值x，将它带入x = ( x )，即可求得x =9 (x )，如此就可反复计算0 1 0x =9 (x ) ,k=0,l,，9 (x)称为迭代函数。如果对于x G a,b,由k +1k0x =9 (x ) ,k=0,l,，得到序列x 有极限lim x = x *，则称迭代方程k+1kkk 8 kx = 9 (x ) ,k=0,l,收敛，且x* = 9 (x*)为9 (x)的不动点，此种方法称为不动点 k +1k迭代法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四斯特芬森迭代法将埃特金加速技巧与不动点迭代结合，则可得到如下的迭代法 :y y、 x1 - x f (Xn yn )k=0，1 称之为特芬森迭代法n + 1ny = y hf (x , y )n+1nn nSteffensen.m 文件functionx_star,k=steffensen(fun,x0,ep,Nmax);if nargin4 Nmax=500; endif narginep&k fun =inline( xW );x_star,k=iterate(fun,0.5);得:kxkkxk01.5000000000000041.3289513241715011.5000200000000051.3248045924895021.4163095909551961.3247179940566631.35566023877734五牛顿迭代法设已知方程f (x) = 0有近似根x (假定f (x 2 0)，将函数f (X)在点x展k k k开，有f (x) q f (x ) + f(x )(x - x ) 于是万程f (x) = 0可近似地表示为k k kf (x ) + f(x )(x - x ) = 0。这是个线性方程，记其根为x ，则x的计算公式为kkkk +1k +1 f (x ) k f(x ) kk=0,l.这就是牛顿法。Newton.m 文件functionx_star,k=Newton(fname,dfname,x0,ep,Nmax)if nargin5 Nmax=500; endif narginep&kfname=inline( xA3+2*xA2+10*x-20);dfname=inline( 3*乂八2+4*乂+10);x_star,k=Newton(fname,dfname,1);迭代 4 次后的计算结果如下表：k1234xk1.411764711.369336741.368808191.36880810六欧拉法一般的在一曲线中，做折线，设折线的顶点P ,过P (x y )依方向场的方向n n n , n在推进到 P ( x y ) ，显然两个顶点 P ， P 的坐标有关系n + 1 n +1, n +1nn +1y 一 y n +1n x 一 xn+1n=f (x , y )nny = y 一 hf (x , y )n +1nn n这就是欧拉法。实际上，这是对常微分方程中的导数用均差近似，即y( xn+1)一 y (叮.y(x ) = f (x , y (x )直接得到的，若初值y已知，则可由hnn n0y 一 y 一 、y = y +hf (x , y )f斤=f (x,y )逐次算出 1 000 x 一 xn n 逐次算出y = y + hf (x , y )n+1n211 1y = y 一 hf (x , y )n+1nn neular1.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;例题：求解初值问题d2 xy = y -, x e 0,1dyfun=inline( y-2*x/y);x,y=eular1(dyfun,0,1,1,0.2)；x,y得到 ans= 01.00000.20001.20000.40001.3733 0.60001.53151.0000 1.82690.8000 1.6811七改进欧拉法先用欧拉公式求的一个初步的近似值y ,称之为预测值，预测值y的 n + 1n +1精度可能很差，在用梯形公式y = y + (f (X , y ) + f (x , y )将它校正一次, n + 1 n 2 n nn +1 n + 1得到的结果 y 称为校正值，则这样的预测校正系统通常称之为该进的欧拉公式 n+1预测yn +1校正yn +1=yn=yn+ f ( x , y )nn+(f (X , y ) + f (X, y )2n nn+1n+1或表示为以下平均化形式y = y + f (x , y )p nn n y = y + f (x , y )e nn +1 p1y = (y + y )n +12 peeular1.m 文件如下functionx,y=eular2(dyfun,xspan,y0,h)x=xspan(1) :h :xspan(2) ; y(1)=y0 ;for n=1 :length(x)-1k1=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;endx=x;y=y;d2 x例题：用改进欧拉法求解初值问题r = y ,x G 0,1 dyfun = inline(y-2*x/y);x,y=eularg2(dyfun,0,1,1,0.2);x,y得 ans = 01.00000.20001.18670.40001.34830.60001.49370.80001.62791.00001.7542八曲线拟合的最小二乘法设q g q是n个线性无关的连续函数，9是由申(i = 0,1,2., n)的所有线 01, nni性组合构成的函数集合，记q =span q ,q . q ， 任取q g q (x)则n01, nx nq(x) = a q (x)，对于已知点(x , y ) (i = 0,1, 2., m)，在q 中求一函数k ki ink=0q (x)，使得 II q ll2=为(q (x ) - y )2 = 2 (为 a q (x ) - y )2 达到最小，这2iik i iii=0i=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+l,j+l)=sum(x/(i+j);endb(i+1)=sum(x.Ai.*y);enda=Ab;p=fliplar(a);例题：已知一组实验数据，求它的拟合曲线i12345xi165123150123141yi187126172125148输入 x=165 123 150 123 141;y=187 126 172 125 148;ployfit(x,y,1)ployfit(x,y,2)得 ans = 1.5138-60.939ans = -0.00222.1326 -104.4234即所求的一次拟合曲线为y = 1.5138 x 60.939二次拟合曲线为y = 0.0022 x 2 + 2.1326 x 104.4234