Matlab实现数值分析插值及积分摘要：数值分析(numerical analysis)是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。

在实际生产实践中，常常将实际问题转化为数学模型来解决，这个过程就是数学建模。

学习数值分析这门课程可以让我们学到很多的数学建模方法。

分别运用matlab数学软件编程来解决插值问题和数值积分问题。

题目中的要求是计算差值和积分，对于问题一，可以分别利用朗格朗日插值公式，牛顿插值公式，埃特金逐次线性插值公式来进行编程求解，具体matlab代码见正文。

编程求解出来的结果为：=+。

其中Aitken插值计算的结果图如下：对于问题二，可以分别利用复化梯形公式，复化的辛卜生公式，复化的柯特斯公式编写程序来进行求解，具体matlab代码见正文。

编程求解出来的结果为： 0.6932其中复化梯形公式计算的结果图如下：问题重述问题一：已知列表函数表格 1分别用拉格朗日，牛顿，埃特金插值方法计算。

问题二：用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式计算积分，使精度小于5。

问题解决问题一：插值方法对于问题一，用三种差值方法：拉格朗日，牛顿，埃特金差值方法来解决。

一、拉格朗日插值法：拉格朗日插值多项式如下：首先构造1+n 个插值节点n x x x ,,,10 上的n 插值基函数，对任一点i x 所对应的插值基函数)(x l i ，由于在所有),,1,1,,1,0(n i i j x j +-=取零值，因此)(x l i 有因子)())(()(110n i i x x x x x x x x ----+- 。

又因)(x l i 是一个次数不超过n 的多项式，所以只可能相差一个常数因子，固)(x l i 可表示成：)())(()()(110n i i i x x x x x x x x A x l ----=+-利用1)(=i i x l 得：)())(()(1110n i i i i i i x x x x x x x x A ----=+-于是),,2,1,0()())(()()())(()()(110110n i x x x x x x x x x x x x x x x x x l n i i i i i i n i i i =--------=+-+-因此满足i i n y x L =)( ),2,1,0(n i =的插值多项式可表示为：∑==nj j j n x l y x L 0)()(从而n 次拉格朗日插值多项式为：),,2,1,0()()(0n i x l y x L nj i j j i n ==∑=matlab 编程：编程思想：主要从上述朗格朗日公式入手：依靠循环，运用poly （）函数和conv （）函数表示拉格朗日公式，其中的poly （i ）函数表示以i 作为根的多项式的系数，例如poly （1）表示x-1的系数，输出为1 -1，而poly （poly （1））表示（x-1）*（x-1）=x^2-2*x+1的系数，输出为1 -2 1；而conv （）表示多项式系数乘积的结果，例如conv （poly （1），poly （1））输出为1 -2 1；所以程序最后结果为x^n+x^n-1+……+x^2+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+1 v=1;for j=1:n+1 if k~=jv=conv(v,poly(x(j)))/(x(k)-x(j)); end endl(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 插值多项式的表达式如下：010011()()()()()n n n N x c c x x c x x x x x x -=+-+⋅⋅⋅+--⋅⋅⋅-其中每一项的系数ci 的表达式如下：12011010[,,,][,,,][,,,]i i i i i f x x x f x x x c f x x x x x -⋅⋅⋅-⋅⋅⋅=⋅⋅⋅=-即为 f (x)在点01,,,i x x x ⋅⋅⋅处的i 阶差商，（[]()i i f x f x =，1,2,,i n =），由差商01[,,,]i f x x x ⋅⋅⋅的性质可知：()0101[,,,]()ii i j j k j k k jf x x x f x x x ==≠⋅⋅⋅=-∑∏matlab 编程：编程思想：主要从上述牛顿插值公式入手：依靠循环，运用poly （）函数和conv （）函数表示拉格朗日公式，其中的poly （i ）函数表示以i 作为根的多项式的系数，例如poly （1）表示x-1的系数，输出为1 -1，而poly（poly（1））表示（x-1）*（x-1）=x^2-2*x+1的系数，输出为1 -2 1；而conv（）表示多项式系数乘积的结果，例如conv（poly（1），poly（1））输出为1 -2 1；所以程序最后结果为x^n+x^n-1+……+x^2+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));endendc=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次嵌套循环将公式表示出来。

在命令窗口输入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));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,y)运行结果为ans =z^4 + 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 =2/x^3结果为：f2 =2/x^3如图表4：图表 4程序2：在命令窗口输入edit tixing 来建立tixing.m文件，文件中的程序如下：function y=tixing()syms x ； %定义自变量xf=inline('1/x','x')；%定义函数f(x)= 1/xf2=inline('2/x^3','x') ；%定义f(x)的二阶导数,输入程序1里求出的f2即可。

