matlab实现数值分析插值及积分.docx

1、matlab实现数值分析插值及积分Matlab实现数值分析插值及积分摘要：数值分析(numerical analysis)是研究分析用计算机求解数学计算问题的数值计算方法及 其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。 在实际生产实践中， 常常将实际问题转化为数学模型来解决， 这个过程就是数学建模。 学习数值分析这门课程可以让我们学到很多的数学建模方法。分别运用matlab数学软件编程来解决插值问题和数值积分问题。 题目中的要求是计算差值和积分，对于问题一，可以分别利用朗格朗日插值公式，牛顿插值公式， 埃特金逐次线性插值公式来进行编程求解，具体 matla

2、b代码见正文。编程求解出来的结果为： =+。其中Aitken插值计算的结果图如下：对于问题二，可以分别利用复化梯形公式，复化的辛卜生公式，复化的柯特斯公式编 写程序来进行求解，具体 matlab代码见正文。编程求解出来的结果为： 0.6932其中复化梯形公式计算的结果图如下：ator Documents MATLAB 匚Dmmand Window clear t ixing ()用真化梯形算法计算的结果Iz 0. 6932等分埶n= 58A 问题重述问题一：已知列表函数表格i01234121782257分别用拉格朗日，牛顿，埃特金插值方法计算。问题二：用复化的梯形公式，复化的辛卜生公式，复化的

