复化梯形公式和复化Simpson公式文档格式.docx
《复化梯形公式和复化Simpson公式文档格式.docx》由会员分享,可在线阅读,更多相关《复化梯形公式和复化Simpson公式文档格式.docx(9页珍藏版)》请在冰豆网上搜索。
2*e*p(*)+**e*p(*)
程序2:
clc
clear
syms*%定义自变量*
f=inline('
**e*p(*)'
'
*'
)%定义函数f(*)=**e*p(*),换函数时只需换该函数表达式即可
f2=inline('
(2*e*p(*)+**e*p(*))'
)%定义f(*)的二阶导数,输入程序1里求出的f2即可。
f3='
-(2*e*p(*)+**e*p(*))'
%因fminbnd〔〕函数求的是表达式的最小值,且要求表达式带引号,故取负号,以便求最大值
e=5*10^(-8)%精度要求值
a=1%积分下限
b=2%积分上限
*1=fminbnd(f3,1,2)%求负的二阶导数的最小值点,也就是求二阶导数的最大值点对应的*值
forn=2:
1000000%求等分数n
Rn=-(b-a)/12*((b-a)/n)^2*f2(*1)%计算余项
ifabs(Rn)<
e%用余项进展判断
break%符合要求时完毕
end
h=(b-a)/n%求h
Tn1=0
fork=1:
n-1%求连加和
*k=a+k*h
Tn1=Tn1+f(*k)
Tn=h/2*((f(a)+2*Tn1+f(b)))
z=e*p
(2)
R=Tn-z%求值与计算值的差
fprintf('
用复化梯形算法计算的结果Tn='
)
disp(Tn)
等分数n='
disp(n)%输出等分数
值与计算值的误差R='
disp(R)
输出结果显示:
用复化梯形算法计算的结果Tn=7.3891
等分数n=7019
值与计算值的误差R=2.8300e-008
2.Simpson公式
程序1:
〔求f〔*〕的n阶导数〕:
输入n=4
4*e*p(*)+**e*p(*)
(4*e*p(*)+**e*p(*))'
)%定义f(*)的四阶导数,输入程序1里求出的f2即可
-(4*e*p(*)+**e*p(*))'
%因fminbnd〔〕函数求的是表达式的最小值,且要求表达式带引号,故取负号,一边求最大值
*1=fminbnd(f3,1,2)%求负的四阶导数的最小值点,也就是求四阶导数的最大值点对应的*值
Rn=-(b-a)/180*((b-a)/(2*n))^4*f2(*1)%计算余项
Sn1=0
Sn2=0
fork=0:
n-1%求两组连加和
*k1=*k+h/2
Sn1=Sn1+f(*k1)
Sn2=Sn2+f(*k)
Sn=h/6*(f(a)+4*Sn1+2*(Sn2-f(a))+f(b))%因Sn2多加了k=0时的值,故减去f(a)
R=Sn-z%求值与计算值的差
用Simpson公式计算的结果Sn='
disp(Sn)
disp(n)
用Simpson公式计算的结果Sn=7.3891
等分数n=24
值与计算值的误差R=2.7284e-008
用复化梯形公式计算的结果为:
7.3891,与准确解的误差为:
2.8300e-008。
等分数n=7019
用复化Simpson公式计算的结果为:
2.7284e-008。
等分数n=24
3、柯斯特公式求积分:
程序代码:
〔1〕function[y,Ck,Ak]=NewtonCotes(fun,a,b,n)
ifnargin==1
[mm,nn]=size(fun);
ifmm>
=8
error('
为了保证NewtonCotes积分的稳定性,最多只能有9个等距节点!
elseifnn~=2
fun构成应为:
第一列为*,第二列为y,并且个数为小于10的等距节点!
end
*k=fun(1,:
);
fk=fun(2,:
a=min(*k);
b=ma*(*k);
n=mm-1;
elseifnargin==4
*k=linspace(a,b,n+1);
ifisa(fun,'
function_handle'
f*=fun(*k);
else
fun积分函数的句柄,且必须能够承受矢量输入!
else
输入参数错误,请参考函数帮助!
end
Ck=cotescoeff(n);
Ak=(b-a)*Ck;
y=Ak*f*'
;
〔2〕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);
〔3〕functionf=intfun(t,n,k)
f=1;
fori=[0:
k-1,k+1:
n]
f=f.*(t-i);
代码解释:
function[y,Ck,Ak]=NewtonCotes(fun,a,b,n)
%y=NewtonCotes(fun,a,b,n)
%牛顿-科特斯数值积分公式
%参数说明:
%fun,积分表达式,这里有两种选择
%
(1)积分函数句柄,必须能够承受矢量输入,比方fun=(*)sin(*).*cos(*)
%
(2)*,y坐标的离散点,第一列为*,第二列为y,必须等距,且节点的个数小于9,比方:
fun=[1:
8;
sin(1:
8)]'
%如果fun的表采用第二种方式,则只需要输入第一个参数即可,否则还要输入a,b,n三个参数
%a,积分下限
%b,积分上限
%n,牛顿-科特斯数公式的阶数,必须满足1<
n<
7,因为n>
=8时不能保证公式的稳定性
%
(1)n=1,即梯形公式
%
(2)n=2,即辛普森公式
%(3)n=4,即科特斯公式
%y,数值积分结果
%Ck,科特斯系数
%Ak,求积系数
%
%E*ample
%fun1=(*)sin(*);
%必须可以承受矢量输入
%fun2=[0:
0.1:
0.5;
sin(0:
0.5)];
%最多8个点,必须等距
%y1=NewtonCotes(fun1,0,0.5,6)
%y2==NewtonCotes(fun2)
%计算积分节点*k和节点函数值f*
*k=linspace(a,b,n+1);
%计算科特斯系数
%计算求积系数
%求和算积分
functionCk=cotescoeff(n)
%由于科特斯系数最多7阶,为了方便我们可以直接使用,省得每次都计算
%A1=[1,1]/2
%A2=[1,4,1]/6
%A3=[1,3,3,1]/8
%A4=[7,32,12,32,1]/90
%A5=[19,75,50,50,75,19]/288
%A6=[41,216,27,272,27,216,41]/840
%A7=[751,3577,1323,2989,2989,1323,3577,751]/17280
%当时为了表达公式,我们使用程序计算n阶科特斯系数
Ck(i)=(-1)^(n-k)/factorial(k)/factorial(n-k)/n*quadl((t)intfun(t,n,k),0,n);
functionf=intfun(t,n,k)
%科特斯系数中的积分表达式
输出结果:
fun=(*)e*p(*);
a=-1;
b=1;
n=4;
NewtonCotes(fun,a,b,n)
ans=2.3505
二、三点数值微分