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

数值分析

学院：计算机
专业：计算机科学与技术
班级：xxx
学号：xxx
姓名：xxx
指导教师：xxx

数值分析
摘要：
数值分析(numerical analysis)是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。在实际生产实践中，常常将实际问题转化为数学模型来解决，这个过程就是数学建模。学习数值分析这门课程可以让我们学到很多的数学建模方法。
分别运用matlab数学软件编程来解决插值问题和数值积分问题。题目中的要求是计算差值和积分，对于问题一，可以分别利用朗格朗日插值公式，牛顿插值公式，埃特金逐次线性插值公式来进行编程求解，具体matlab代码见正文。编程求解出来的结果为：P4x=x4+1。
其中Aitken插值计算的结果图如下：

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

问题重述
问题一：已知列表函数

表格 1
x
0
1
2
3
4
fx
1
2
17
82
257

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

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

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

首先构造个插值节点上的插值基函数，对任一点所对应的插值基函数，由于在所有取零值，因此有因子。又因是一个次数不超过的多项式，所以只可能相差一个常数因子，固可表示成：

利用得：

于是

因此满足 的插值多项式可表示为：

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

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
结果为：P4x=x4+1。
如图表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）=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插值公式如下：
递推表达式为：
P0，1，…，nx = x-xnxn-1-xnP0，1，…，n-1(x) + x-xn-1xn-xn-1P0，1，…，n-2，n(x)

当n=1时，
P0，1x = x-x1x0-x1fx0 + x-x0x1-x0fx1
当n=2时，
P0，1，2x = x-x2x1-x2P0，1(x) + x-x1x2-x1P0，2(x)
其中的P0，2(x)带入递推表达式求得。
由此递推下去，最终得到P0，1，…，nx的结果。
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
问题二：复化积分
对于问题二来说，用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式结局问题（计算积分121xdx，使精度小于5×10-5）。
一 复化的梯形公式：

复化梯形的迭代公式为：
abf(x)dx=h/2fa+2j=1n-1fxj+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)<e %用余项进行判断
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迭代公式为：
abf(x)dx=h/3fa+2j=1(n2)-1fx2j+4j=1(n2)fx2j-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)<e
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=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<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