圆心距离公式复化梯形公式和复化Simpson公式

2020年10月2日发(作者：徐霞客)
数值计算方法上机题目3
一、计算定积分得近似值：

要求：
（1）若用复化梯形公式与复化Simpson公式计算,要求误差限，分
别利用她们得余项估计对每种 算法做出步长得事前估计；
（2）分别利用复化梯形公式与复化Sｉmpｓon公式计算定积分；
（3）将计算结果与精确解比较，并比较两种算法得计算量.
1、复化梯形公式
程序：
程序1（求ｆ(x）得ｎ阶导数:
ｓｙms x
ｆ=x＊exp（x) %定义函数f(x)
n=inｐｕt('输入所求导数阶数:’）
f2=ｄiff（f,x,n） %求f(x）得ｎ阶导数
结果1
输入n=2
f２ =
２*ｅxp（ｘ) + x＊exp（x)
程序2：
cｌc
cｌear
ｓyｍs ｘ %定义自变量ｘ
ｆ=ｉnlｉnｅ（'x＊exp(x）＇，'x’) ％定义函数f(x)=x*exp（x）,换函数时只需换
该函数表达式即可
ｆ２=inlinｅ（’（2＊exｐ(x） + x*exp(x))'，'x’） ％定义f（ｘ)得二阶导数，输
入程序1里求出得f2即可。
ｆ3='－（2＊ｅxｐ（x） ＋ ｘ*exp（x））＇ %因ｆminbnd()函数求得就是表达式
得最小值，且要求表达式带引号,故取负号，以便求最大值
ｅ=5＊1０^（—8） ％精度要求值
a＝１ %积分下限
b=2 ％积分上限
x1=fminbｎd（f3,1，２) %求负得二阶导数得最小值点，也就就是求二阶导数得
最大值点对应得x值
foｒ n＝2：１0000００ ％求等分数n
Rn＝—(b—a)12＊(（b－a）n)^2＊f2(x1) %计算余项
if abs（Rn） brｅaｋ ％ 符合要求时结束
eｎd
ｅnd
h＝（b－a）n %求ｈ
Tn1＝0
fｏr k=1：n—1 ％求连加与
ｘｋ=a+k＊h
Ｔｎ1=Tｎ1+ｆ（xk）
end
Ｔｎ=ｈ2＊((f(a）+2*Tｎ1+f（b)))
ｚ=ｅxp（2)
R=Tn-z %求已知值与计算值得差
fpｒintf（＇用复化梯形算法计算得结果 Tn=＇）
ｄiｓｐ（Ｔn）
fｐrintｆ('等分数 n=’)
disp（ｎ) %输出等分数
fprintf('已知值与计算值得误差 R＝'）
dｉsｐ（Ｒ)
输出结果显示:
用复化梯形算法计算得结果 Tｎ= 7、3８91
等分数 n=7０１９
已知值与计算值得误差 Ｒ= 2、８300e—00８
2、 Simpson公式
程序:
程序1：（求f（x)得n阶导数）：
syｍs x

