《matlab实现数值分析报告插值及积分(共19页).doc》由会员分享，可在线阅读，更多相关《matlab实现数值分析报告插值及积分(共19页).doc（19页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上 Matlab实现数值分析插值及积分摘要： (numerical analysis)是研究分析用求解计算问题的及其理论的学科，是的一个分支，它以求解数学问题的理论和方法为研究对象。在实际生产实践中，常常将实际问题转化为数学模型来解决，这个过程就是数学建模。学习数值分析这门课程可以让我们学到很多的数学建模方法。分别运用matlab数学软件编程来解决插值问题和数值积分问题。题目中的要求是计算差值和积分，对于问题一，可以分别利用朗格朗日插值公式，牛顿插值公式，埃特金逐次线性插值公式来进行编程求解，具体matlab代码见正文。编程求解出来的结果为：=+。其中Aitken插值计
2、算的结果图如下： 对于问题二，可以分别利用复化梯形公式，复化的辛卜生公式，复化的柯特斯公式编写程序来进行求解，具体matlab代码见正文。编程求解出来的结果为： 0.6932其中复化梯形公式计算的结果图如下： 问题重述问题一：已知列表函数表格 101234121782257分别用拉格朗日，牛顿，埃特金插值方法计算。问题二：用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式计算积分，使精度小于5。问题解决问题一：插值方法对于问题一，用三种差值方法：拉格朗日，牛顿，埃特金差值方法来解决。一、拉格朗日插值法：拉格朗日插值多项式如下：首先构造个插值节点上的插值基函数，对任一点所对应的插值基函数，由于
3、在所有取零值，因此有因子。又因是一个次数不超过的多项式，所以只可能相差一个常数因子，固可表示成：利用得：于是因此满足 的插值多项式可表示为：从而次拉格朗日插值多项式为：matlab编程：编程思想：主要从上述朗格朗日公式入手：依靠循环，运用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；所
4、以程序最后结果为xn+xn-1+x2+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=j v=conv(v,poly(x(j)/(x(k)-x(j); end end l(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
5、-0.0000 0 1.0000结果为：=+。如图表1：图表 1二牛顿插值法newton插值多项式的表达式如下：其中每一项的系数ci的表达式如下：即为f (x)在点处的i阶差商，（，），由差商的性质可知：matlab编程：编程思想：主要从上述牛顿插值公式入手：依靠循环，运用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），
6、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:n for k=j:n d(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:1 c=conv(c,poly(x(k); m=length(c); c(m)=c(m)+
7、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 = Aitk
8、en(x,y) syms z; n = length(x); y1(1:n) = z; for i=1:n-1 for j=i+1:n y1(j) = y(j)*(z-x(i)/(x(j)-x(i)+y(i)*(z-x(j)/(x(i)-x(j); end y = y1; simplify(y1);end simplify(y1(n); f = collect(y1(n); 输入： x=0 1 2 3 4; y=1 2 17 82 257; Aitken(x,y)运行结果为ans = z4 + 1所以=+。如图表3：图表 3问题二：复化积分对于问题二来说，用复化的梯形公式，复化的辛卜生公式，复
9、化的柯特斯公式结局问题（计算积分，使精度小于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/x3结果为：f2 =2/x3如图表4：图表 4程序2：在命令窗口输入edit tixing 来建立tixing.m文件，文件中的程序如下：function
10、 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=fminbnd(f3,1,2)； %求负的二阶导数的最小值点，也就是求二阶导数的最大值点对应的x值for n=2:； %求等分数n Rn=-(b-a)/12*(b-a)/n)2*f2(x1)；
11、%计算余项 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（求
12、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/x5结果为：f2 = 24/x5如图表6：图表 6程序2：在命令窗口输入edit xinpusheng 来建立xinpusheng.m文件，文件中的程序如下：function y=xinpusheng()syms x ; f=inline(1/x,x); f2=inline(
13、24/x5,x); f3=-24/x5; e=5*10(-5); a=1; b=2; x1=fminbnd(f3,1,2); for n=2: Rn=-(b-a)/180*(b-a)/(2*n)4*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=NewtonCotes(
14、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 xk=fun(1,:); fk=fun(2,:); a=min(xk); b=max(xk); n=mm-1; elseif nargin=4 xk=linspace(a,b,n+1); if isa(fun,function_handle) fx=fun(xk); else error(fun积分函
15、数的句柄，且必须能够接受矢量输入！) end else error(输入参数错误，请参考函数帮助！) 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)/n*quadl(t)intfun(t,n,k),0,n); end （3）在命令窗口输入edit intfun来建立intfun
16、.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，牛顿-科特斯数公式的阶数，必须满足1n=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专心-专注-专业