f3='-2/x^3'；%因fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，以便求最大值e=5*10^(-5)； %精度要求值a=1；%积分下限b=2；%积分上限x1=fminbnd(f3,1,2)； %求负的二阶导数的最小值点，也就是求二阶导数的最大值点对应的x值for n=2:1000000；%求等分数nRn=-(b-a)/12*((b-a)/n)^2*f2(x1)； %计算余项if abs(Rn)<e %用余项进行判断break % 符合要求时结束endendh=(b-a)/n； %求hTn1=0；for k=1:n-1 %求连加和xk=a+k*h；Tn1=Tn1+f(xk)；endTn=h/2*((f(a)+2*Tn1+f(b)))；fprintf('用复化梯形算法计算的结果 Tn=')disp(Tn)fprintf('等分数 n=')disp(n) %输出等分数输入：tixing()运行结果为:用复化梯形计算的结果 Tn= 0.6932等分数 n= 58结果如图表:5图表 5二复化的辛卜生公式：复化simpson迭代公式为：;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）输入所求导数阶数:4显示：n = 4ans =24/x^5结果为：f2 = 24/x^5如图表6：图表 6程序2：在命令窗口输入edit xinpusheng 来建立xinpusheng.m文件，文件中的程序如下：function y=xinpusheng()syms x ;f=inline('1/x','x');f2=inline('24/x^5','x');f3='-24/x^5';e=5*10^(-5);a=1;b=2;x1=fminbnd(f3,1,2);for n=2:1000000Rn=-(b-a)/180*((b-a)/(2*n))^4*f2(x1);if abs(Rn)<ebreakendendh=(b-a)/n;Sn1=0;Sn2=0;for k=0:n-1xk=a+k*h;xk1=xk+h/2;Sn1=Sn1+f(xk1);Sn2=Sn2+f(xk);endSn=h/6*(f(a)+4*Sn1+2*(Sn2-f(a))+f(b));fprintf('用Simpson公式计算的结果 Sn=')disp(Sn)fprintf('等分数 n=')disp(n)调用xinpusheng.m中xinpusheng函数：输入：>> clear>> xinpusheng()运行结果为用Simpson 公式计算的结果 Sn= 0.6932等分数 n= 4结果如图表7图表 7三 复化的柯特斯公式：牛顿-柯特斯公式如下：matlab 编程：(1)在命令窗口输入edit NewtonCotes 来建立NewtonCotes.m 文件，文件中的程序如下： function [y,Ck,Ak]=NewtonCotes(fun,a,b,n)if nargin==1[mm,nn]=size(fun);if mm>=8error('为了保证NewtonCotes 积分的稳定性，最多只能有9个等距节点！') elseif nn~=2error('fun 构成应为：第一列为x ，第二列为y ，并且个数为小于10的等距节点！')endxk=fun(1,:);fk=fun(2,:);a=min(xk);0()()n b k k a k f x dx A f x =≈∑⎰011011()()()()()()()()()b b k k n k k a a k k k k k k n x x x x x x x x A l x dx dx x x x x x x x x -+-+----==----⎰⎰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来建立cotescoeff.m文件，文件中的程序如下：function Ck=cotescoeff(n)for i=1:n+1k=i-1;Ck(i)=(-1)^(n-k)/factorial(k)/factorial(n-k)/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<n<7，因为n>=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。