3、柯特斯公式计算积分， 使精度小于5。问题解决问题一：插值方法对于问题一，用三种差值方法：拉格朗日，牛顿，埃特金差值方法来解决。一、拉格朗日插值法：拉格朗日插值多项式如下：首先构造n 1个插值节点XoXj , xn上的n插值基函数，对任一点人所对应的插值基函数 li(x)，由于在所有Xj(j 0,1, ,i 1,i 1, ,n)取零值，因此li(x)有因子(x x0) (x XiJ(x Xi J (x xn)。又因li (x)是一个次数不超过 n的多项式，所以只可能相差一个常数因子，固 lj(x)可表示成：li(x) A(X Xo) (X XiJ(X X 1) (x Xn)利用li(Xi)1得：

4、A 1(Xi X0) (X Xi J(X Xi 1) (Xi Xn)于是(i 0,1,2, ,n)(X Xo) (X Xi 1)(X Xi1) (X Xn)(Xi Xo) (Xi Xi 1)( Xi Xi 1) (Xi Xn)Ln(x) yjl j(x)j0从而 n 次拉格朗日插值多项式为：nLn(xi) yjl j(xi)j0(i 0,1,2, ,n)matlab 编程：编程思想：主要从上述朗格朗日公式入手：依靠循环，运用 poly ()函数和conv ()函数表示拉格朗日公式，其中的 poly( i)函数表示以i作为根的多项式的系数，例如 poly( 1)表示 x-1 的系数，输出为 1

5、-1,而 poly( poly( 1)表示(x-1)*( x-1) =xA2-2*x+1 的系数，输出为1 -2 1;而conv ()表示多项式系数乘积的结果， 例如conv (poly( 1 )，poly( 1)输出为1 -2 1;所以程序最后结果为 xAn+xAn-1 +xA2+x+1( n的值据结果的长度为准)的对应系数。在命令窗口输入 edit lagran 来建立 lagran.m 文件,文件中的程序如下：function c,l=lagran(x,y)w=length(x);n=w-1;l=zeros(w,w);for k=1:n+1v=1;for j=1:n+1if k=jv=c

6、onv(v,poly(x(j)/(x(k) -x(j);endl(k,:)=v;end c=y*l;输入： x=0 1 2 3 4; y=1 2 17 82 257; lagra n( x,y)运行结果为ans =1.0000 -0.0000 -0.0000 0 1.0000结果为：=+。如图表1:iand Wi 1 dow c Lear 3E= O1231j； 尸J 2 17 82 257; y)axis =一CL 0000 -0. 00001. OtlOU图表1.牛顿插值法newt on插值多项式的表达式如下:Nn(x) C0 6(x X。)Cn(x X0)(X X1) (xxn 1 )其

7、中每一项的系数 ci的表达式如下:Ci fXo,X1, ,Xif 知 X2, ,XiXfX0,X1,X0,Xi 1即为f (x)在点X0,X1, ,Xi处的i阶差商，（f Xif (Xi) , i1,2,L , n），由差商fXo,X1,x的性质可知：fXo,X1,i i,Xi f(Xj)1j 0 k 0k jXj xkmatlab 编程：编程思想：主要从上述牛顿插值公式入手：依靠循环，运用 poly （）函数和conv （）函数表示拉格朗日公式，其中的 poly（ i）函数表示以i作为根的多项式的系数，例如 poly（ 1）表示 x-1 的系数，输出为 1 -1,而 poly（ poly（

8、1）表示（x-1）*（ x-1） =xA2-2*x+1 的系数，输出为1 -2 1;而conv （）表示多项式系数乘积的结果， 例如conv （poly（ 1 ），poly（ 1）输出为1 -2 1;所以程序最后结果为 xAn+xAn-1 +xA2+x+1（ n的值据结果的长度为准）的对应系数。在命令窗口输入 edit nowpoly来建立newpoly.m文件，文件中的程序如下: function c,d=n ewpoly(x,y)n=len gth(x);d=zeros (n,n);d(：,1)=y;for j=2: nfor k=j:nd(k,j)=(d(k,j-1卜d(k-1,j-1)

9、/(x(k) -x(k-j+1);end c=d( n,n);for k=( n-1):-1:1c=c on v(c,poly(x(k); m=le ngth(c); c(m)=c(m)+d(k,k);end输入： x=0 1 2 3 4; y=1 2 17 82 257; n ewpoly(x,y)运行结果为ans =1 0 0 0 1所以=+。如图表2:C m m n dl LVLndouvclearCO 1 2 3 41y= LI 2 17 82257Jne wp o l.y kxj.皆JSJTLS =1 0图表2三.埃特金插值法:Aitken插值公式如下:递推表达式为：当 n=1 时，

10、= +当 n=2 时，= +其中的带入递推表达式求得。由此递推下去，最终得到的结果。matlab 编程：编程思想： 埃特金插值多项式又称作 Aitken 逐次线性插值多项式， 根据公式的特点，可以 利用 2 次嵌套循环将公式表示出来。在命令窗口输入 edit Aitken 来建立 Aitken.m 文件，文件中的程序如下：function f = Aitken(x,y)syms z;n = length(x);y1(1:n) = z;for i=1:n -1for j=i+1:ny1(j) = y(j)*(z -x(i)/(x(j) -x(i)+y(i)*(z -x(j)/(x(i) -x(j

11、);endy = y1;simplify(y1);simplify(y1(n);f = collect(y1( n);输入： x=0 1 2 3 4; y=1 2 17 82 257; Aitke n( x,y)运行结果为ans = zA4 + 1所以=+。如图表3:图表3问题二：复化积分对于问题二来说，用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式结局问 题（计算积分，使精度小于 5 ）。复化的梯形公式:复化梯形的迭代公式为:matlab 编程：程序1 (求f (x)的n阶导数：)在命令窗口输入 edit qiudao 来建立 qiudao.m 文件，文件中的程序如下： functio

12、n d=qiudao(x,n)syms x;f=1/x;n=input( 输入导数阶数： );d=diff(f,x,n);输入 :qiudao ( x,n)输入所求导数阶数 :2显示：n = 2ans =2/xA3结果为 ：f2 =2/xA3如图表 4：程序2:在命令窗口输入 edit tixi ng来建立tixing.m文件，文件中的程序如下：function y=tixing()syms x ；%定义自变量xf=inline(1/x,x)；%定义函数f(x)= 1/xf2=i nlin e(2/xA3,x)；%定义f(x)的二阶导数，输入程序1里求出的f2即可。f3=-2/xA3；%因fm

13、inbnd ()函数求的是表达式的最小值，且要求表 达式带引号，故取负号，以便求最大值e=5*10A( -5)；%精度要求值a=1 ；%积分下限b=2 ；%积分上限x1=fminbnd(f3,1,2)；%求负的二阶导数的最小值点，也就是求二阶导数的最大 值点对应的X值for n=2:1000000 ；%求等分数nRn二(b-a)/12*(b-a)/n)A2*f2(x1) ; %计算余项if abs(R n) t i笛金怦它CJ0 6 9 32铮汨颛n= 53A 图表5复化的辛卜生公式：复化 simpson 迭代公式为：matlab 编程：程序1 (求f (x )的n阶导数：)在命令窗口输入 e

14、dit qiudao 来建立 qiudao.m 文件，文件中的程序如下： function d=qiudao(x,n)syms x;f=1/x;n=input( 输入导数阶数： );d=diff(f,x,n);输入 :qiudao ( x,n)输入所求导数阶数 :4显示：n = 4ans =24/xA5结果为 ：f2 = 24/xA5如图表 6：榆入加术吕姜比仍如丄n. 4=am.s =5图表6程序2:在命令窗口输入 edit xinpusheng 来建立xinpusheng.m文件，文件中的程序如下:function y=x in pushe ng()syms x ;f=in li ne(1

15、/x,x);f2=i nlin e(24/xA5,x);f3=-24/xA5;e=5*10A( -5);a=1;b=2;x1=fmi nbn d(f3,1,2);for n=2:1000000Rn二(b-a)/180*(b -a)/(2*n)A4*f2(x1);if abs(R n) clear xinpusheng()运行结果为用Simpso n公式计算的结果 Sn= 0.6932等分数 n= 4结果如图表 7 edit kinpusti clearSn=0- 5932 *用S ZLifcp h on处式计尊的结呆图表7三复化的柯特斯公式:牛顿-柯特斯公式如下：(x )dxAkf (x k

16、)lk(x)dxb (x x)L(x Xki)(X Xki)L (x xn) dx a (Xk Xo)L (Xk Xk i)(Xk Xk i)L (Xk Xn)matlab 编程：(1)在命令窗口输入 edit Newt on Cotes 来建立 Newto nCotes.m文件，文件中的程序如下: fun ctio n y,Ck,Ak=Newton Cotes(fu n, a,b, n)if nargin=1mm,nn=size(fu n);if mm=8error(为了保证NewtonCotes积分的稳定性，最多只能有 9个等距节点！)elseif nn=2error（fun构成应为：第一

17、列为 x,第二列为y,并且个数为小于10的等距节点！）xk=fun(1,:);fk=fun(2,:);a=min(xk);b=max(xk);n=mm -1;elseif nargin=4xk=linspace(a,b,n+1);if isa(fun,function_handle)fx=fun(xk);elseerror(fun 积分函数的句柄，且必须能够接受矢量输入！ )endelseerror( 输入参数错误，请参考函数帮助！ )endCk=cotescoeff(n);Ak=(b-a)*Ck;y=Ak*fx;(2) 在命令窗口输入 edit cotescoeff 来建立 cotescoe

18、ff.m 文件，文件中的程序如下：function Ck=cotescoeff(n)for i=1:n+1k=i-1;Ck(i)=(-1)A(n-k)/factorial(k)/factorial(n -k)/n*quadl(t)intfun(t,n,k),O,n);3)在命令窗口输入 edit intfun 来建立 intfun.m 文件，文件中的程序如下：function f=intfun(t,n,k)f=1;for i=0:k -1,k+1:nf=f.*(t -i);end% fun ，积分表达式，这里有两种选择%(1) 积分函数句柄，必须能够接受矢量输入，比如 fun=(x)1./x%

19、 (2)x,y 坐标的离散点， 第一列为 x ， 第二列为 y ， 必须等距， 且节点的个数小于 9， 比如： fun=1:8;1./(1:8)% 如果 fun 的表采用第二种方式，那么只需要输入第一个参数即可，否则还要输入 a,b,n 三个参数% a,积分下限% b,积分上限% n,牛顿-科特斯数公式的阶数，必须满足 1n=8时不能保证公式的稳定性% (1)n=1 ,即梯形公式% (2)n=2 ,即辛普森公式% (3)n=4 ,即科特斯公式在显示窗口输入：fun=(x)1./x;a=-1;b=1;n=4;Newt on Cotes(fu n, a,b ,n)运行结果为：ans =0.6932结果如图表8 f i, :a 1 ub-2n=4 :Umrt 口itC 口 1: e-x (-Eunj a. bj n) ans 0a d9J2图表8