1、 Matlab 实现数值分析插值及积分摘要：数值分析(numerical analysis)是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。在实际生产实践中，常常将实际问题转化为数学模型来解决，这个过程就是数学建模。学习数值分析这门课程可以让我们学到很多的数学建模方法。分别运用matlab数学软件编程来解决插值问题和数值积分问题。题目中的要求是计算差值和积分，对于问题一，可以分别利用朗格朗日插值公式，牛顿插值公式，埃特金逐次线性插值公式来进行编程求解，具体matlab代码见正文。编程求解出来的结果为：=+。其中
2、Aitken插值计算的结果图如下：对于问题二，可以分别利用复化梯形公式，复化的辛卜生公式，复化的柯特斯公式编写程序来进行求解，具体matlab代码见正文。编程求解出来的结果为： 0.6932其中复化梯形公式计算的结果图如下：问题重述问题一：已知列表函数表格 10 1 2 3 41 2 17 82 257分别用拉格朗日，牛顿，埃特金插值方法计算。问题二：用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式计算积分，使精度小于 5。问题解决问题一：插值方法对于问题一，用三种差值方法：拉格朗日，牛顿，埃特金差值方法来解决。一、拉格朗日插值法：拉格朗日插值多项式如下： 首先构造 个插值节点 上的 插值
3、基函数，对任一点 所对应的插值基函1nnx,10 ix数 ，由于在所有 取零值，因此 有因子)(xli ),(ijx )(li。又因 是一个次数不超过 的多项式，所以)110 nii x (xli n只可能相差一个常数因子，固 可表示成：(li )()() 110 niii xxxAl 利用 得：1)(ixl )()()(110 niiii xxxA于是 ),210()()()() 110 nixxxxl niiiii 因此满足 的插值多项式可表示为：iinyxL)(),210(nnjjxlyxL0)()(从而 次拉格朗日插值多项式为：n ),210()()(0 nixlyxLnjijin m
4、atlab 编程：编程思想：主要从上述朗格朗日公式入手：依靠循环，运用 poly（）函数和 conv（）函数表示拉格朗日公式，其中的 poly（i）函数表示以 i 作为根的多项式的系数，例如 poly（1）表示 x-1 的系数，输出为 1 -1，而 poly（poly （1） ）表示（ x-1）*（x-1 ）=x2-2*x+1 的系数，输出为 1 -2 1；而 conv（）表示多项式系数乘积的结果，例如 conv（poly（1 ） ，poly（1） ）输出为 1 -2 1；所以程序最后结果为 xn+xn-1+x2+x+1（n 的值据结果的长度为准）的对应系数。在命令窗口输入 edit lagr
5、an 来建立 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=conv(v,poly(x(j)/(x(k)-x(j);endendl(k,:)=v;endc=y*l;输入： x=0 1 2 3 4; y=1 2 17 82 257; lagran(x,y)运行结果为ans =1.0000 -0.0000 -0.0000 0 1.0000结果为：=+。如图表1：图表 1二牛顿插值法newton 插值多项式的表达式如下： 010
6、011()()()()n n nNxcxcxx其中每一项的系数 ci 的表达式如下：120101,i ii ifxfxcfx 即为 f (x)在 点 处的 i 阶差商， （ ， ） ，由差商01,i ()iifxf,2n的性质可知：01,i 01001,()iiijjkjkfxfxx matlab 编程：编程思想：主要从上述牛顿插值公式入手：依靠循环，运用 poly（）函数和 conv（）函数表示拉格朗日公式，其中的 poly（i）函数表示以 i 作为根的多项式的系数，例如 poly（1）表示 x-1 的系数，输出为 1 -1，而 poly（poly （1） ）表示（ x-1）*（x-1 ）=
7、x2-2*x+1 的系数，输出为 1 -2 1；而 conv（）表示多项式系数乘积的结果，例如 conv（poly（1 ） ，poly（1） ）输出为 1 -2 1；所以程序最后结果为 xn+xn-1+x2+x+1（n 的值据结果的长度为准）的对应系数。在命令窗口输入 edit nowpoly 来建立 newpoly.m 文件，文件中的程序如下：function c,d=newpoly(x,y)n=length(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)/(x(k)-x(k-j+1);endend
8、c=d(n,n);for k=(n-1):-1:1c=conv(c,poly(x(k);m=length(c);c(m)=c(m)+d(k,k);end输入： x=0 1 2 3 4; y=1 2 17 82 257; newpoly(x,y)运行结果为ans =1 0 0 0 1所以=+。如图表2：图表 2三埃特金插值法：Aitken 插值公式如下：递推表达式为：= + 当 n=1 时，= + 当 n=2 时，= + 其中的带入递推表达式求得。由此递推下去，最终得到的结果。matlab 编程：编程思想：埃特金插值多项式又称作 Aitken 逐次线性插值多项式， 根据公式的特点，可以利用 2
9、次嵌套循环将公式表示出来。在命令窗口输入 edit Aitken 来建立 Aitken.m 文件，文件中的程序如下：function f = Aitken(x,y)syms z;n = length(x); y1(1:n) = z; for i=1:n-1 for j=i+1:ny1(j) = y(j)*(z-x(i)/(x(j)-x(i)+y(i)*(z-x(j)/(x(i)-x(j); endy = y1;simplify(y1);endsimplify(y1(n); f = collect(y1(n); 输入： x=0 1 2 3 4; y=1 2 17 82 257; Aitken(x
10、,y)运行结果为ans =z4 + 1所以=+。如图表 3：图表 3问题二：复化积分对于问题二来说，用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式结局问题（计算积分，使精度小于 5） 。一 复化的梯形公式：复化梯形的迭代公式为：;matlab 编程：程序 1（求 f（x）的 n 阶导数：）在命令窗口输入 edit qiudao 来建立 qiudao.m 文件，文件中的程序如下：function d=qiudao(x,n)syms x;f=1/x; n=input(输入导数阶数： ); d=diff(f,x,n);输入:qiudao （x,n）输入所求导数阶数:2显示：n = 2ans =
11、2/x3结果为：f2 =2/x3如图表4：图表 4程序 2：在命令窗口输入 edit tixing 来建立 tixing.m 文件，文件中的程序如下：function y=tixing()syms x ； %定义自变量 xf=inline(1/x,x)； %定义函数 f(x)= 1/x f2=inline(2/x3,x) ； %定义 f(x)的二阶导数,输入程序 1 里求出的 f2 即可。f3=-2/x3； %因 fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，以便求最大值e=5*10(-5)； %精度要求值 a=1； %积分下限b=2； %积分上限x1=fminbn
12、d(f3,1,2)； %求负的二阶导数的最小值点，也就是求二阶导数的最大值点对应的 x 值for n=2:1000000； %求等分数 nRn=-(b-a)/12*(b-a)/n)2*f2(x1)； %计算余项if abs(Rn) clear xinpusheng()运行结果为用Simpson公式计算的结果 Sn= 0.6932等分数 n= 4结果如图表 7图表 7三 复化的柯特斯公式：牛顿-柯特斯公式如下：matlab 编程：(1)在命令窗口输入 edit NewtonCotes 来建立 NewtonCotes.m 文件，文件中的程序如下：function y,Ck,Ak=NewtonCot
13、es(fun,a,b,n) if nargin=1 mm,nn=size(fun);if mm=8 error(为了保证 NewtonCotes 积分的稳定性，最多只能有 9 个等距节点！) elseif nn=2 error(fun 构成应为：第一列为 x，第二列为 y，并且个数为小于 10 的等距节点！) end 0()()nbkafxdAfx011()()()()()bbkknkaak kxAlxd dxx xk=fun(1,:); fk=fun(2,:); a=min(xk); b=max(xk); n=mm-1; elseif nargin=4 xk=linspace(a,b,n+1
14、); if isa(fun,function_handle) fx=fun(xk); else error(fun 积分函数的句柄，且必须能够接受矢量输入！) end elseerror(输入参数错误，请参考函数帮助！) end Ck=cotescoeff(n); Ak=(b-a)*Ck; y=Ak*fx; （2）在命令窗口输入 edit cotescoeff 来建立 cotescoeff.m 文件，文件中的程序如下：function Ck=cotescoeff(n) for i=1:n+1 k=i-1; Ck(i)=(-1)(n-k)/factorial(k)/factorial(n-k)/
15、n*quadl(t)intfun(t,n,k),0,n); end （3 ）在命令窗口输入 edit intfun 来建立 intfun.m 文件，文件中的程序如下：function f=intfun(t,n,k) f=1; for i=0:k-1,k+1:n f=f.*(t-i); end% fun，积分表达式，这里有两种选择 %(1)积分函数句柄，必须能够接受矢量输入，比如 fun=(x)1./x % (2)x,y 坐标的离散点， 第一列为 x， 第二列为 y， 必须等距， 且节点的个数小于 9， 比如： fun=1:8;1./(1:8)% 如果 fun 的表采用第二种方式，那么只需要输入第一个参数即可，否则还要输入 a,b,n三个参数 % a，积分下限 % b，积分上限 % n，牛顿-科特斯数公式的阶数，必须满足 1=8 时不能保证公式的稳定性 % (1)n=1，即梯形公式 % (2)n=2，即辛普森公式% (3)n=4，即科特斯公式 在显示窗口输入：fun=(x)1./x;a=-1;b=1;n=4;NewtonCotes(fun,a,b,n) 运行结果为：ans =0.6932结果如图表 8图表 8
