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

Matlab实现数值分析插值及积分摘要：数值分析(numericalanalysis)是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。在实际生产实践中，常常将实际问题转化为数学模型来解决，这个过程就是数学建模。学习数值分析这门课程可以让我们学到很多的数学建模方法。分别运用matlab数学软件编程来解决插值问题和数值积分问题。题目中的要求是计算差值和积分，对于问题一，可以分别利用朗格朗日插值公式，牛顿插值公式，埃特金逐次线性插值公式来进行编程求解，具体matlab代码见正文。编程求解出来的结果为：=+。其中Aitken插值计算的结果图如下：对于问题二，可以分别利用复化梯形公式，复化的辛卜生公式，复化的柯特斯公式编写程序来进行求解，具体matlab代码见正文。编程求解出来的结果为：0.6932其中复化梯形公式计算的结果图如下：问题重述问题一：已知列表函数表格101234121782257分别用拉格朗日，牛顿，埃特金插值方法计算。问题二：用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式计算积分，使精度小于5。问题解决问题一：插值方法对于问题一，用三种差值方法：拉格朗日，牛顿，埃特金差值方法来解决。一、拉格朗日插值法：拉格朗日插值多项式如下：首先构造1n个插值节点nxxx,,,10上的n插值基函数，对任一点ix所对应的插值基函数)(xli，由于在所有),,1,1,,1,0(niijxj取零值，因此)(xli有因子)())(()(110niixxxxxxxx。又因)(xli是一个次数不超过n的多项式，所以只可能相差一个常数因子，固)(xli可表示成：)())(()()(110niiixxxxxxxxAxl利用1)(iixl得：)())(()(1110niiiiiixxxxxxxxA于是),,2,1,0()())(()()())(()()(110110nixxxxxxxxxxxxxxxxxlniiiiiiniii因此满足iinyxL)(),2,1,0(ni的插值多项式可表示为：njjjnxlyxL0)()(从而n次拉格朗日插值多项式为：),,2,1,0()()(0nixlyxLnjijjin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-21；而conv（）表示多项式系数乘积的结果，例如conv（poly（1），poly（1））输出为1-21；所以程序最后结果为x^n+x^n-1+……+x^2+x+1（n的值据结果的长度为准）的对应系数。在命令窗口输入editlagran来建立lagran.m文件，文件中的程序如下：function[c,l]=lagran(x,y)w=length(x);n=w-1;l=zeros(w,w);fork=1:n+1v=1;forj=1:n+1ifk~=jv=conv(v,poly(x(j)))/(x(k)-x(j));endendl(k,:)=v;endc=y*l;输入：x=[01234];y=[121782257];lagran(x,y)运行结果为ans=1.0000-0.0000-0.000001.0000结果为：=+。如图表1：图表1二．牛顿插值法newton插值多项式的表达式如下：010011()()()()()nnnNxccxxcxxxxxx其中每一项的系数ci的表达式如下：12011010[,,,][,,,][,,,]iiiiifxxxfxxxcfxxxxx即为f(x)在点01,,,ixxx处的i阶差商，（()iifxfx，1,2,,in），由差商01[,,,]ifxxx的性质可知：01001[,,,]()iiijjkjkkjfxxxfxxx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-21；而conv（）表示多项式系数乘积的结果，例如conv（poly（1），poly（1））输出为1-21；所以程序最后结果为x^n+x^n-1+……+x^2+x+1（n的值据结果的长度为准）的对应系数。在命令窗口输入editnowpoly来建立newpoly.m文件，文件中的程序如下：function[c,d]=newpoly(x,y)n=length(x);d=zeros(n,n);d(:,1)=y';forj=2:nfork=j:nd(k,j)=(d(k,j-1)-d(k-1,j-1))/(x(k)-x(k-j+1));endendc=d(n,n);fork=(n-1):-1:1c=conv(c,poly(x(k)));m=length(c);c(m)=c(m)+d(k,k);end输入：x=[01234];y=[121782257];newpoly(x,y)运行结果为ans=10001所以=+。如图表2：图表2三．埃特金插值法：Aitken插值公式如下：递推表达式为：=+当n=1时，=+当n=2时，=+其中的带入递推表达式求得。由此递推下去，最终得到的结果。matlab编程：编程思想：埃特金插值多项式又称作Aitken逐次线性插值多项式，根据公式的特点，可以利用2次嵌套循环将公式表示出来。在命令窗口输入editAitken来建立Aitken.m文件，文件中的程序如下：functionf=Aitken(x,y)symsz;n=length(x);y1(1:n)=z;fori=1:n-1forj=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=[01234];y=[121782257];Aitken(x,y)运行结果为ans=z^4+1所以=+。如图表3：图表3问题二：复化积分对于问题二来说，用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式结局问题（计算积分，使精度小于5）。一复化的梯形公式：复化梯形的迭代公式为：;matlab编程：程序1（求f（x）的n阶导数：）在命令窗口输入editqiudao来建立qiudao.m文件，文件中的程序如下：functiond=qiudao(x,n)symsx;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：在命令窗口输入edittixing来建立tixing.m文件，文件中的程序如下：functiony=tixing()symsx；%定义自变量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值forn=2:1000000；%求等分数nRn=-(b-a)/12*((b-a)/n)^2*f2(x1)；%计算余项ifabs(Rn)e%用余项进行判断break%符合要求时结束endendh=(b-a)/n；%求hTn1=0；fork=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阶导数：）在命令窗口输入editqiudao来建立qiudao.m文件，文件中的程序如下：functiond=qiudao(x,n)symsx;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：在命令窗口输入editxinpusheng来建立xinpusheng.m文件，文件中的程序如下：functiony=xinpusheng()symsx;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);forn=2:1000000Rn=-(b-a)/180*((b-a)/(2*n))^4*f2(x1);ifabs(Rn)ebreakendendh=(b-a)/n;Sn1=0;Sn2=0;fork=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函数：输入：clearxinpusheng()运行结果为用Simpson公式计算的结果Sn=0.6932等分数n=4结果如图表7图表7三复化的柯特斯公式：牛顿-柯特斯公式如下：matlab编程：(1)在命令窗口输入editNewtonCotes来建立NewtonCotes.m文件，文件中的程序如下：function[y,Ck,Ak]=NewtonCotes(fun,a,b,n)ifnargin==1[mm,nn]=size(fun);ifmm=8error('为了保证NewtonCotes积分的稳定性，最多只能有9个等距节点！')elseifnn~=2error('fun构成应为：第一列为x，第二列为y，并且个数为小于10的等距节点！')endxk=fun(1,:);fk=fun(2,:);a=min(xk);b=max(xk);0()()nbkkakfxdxAfx011011()()()()()()()()()bbkknkkaakkkkkknxxxxxxxxAlxdxdxxxxxxxxxn=mm-1;elseifnargin==4xk=linspace(a,b,n+1);ifisa(fun,'function_handle')fx=fun(xk);elseerror('fun积分函数的句柄，且必须能够接受矢量输入！')endelseerror('输入参数错误，请参考函数帮助！