数值积分与数值微分.docx
《数值积分与数值微分.docx》由会员分享，可在线阅读，更多相关《数值积分与数值微分.docx（20页珍藏版）》请在冰点文库上搜索。
数值积分与数值微分
8数值积分与数值微分
8.1例题解答
例8.1给定积分
分别用梯形公式、
公式、
公式作近似计算.
解:
先输入主要初始参数
>>a=0.5;
>>b=1;
>>f=inline（'x^（1/2）'）;
%梯形公式
>>I1=（b-a）/2*（feval（f,a）+feval（f,b））
I1=
0.426776695296637
%simpson公式
>>I2=（b-a）/6*（feval（f,a）+4*feval（f,（a+b）/2）+feval（f,b））
I2=
0.430934033027025
%Cotes公式（n=4）
>>tc=0;
>>C0=[73212327];
>>fori=0:
4
tc=tc+C0（i+1）*feval（f,a+i*（b-a）/4）;
end
>>I3=（b-a）/90*tc
I3=
0.430964070495876
%准确值
>>I=int（char（f）,a,b）
>>vpa（I）
I=
-1/6*2^（1/2）+2/3
ans=
0.43096440627115082519971854596505
例8.2对积分
为使其精度达到
.
若用复化梯形公式,应将[0,1]多少等分?
若用复化
公式,应将[0,1]多少等分?
解:
直接按余项计算即可.
复化梯形公式的余项为:
复化
公式余项为:
对于
在课本中我们已证得以下不等式成立:
直接利用上述不等式关系解答本题.
先输入误差精度:
>>Eps=1E-4
Eps=
1.000000000000000e-004
（1）复化梯形公式
>>h1=sqrt（Eps/abs（-（1-0）/12*1/（2+1）））%先求出步长
h1=
0.060000000000000
>>N1=ceil（1/h1）%向上取整,得到等分区间数
N1=
17
故可将区间17等分即可达到所要求的精度.
（2）复化
公式
>>h2=power（Eps/abs（-（1-0）/180*1/（1+4））,1/4）%先求出步长
h2=
0.547722557505166
>>N2=ceil（1/h2）%向上取整,得到等分区间数
N2=
2
故可将区间2等分即可达到所要求的精度.
Ø扩展:
1）Matlab中复化梯形公式命令为I=trapz（x,y）,复化
公式命令为quad（）.
2）Matlab中有四个取整函数,分别为ceil（）,floor（）,fix（）,round（）,分别表示向正无穷大方向取整、向负无穷大方向取整、向靠近零方向舍入和四舍五入.
例8.3对积分
利用变步长方法求其近似值,使其精度达到
.
解：
利用变步长法前先建立三种变步长复化积分公式的函数.注意在Matlab中直接用sin（0）/0得不到1，
，因此解此题时我们改用求极限的方法得到函数值,此函数名为limit（）.
先建立三种复化公式的函数文件,它们分别为复化梯形公式trap.m、复化
公式为simpson.m、
公式为cotes.m,三个函数的源程序如下:
（1）复化梯形公式trap.m
functionT=trap（f,a,b,n）
%trap.m
%复化梯形公式求积分值
%f为积分函数
%[a,b]为积分区间
%n是等分区间份数
h=（b-a）/n;%步长
T=0;
fork=1:
（n-1）
x0=a+h*k;
T=T+limit（f,x0）;
end
T=h*（limit（f,a）+limit（f,b））/2+h*T;
T=double（T）;
（2）复化
公式simpson.m:
functionS=simpson（f,a,b,n）
%simpson.m
%Simpson公式求积分值
%f为积分函数
%[a,b]为积分区间
%n是等分区间份数
h=（b-a）/（2*n）;%步长
s1=0;
s2=0;
fork=1:
n
x0=a+h*（2*k-1）;
s1=s1+limit（f,x0）;
end
fork=1:
（n-1）
x0=a+h*2*k;
s2=s2+limit（f,x0）;
end
S=h*（limit（f,a）+limit（f,b）+4*s1+2*s2）/3;
S=double（S）;
（3）复化
公式cotes.m:
functionC=cote（f,a,b,n）
%cote.m
%复化cotes公式求积分值
%f为积分函数
%[a,b]为积分区间
%n是等分区间份数
h=（b-a）/n;%步长
C=0;
fori=1:
（n-1）
x0=a+i*h;
C=C+14*limit（f,x0）;
end
fork=0:
（n-1）
x0=a+h*k;
s=32*limit（f,x0+h*1/4）+12*limit（f,x0+h*1/2）+32*limit（f,x0+h*3/4）;
C=C+s;
end
C=C+7*（limit（f,a）+limit（f,b））;
C=C*h/90;
C=double（C）;
再编写主程序调用这三个函数,主程序名为ex8_3.m,源程序如下:
%ex8_3.m
clc;
symsx;
f=sym（'sin（x）/x'）;
a=0;b=1;%积分上下限
n=20;%作1,2,3,…,20次区间等分
%复化梯形公式
T=zeros（n,1）;
fori=1:
n
T（i）=trap（f,a,b,i）;
end
%复化Simpson公式;
S=zeros（n,1）;
fori=1:
n
S（i）=simpson（f,a,b,i）;
end
%复化Cotes公式
C=zeros（n,1）;
fori=1:
n
C（i）=cote（f,a,b,i）;
end
%准确值
I=int（f,a,b）;
I=double（I）;
%画图作出直观观察
x=[];
x=1:
n;
figure;
plot（x,ones（1,n）*I,'-'）;
holdon;
plot（x,T','r--','LineWidth',2）;
plot（x,S','m.-','LineWidth',1）;
plot（x,C','c:
','LineWidth',1.5）;
gridon;
title（'三种复化公式积分效果对比图'）;
legend（'准确值曲线','复化梯形公式','复化Simpson公式','复化Cotes公式'）;
holdoff;
disp（'复化梯形公式复化Simpson公式复化Cotes公式'）;
disp（[TSC]）;
在Matlab命令窗口中输入ex8_3,得到如下积分结果:
复化梯形公式复化Simpson公式复化Cotes公式
0.9207354924039480.9461458822735870.946083004063674
0.9397932848061770.9460869339517940.946083069350917
0.9432914291323370.9460838313116990.946083070278278
0.9445135216653900.9460833108884720.946083070351379
0.9450787809534020.9460831688380730.946083070363043
0.9453857307668590.9460831178428670.946083070365797
0.9455707762562460.9460830959894030.946083070366633
0.9456908635827010.9460830853849480.946083070366936
0.9457731885497520.9460830797420530.946083070367061
0.9458320718669050.9460830765177320.946083070367118
0.9458756371191980.9460830745679370.946083070367147
0.9459087710738650.9460830733331140.946083070367161
0.9459345564884750.9460830725204760.946083070367170
0.9459550160561140.9460830719680570.946083070367174
0.9459715215529850.9460830715819640.946083070367177
0.9459850299343860.9460830713055620.946083070367179
0.9459962252423760.9460830711034890.946083070367180
0.9460056069439780.9460830709529980.946083070367181
0.9460135466239660.9460830708390650.946083070367182
0.9460203253550250.9460830707515320.946083070367182
由以上结果可看到复化梯形公式有一个上升接近准确值过程,而复化
公式积分结果和复化
公式积分的结果基本上和准确值的曲线重叠在一块,可见它们的精度是相当高的.
例8.4用
积分法计算
精度
.
解：
编写
积分法的函数M文件,源程序如下（romberg.m）:
function[I,T]=romberg（f,a,b,n,Eps）
%Romberg积分计算
%f为积分函数
%[a,b]为积分区间
%n+1是T数表的列数目
%Eps为迭代精度
%返回值中I为积分结果,T是积分表
ifnargin<5
Eps=1E-6;
end
m=1;
h=（b-a）;
err=1;
j=0;
T=zeros（4,4）;
T（1,1）=h*（limit（f,a）+limit（f,b））/2;
while（（err>Eps）&（jj=j+1;
h=h/2;
s=0;
forp=1:
m
x0=a+h*（2*p-1）;
s=s+limit（f,x0）;
end
T（j+1,1）=T（j,1）/2+h*s;
m=2*m;
fork=1:
j
T（j+1,k+1）=T（j+1,k）+（T（j+1,k）-T（j,k））/（4^k-1）;
end
err=abs（T（j,j）-T（j+1,k+1））;
end
I=T（j+1,j+1）;
ifnargout==1
T=[];
end
将上述源程序另存为romberg.m后,进入计算:
>>symsx;%创建符号变量
>>f=sym（'sin（x）/x'）%符号函数
f=
sin（x）/x
>>[I,T]=romberg（f,0,1,3,1E-6）%积分计算
I=
0.9461
T=
0.92070000
0.93980.9461000
0.94450.94610.946100
0.94570.94610.94610.94610
0.94600.94610.94610.94610.9461
其中T为
积分表,由输出结果可知计算结果为
.
例8.5利用
积分法求
.
解:
直接利用例8.4的函数即可.
>>symsx;
>>f=sym（'4/（1+x^2）'）
f=
4/（1+x^2）
>>[I,T]=romberg（f,0,1,5,1E-6）
I=
3.1416
T=
3.000000000
3.10003.13330000
3.13123.14163.1421000
3.13903.14163.14163.141600
3.14093.14163.14163.14163.14160
3.14143.14163.14163.14163.14163.1416
故结果为
.
例8.6对积分
构造其
型求积公式.
解:
取
得到方程组:
这是一个抽象方程组,可以利用Matlab的符号法来解之,该函数名为solve（）:
%直接输入解之:
>>x=solve（'A0+A1=2','A0*x0+A1*x1=0','A0*x0^2+A1*x1^2=2/3','A0*x0^3+A1*x1^3=0','A0,A1,x0,x1'）
x=
A0:
[2x1sym]
A1:
[2x1sym]
x0:
[2x1sym]
x1:
[2x1sym]
%显示结果
>>x.A0,x.A1,x.x0,x.x1
ans=
1
1
ans=
1
1
ans=
1/3*3^（1/2）
-1/3*3^（1/2）
ans=
-1/3*3^（1/2）
1/3*3^（1/2）
因此得到两组解为:
或:
求积公式为
例8.7分别利用
公式及
公式计算积分:
解:
（1）先计算准确值:
>>I=int（'（x+1.5）^（1/2）',-1,1）
I=
2.3995291230781336018654631663256
（2）两点
公式:
由课本表中查得两点
公式的求积系数为1,求积节点为
>>f=inline（'（x+1.5）^（1/2）'）
f=
Inlinefunction:
f（x）=（x+1.5）^（1/2）
>>I1=1*feval（f,-0.57735）+1*feval（f,0.57735）
I1=
2.401848214499055
（3）两个节点梯形公式:
>>f=inline（'（x+1.5）^（1/2）'）
>>a=-1;%求积区间
>>b=1;
>>I2=（b-a）/2*（feval（f,a）+feval（f,b））%两点梯形公式
f=
Inlinefunction:
f（x）=（x+1.5）^（1/2）
I2=
2.288245611270737
（4）三点
公式:
查表有
:
>>f=inline（'（x+1.5）^（1/2）'）
>>x=[0.7745966692,-0.7745966692,0]%节点系数
>>A=[0.555555556,0.555555556,0.888888889]%求积系数
>>I3=0;
>>fori=1:
length（x）
I3=I3+feval（f,x（i））*A（i）;%三点Gauss-Legendre公式
end
>>I3%显示结果
f=
Inlinefunction:
f（x）=（x+1.5）^（1/2）
x=
0.774596669200000-0.7745966692000000
A=
0.5555555560000000.5555555560000000.888888889000000
I3=
2.399708072133707
（5）三个节点
公式:
>>f=inline（'（x+1.5）^（1/2）'）
>>a=-1;%积分区间
>>b=1;
>>x=[a,（a+b）/2,b]%积分节点
>>A=[141]%积分节点系数
>>I4=0;
>>fori=1:
length（x）
I4=I4+feval（f,x（i））*A（i）;%三点Simpson公式
end
>>I4=I4*（b-a）/6%显示计算结果
f=
Inlinefunction:
f（x）=（x+1.5）^（1/2）
x=
-101
A=
141
I4=
2.395741698945698
例8.8利用
求积公式求积分
.（准确值为1.）
解:
先作变换,将积分区间变换到
上,令
则有:
于是就可以利用
求积公式来解题.
（1）两点
求积公式:
>>f=inline（'sin（pi*（t+1）/4）'）
f=
Inlinefunction:
f（t）=sin（pi*（t+1）/4）
>>I1=1*feval（f,-0.57735）+1*feval（f,0.57735）;
>>I1=I1*pi/4
I1=
0.998472716275797
（2）三点
求积公式:
>>f=inline（'sin（pi*（t+1）/4）'）
f=
Inlinefunction:
f（t）=sin（pi*（t+1）/4）
>>x=[0.7745966692,-0.7745966692,0]%节点系数
x=
0.774596669200000-0.7745966692000000
>>A=[0.555555556,0.555555556,0.888888889]%求积系数
A=
0.5555555560000000.5555555560000000.888888889000000
>>I2=0;
>>fori=1:
length（x）
I2=I2+feval（f,x（i））*pi/4*A（i）;
end
>>I2
I2=
1.000008122033779
例8.9计算积分
精度要求
.
解:
先利用余项来估计误差精度,再进行计算即可.
这里利用
求积公式余项作误差估计.
余项为:
对右边进行放大得到:
%先搜索满足上述不等式的项数n
>>i=1;
>>E=pi*exp
（1）/（power（2,2*i+1）*factorial（2*i+2））;
>>whileE>1E-6
i=i+1;
E=pi*exp
（1）/（power（2,2*i+1）*factorial（2*i+2））;
end
>>i,E%显示结果
i=
4
E=
4.596331680902589e-009
可见只要计算前4项即可.
>>f=inline（'exp（x）'）
f=
Inlinefunction:
f（x）=exp（x）
>>n=i;
>>I=0;
>>forj=0:
n
x=cos（pi*（2*j+1）/（2*n+2））;
I=I+feval（f,x）;
end
>>I=I*pi/（n+1）
I=
3.977463258776695
升级会员 