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

数值分析

学院：计算机

专业：计算机科学与技术

班级：xxx

学号：xxx

姓名：xxx

指导教师：xxx

数值分析

摘要：

数值分析(numerical analysis)是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。在实际生产实践中，常常将实际问题转化为数学模型来解决，这个过程就是数学建模。学习数值分析这门课程可以让我们学到很多的数学建模方法。

分别运用matlab数学软件编程来解决插值问题和数值积分问题。题目中的要求是计算差值和积分，对于问题一，可以分别利用朗格朗日插值公式，牛顿插值公式，埃特金逐次线性插值公式来进行编程求解，具体matlab代码见正文。编程求解出来的结果为：P4x=x4+1。其中Aitken插值计算的结果图如下：

对于问题二，可以分别利用复化梯形公式，复化的辛卜生公式，复化的柯特斯公式编写程序来进行求解，具体matlab代码见正文。编程求解出来的结果为：0.6932

其中复化梯形公式计算的结果图如下：

问题重述

问题一：已知列表函数

表格 1

分别用拉格朗日，牛顿，埃特金插值方法计算P 4 x 。

问题二：用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式计算积分 1

x 21dx ，使精度小于5×10?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 得：

)

())(()(1

110n 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 =的插值多项式可表示为：

∑==n

j j j n x l y x L 0

)()(

从而n 次拉格朗日插值多项式为：

),,2,1,0()

()(0

n i x l y x L n

j 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~=j

v=conv(v,poly(x(j)))/(x(k)-x(j)); end end l(k,:)=v; end c=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 结果为：P 4 x =x 4+1。 如图表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 ???的性质可知：

()

010

1

[,,,]()i

i

i j j k j k k j

f 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:n

for k=j:n

d(k,j)=(d(k,j-1)-d(k-1,j-1))/(x(k)-x(k-j+1));

end

end

c=d(n,n);

for k=(n-1):-1:1

c=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

所以P4x=x4+1。

如图表2：

图表2

三．埃特金插值法：

Aitken插值公式如下：

递推表达式为：

P

0，1，…，n x= x?x n

x n?1?x n

P

0，1，…，n?1

(x)+ x?x n?1

x n?x n?1

P

0，1，…，n?2，n

(x)

当n=1时，

P

0，1x= x?x1

x0?x1

f x0+ x?x0

x1?x0

f x1

当n=2时，

P

0，1，2x= x?x2

x1?x2

P

0，1

(x)+ x?x1

x2?x1

P

0，2

(x)

其中的P

0，2

(x)带入递推表达式求得。

由此递推下去，最终得到P

0，1，…，n

x的结果。

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-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 =

z^4 + 1

所以P4x=x4+1。

如图表3：

图表 3

问题二：复化积分

对于问题二来说，用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式结局问题（计算积分 1

x 21dx ，使精度小于5×10?5）。

一 复化的梯形公式：

复化梯形的迭代公式为：

f(x)dx b

a =h/2 f a +2 f x j n ?1

j=1+f(b) ; 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 = 2 ans = 2/x^3 结果为：

f2 =2/x^3

如图表4：

图表4

程序2：

在命令窗口输入edit tixing 来建立tixing.m文件，文件中的程序如下：

function y=tixing()

syms x ；%定义自变量x

f=inline('1/x','x')；%定义函数f(x)= 1/x

f2=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；%求等分数n

Rn=-(b-a)/12*((b-a)/n)^2*f2(x1)；%计算余项

if abs(Rn)

break % 符合要求时结束

end

end

h=(b-a)/n；%求h

Tn1=0；

for k=1:n-1 %求连加和

xk=a+k*h；

Tn1=Tn1+f(xk)；

end

Tn=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 迭代公式为：

f(x)dx

b

a =h/3 f a +2 f x 2j (n 2

)?1j=1+4 f x 2j ?1 (n 2

)

j=1+f(b) ;

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 = 4 ans = 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:1000000

Rn=-(b-a)/180*((b-a)/(2*n))^4*f2(x1);

if abs(Rn)

break

end

end

h=(b-a)/n;

Sn1=0;

Sn2=0;

for k=0:n-1

xk=a+k*h;

xk1=xk+h/2;

Sn1=Sn1+f(xk1);

Sn2=Sn2+f(xk);

end

Sn=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>=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 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 -+-+----==----??

n=mm-1;

elseif nargin==4

xk=linspace(a,b,n+1);

if isa(fun,'function_handle')

fx=fun(xk);

else

error('fun积分函数的句柄，且必须能够接受矢量输入！')

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.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