复化梯形公式和复化Simpson公式Word文档下载推荐.docx

复化梯形公式和复化Simpson公式Word文档下载推荐.docx

《复化梯形公式和复化Simpson公式Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《复化梯形公式和复化Simpson公式Word文档下载推荐.docx（9页珍藏版）》请在冰豆网上搜索。

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））'

）%定义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）

Tn=h/2*（（f（a）+2*Tn1+f（b）））

z=exp

（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（x）的n阶导数）：

输入n=4

4*exp（x）+x*exp（x）

（4*exp（x）+x*exp（x））'

）%定义f（x）的四阶导数，输入程序1里求出的f2即可

-（4*exp（x）+x*exp（x））'

%因fminbnd（）函数求的是表达式的最小值，且要求表达式带引号，故取负号，一边求最大值

e=5*10^（-8）%精度要求值

x1=fminbnd（f3,1,2）%求负的四阶导数的最小值点，也就是求四阶导数的最大值点对应的x值

Rn=-（b-a）/180*（（b-a）/（2*n））^4*f2（x1）%计算余项

Sn1=0

Sn2=0

fork=0:

n-1%求两组连加和

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）

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构成应为：

第一列为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

fun积分函数的句柄，且必须能够接受矢量输入！

else

输入参数错误，请参考函数帮助！

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）;

（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=@（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:

0.1:

0.5;

sin（0:

0.5）];

%最多8个点，必须等距

%y1=NewtonCotes（fun1,0,0.5,6）

%y2==NewtonCotes（fun2）

%计算积分节点xk和节点函数值fx

xk=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=@（x）exp（x）;

a=-1;

b=1;

n=4;

NewtonCotes（fun,a,b,n）

ans=2.3505

二、三点数值微分