matlab实现数值分析报告插值及积分.docx
《matlab实现数值分析报告插值及积分.docx》由会员分享，可在线阅读，更多相关《matlab实现数值分析报告插值及积分.docx（18页珍藏版）》请在冰点文库上搜索。
matlab实现数值分析报告插值及积分
Matlab实现数值分析插值及积分
摘要：
数值分析（numericalanalysis）是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。
在实际生产实践中，常常将实际问题转化为数学模型来解决，这个过程就是数学建模。
学习数值分析这门课程可以让我们学到很多的数学建模方法。
分别运用matlab数学软件编程来解决插值问题和数值积分问题。
题目中的要求是计算差值和积分，对于问题一，可以分别利用朗格朗日插值公式，牛顿插值公式，埃特金逐次线性插值公式来进行编程求解，具体matlab代码见正文。
编程求解出来的结果为：
=+。
其中Aitken插值计算的结果图如下：
对于问题二，可以分别利用复化梯形公式，复化的辛卜生公式，复化的柯特斯公式编写程序来进行求解，具体matlab代码见正文。
编程求解出来的结果为：
0.6932
其中复化梯形公式计算的结果图如下：
问题重述
问题一：
已知列表函数
表格1
0
1
2
3
4
1
2
17
82
257
分别用拉格朗日，牛顿，埃特金插值方法计算。
问题二：
用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式计算积分，使精度小于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-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+1
v=1;
forj=1:
n+1
ifk~=j
v=conv（v,poly（x（j）））/（x（k）-x（j））;
end
end
l（k,:
）=v;
end
c=y*l;
输入：
>>x=[01234];
>>y=[121782257];
>>lagran（x,y）
运行结果为
ans=
1.0000-0.0000-0.000001.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）=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:
n
fork=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）;
fork=（n-1）:
-1:
1
c=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-1
forj=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=[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=2
ans=
2/x^3
结果为：
f2=2/x^3
如图表4：
图表4
程序2：
在命令窗口输入edittixing来建立tixing.m文件，文件中的程序如下：
functiony=tixing（）
symsx；%定义自变量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值
forn=2:
1000000；%求等分数n
Rn=-（b-a）/12*（（b-a）/n）^2*f2（x1）；%计算余项
ifabs（Rn）break%符合要求时结束
end
end
h=（b-a）/n；%求h
Tn1=0；
fork=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迭代公式为：
;
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=4
ans=
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:
1000000
Rn=-（b-a）/180*（（b-a）/（2*n））^4*f2（x1）;
ifabs（Rn）break
end
end
h=（b-a）/n;
Sn1=0;
Sn2=0;
fork=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）在命令窗口输入editNewtonCotes来建立NewtonCotes.m文件，文件中的程序如下：
function[y,Ck,Ak]=NewtonCotes（fun,a,b,n）
ifnargin==1
[mm,nn]=size（fun）;
ifmm>=8
error（'为了保证NewtonCotes积分的稳定性，最多只能有9个等距节点！
'）
elseifnn~=2
error（'fun构成应为：
第一列为x，第二列为y，并且个数为小于10的等距节点！
'）
end
xk=fun（1,:
）;
fk=fun（2,:
）;
a=min（xk）;
b=max（xk）;
n=mm-1;
elseifnargin==4
xk=linspace（a,b,n+1）;
ifisa（fun,'function_handle'）
fx=fun（xk）;
else
error（'fun积分函数的句柄，且必须能够接受矢量输入！
'）
end
else
error（'输入参数错误，请参考函数帮助！
'）
end
Ck=cotescoeff（n）;
Ak=（b-a）*Ck;
y=Ak*fx';
（2）在命令窗口输入editcotescoeff来建立cotescoeff.m文件，文件中的程序如下：
functionCk=cotescoeff（n）
fori=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）在命令窗口输入editintfun来建立intfun.m文件，文件中的程序如下：
functionf=intfun（t,n,k）
f=1;
fori=[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
升级会员 