f=x＊eｘｐ（x） ％定义函数f(x）

ｎ=inｐuｔ（'输入所求导数阶数：＇）

f2=diff（f，x，n） ％求f(ｘ）得n阶导数
结果1
输入n＝4
f2 =
4＊eｘp（ｘ) + x*exp(x)
程序2:
clc
cleaｒ
sｙmｓ x ％定义自变量ｘ
ｆ=inline(’x＊exp(ｘ）’,＇x’） ％定义函数f(x）=ｘ＊exp（ｘ），换函数时只需
换该函数表达式即可
f2=ｉnline('(４＊eｘp（x) + x*ｅxp(x）)’,’x'） ％定义f(x）得四阶导数，输
入程序1里求出得f２即可
f3='—（４＊exp（x) + x＊ｅxp（x））' %因ｆmiｎbnd（)函数求得就是表达式
得最小值,且要求表达式带引号，故取负号,一边求最大值
e=５*１０^（-８) ％精度要求值
a＝1 %积分下限
b＝2 ％积分上限
ｘ１=fｍiｎbnd(f3，1,2) ％求负得四阶导数得最小值点，也就就是求
四阶导数得最大值点对应得x值
for n=２:１0０0000 ％求等分数n
Rn＝—（b—a）1８0*（(b-ａ）(2＊n）)^４＊ｆ2(x１) ％计算余项
if aｂs(Rn)〈e ％用余项进行判断
break ％ 符合要求时结束
ｅnd
ｅnd
h=（b－ａ）ｎ ％求ｈ
Ｓn1=０
Ｓn2=0
for k=０:ｎ－1 ％求两组连加与
xk=a+k＊h
ｘk1=ｘk+h2
Ｓn1=Sn1+f（xk1）
Sｎ2=Sn２+f（xk）
ｅnd
Sn=h6*（f(a）+4＊Sn1＋２*（Sn２-ｆ（a)）+ｆ(b)） ％因Ｓn2多加了k=0时得值,
故减去f（a)
z=ｅxp(２)
R=Sn—z %求已知值与计算值得差
fpｒintf(’用Simpsoｎ公式计算得结果 Sn＝＇）
dｉsp（Sn)
fpｒｉntｆ('等分数 n＝’）
disp(n)
fprｉntf('已知值与计算值得误差 R=’）
disp（Ｒ)
输出结果显示:
用Siｍpsoｎ公式计算得结果 Sn= 7、3891
等分数 n=２４
已知值与计算值得误差 R＝ 2、7２84ｅ—008
用复化梯形公式计算得 结果为:7、3８91，与精确解得误差为:2、８3
０0ｅ-00８。等分数n＝70１9
用复化Simpｓon公式计算得结果为：7、38９1，与精确解得误差
为:2、7２8４e-008 .等分数ｎ＝2４
3、柯斯特公式求积分：
程序代码：
(1)fuｎctiｏｎ ［y，Ck，Ak］=NewtｏnCotes（fun,a，b,n）
if nａrgin==1
［mm，nn］=size（fun);
ｉｆ mm〉＝8
erroｒ（'为了保证NeｗtoｎCｏtes积分得稳定性,
最多只能有9个等距节点！'）
ｅｌseｉf nn~=2
errｏr（’fｕｎ构成应为：第一列为x，第二列为y，
并且个数为小于10得等距节点!'）
end
xk=fun（1,：);
fk=fun（2,:)；
a＝min（xk)；
b＝max（xk）；
n＝mm－1；
ｅlｓｅif nargiｎ==4
xk=ｌinｓpａce（a,b，ｎ+1）;
ｉｆ isa(fuｎ，’fuｎｃtiｏn＿hanｄle')
fｘ=fuｎ(xk）；
elｓｅ
eｒror（'fuｎ积分函数得句柄,且必须能够接受矢量输
入！')
eｎｄ
ｅlsｅ
erｒoｒ（’输入参数错误,请参考函数帮助!’)
end
Ck＝cｏtｅscoefｆ(ｎ);
Ak=（ｂ-a)＊Ｃk；
y＝Ak*fx＇;
（２)fuｎｃtion Ｃk=cｏｔescoeff(n）
fｏr ｉ=１：n+１
k=i-1;
Ｃｋ(i)=(—１)＾（n -k)faｃtorial（k)factorial（ｎ—k)n*ｑ
uadl(（t）ｉntfｕn (t，n，ｋ),0,n);
end
（3)ｆuncｔioｎ f=intｆｕｎ（t，n，k）
ｆ=1；
for i=［0:k—1,k+1：n］
f＝ｆ、*(t-i）;
end
代码解释:
fｕnｃtion ［y,Ck，Ak]=ＮewｔonCoｔｅｓ(ｆun,a,b，n)
％ y＝NewtoｎCotｅs(fｕn,a，ｂ,n)
% 牛顿－科特斯数值积分公式
% 参数说明:
％ ｆｕn，积分表达式,这里有两种选择
％（1）积分函数 句柄，必须能够接受矢量输入，比如ｆun＝（x）
sin（x）、*coｓ（x)
% （2)x，ｙ坐标得离散点， 第一列为ｘ, 第二列为y, 必须等距, 且
节点得个数小于9, 比如： fun＝[1:8；sｉn（1：8)］'
% 如果fｕn得表采用第二种方式，那么只需要输入第一个参数即可,
否则还要输入ａ，b，ｎ三个参数
% ａ，积分下限
% b，积分上限
％ ｎ，牛顿- 科特斯数公式得阶数，必须满足1〈ｎ〈7，因为ｎ〉
＝8时不能保证公式得稳定性
％ （１）n＝1,即梯形公式
% (2）n＝2，即辛普森公式
% （３）n=４，即科特斯公式
％ y，数值积分结果
％ Ck,科特斯系数
％ Ak，求积系数
％
％ Ｅxamｐle
% fｕn1=（x)sin（ｘ);％必须可以接受矢量输入
% fｕn2＝[0：0、1:０、5;sin（0：０、1:0、５）］;%最多8个点，
必须等距
％ y1=NewtoｎCoteｓ(fun1，0，0、５，6）
% y２=＝ＮeｗtonCotｅs（fun2)
ｉｆ nａrgｉn==1
[mm，nｎ]=sｉzｅ(fuｎ);
ｉf mm〉=8
eｒrｏｒ('为了保证NewtonCｏtes积分得稳定性,最
多只能有９个等距节点!')
eｌseif nn~=2
eｒror（’fun构成应为:第一列为x,第二列为y,并且
个数为小于10得等距节点！＇)
end
ｘk=fun(1,：）；
ｆｋ＝fｕn(2,：）;
a=miｎ(xk)；
b＝max（xｋ）;
n=mm－1;
elseｉf ｎaｒｇiｎ==4
％ 计算积分节点xk与节点函数值fx
xｋ=linspace（ａ，ｂ，n+１);
ｉf isａ(fｕn，’fｕncｔiｏn_handlｅ’）
fx＝ｆｕn（ｘk）;
elsｅ
eｒroｒ（＇fun积分函数得句柄,且必须能够接受矢
量输入！’）
end
else
erroｒ（＇输入参数错误，请参考函数帮助！')
end
% 计算科特斯系数
Ｃk=coｔeｓcoｅff（n)；
% 计算求积系数
Ａk=(ｂ-ａ）＊Ｃk；
％ 求与算积分
y=Ak＊fx’;
fuｎctｉon Ｃk=ｃotescoｅfｆ(n)
％ 由于科特斯系数最多7阶，为了方便我们可以直接使用，省
得每次都计算
％ A1=［1,1]2
% A2=［1,4，１]６
% Ａ3=［1，３，3,1］８
％ A４=［7，３２，12,32，1］／９0
% Ａ5＝[1９,7５，50，５0,75，１９］288
％ A６=[41,2１6，27,272，27,216，41］840
% A7=[751,３57 ７，1３2３，2９8９，2989，1323，357７,
７5１］17２8０
％ 当时为了体现公式，我们使用程序计算n阶科特斯系数
for i=1:n+1
k=i-１；
Cｋ(i）＝（-１）＾（n－ｋ)／fａｃｔoriａl(k)factｏ rｉ
al(n—k)n*qｕadl(（t)ｉntfuｎ（ｔ,n,k），0，n);
ｅｎｄ
funｃtion f=iｎtfuｎ（ｔ,ｎ，k）
％ 科特斯系数中得积分表达式
ｆ=1；
ｆｏｒ i=[0：k－1，k+１：n］
f＝ｆ、＊（t—i）；
end
输出结果:
fun=（x)exp(x）;
a=-1;
b=１；
n＝4;
NewtonＣoteｓ（fun，a，b，ｎ)
ans =2、３505
二、三点数值微分