复化梯形公式和复化Simpson公式.docx

文档介绍：该【复化梯形公式和复化Simpson公式 】是由【HJ01】上传分享，文档一共【12】页，该文档可以免费在线阅读，需要了解更多关于【复化梯形公式和复化Simpson公式 】的内容，可以使用淘豆网的站内搜索功能，选择自己适合的文档，以下文字是截取该文章内的部分文字，如需要获得完整电子版，请下载此文档到您的设备，方便您编辑和打印。Documentnumber:NOCG-YUNOO-BUYTT-UU986-1986UT
复化梯形公式和复化Simpson公式
数值计算方法上机题目3
一、计算定积分的近似值:
要求:
若用复化梯形公式和复化Simpson公式计算,要求误差限,分别利用他们的余项估计对每种算法做出步长的事前估计;
分别利用复化梯形公式和复化Simpson公式计算定积分;
将计算结果与精确解比较,并比较两种算法的计算量。

程序:
程序1(求f(x)的n阶导数:
symsx
f=x*exp(x)%定义函数f(x)
n=input('输入所求导数阶数:')
f2=diff(f,x,n)%求f(x)的n阶导数
结果1
输入n=2
f2=
2*exp(x)+x*exp(x)
程序2:
clc
clear
symsx%定义自变量x
f=inline('x*exp(x)','x')%定义函数f(x)=x*exp(x),换函数时只需换该函数表达式即可
f2=inline('(2*exp(x)+x*exp(x))','x')%定义f(x)的二阶导数,输入程序1里求出的f2即可。
f3='-(2*exp(x)+x*exp(x))'%因fminbnd()函数求的是表达式的最小值,且要求表达式带引号,故取负号,以便求最大值
e=5*10^(-8)%精度要求值
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)<e%用余项进行判断
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)))
z=exp(2)
R=Tn-z%求已知值与计算值的差
fprintf('用复化梯形算法计算的结果Tn=')
disp(Tn)
fprintf('等分数n=')
disp(n)%输出等分数
fprintf('已知值与计算值的误差R=')
disp(R)
输出结果显示:
用复化梯形算法计算的结果Tn=
等分数n=7019
已知值与计算值的误差R=

程序:
程序1:(求f(x)的n阶导数):
symsx
f=x*exp(x)%定义函数f(x)
n=input('输入所求导数阶数:')
f2=diff(f,x,n)%求f(x)的n阶导数
结果1
输入n=4
f2=
4*exp(x)+x*exp(x)
程序2:
clc
clear
symsx%定义自变量x
f=inline('x*exp(x)','x')%定义函数f(x)=x*exp(x),换函数时只需换该函数表达式即可
f2=inline('(4*exp(x)+x*exp(x))','x')%定义f(x)的四阶导数,输入程序1里求出的f2即可
f3='-(4*exp(x)+x*exp(x))'%因fminbnd()函数求的是表达式的最小值,且要求表达式带引号,故取负号,一边求最大值
e=5*10^(-8)%精度要求值
a=1%积分下限
b=2%积分上限
x1=fminbnd(f3,1,2)%求负的四阶导数的最小值点,也就是求四阶导数的最大值点对应的x值
forn=2:1000000%求等分数n
Rn=-(b-a)/180*((b-a)/(2*n))^4*f2(x1)%计算余项
ifabs(Rn)<e%用余项进行判断
break%符合要求时结束
end
end
h=(b-a)/n%求h
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))%因Sn2多加了k=0时的值,故减去f(a)
z=exp(2)
R=Sn-z%求已知值与计算值的差
fprintf('用Simpson公式计算的结果Sn=')
disp(Sn)
fprintf('等分数n=')
disp(n)
fprintf('已知值与计算值的误差R=')
disp(R)
输出结果显示:
用Simpson公式计算的结果Sn=
等分数n=24
已知值与计算值的误差R=
用复化梯形公式计算的结果为:,与精确解的误差为:。等分数n=7019
用复化Simpson公式计算的结果为:,与精确解的误差为:。等分数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
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)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)functionf=intfun(t,n,k)
f=1;
fori=[0:k-1,k+1:n]
f=f.*(t-i);
end
代码解释:
function[y,Ck,Ak]=NewtonCotes(fun,a,b,n)
%y=NewtonCotes(fun,a,b,n)
%牛顿-科特斯数值积分公式
%参数说明:
%fun,积分表达式,这里有两种选择
%(1)积分函数句柄,必须能够接受矢量输入,比如fun=@(x)sin(x).*cos(x)
%(2)x,y坐标的离散点,第一列为x,第二列为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,求积系数
%
%Example
%fun1=@(x)sin(x);%必须可以接受矢量输入
%fun2=[0::;sin(0::];%最多8个点,必须等距
%y1=NewtonCotes(fun1,0,,6)
%y2==NewtonCotes(fun2)
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和节点函数值fx
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';
functionCk=cotescoeff(n)