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

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

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

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

数值计算方法上机题目3

一、计算定积分的近似值：

要求：

（1）若用复化梯形公式和复化Simpson公式计算，要求误差限

，分别利用他们的余项估计对每种算法做出步长的事前估计；

（2）分别利用复化梯形公式和复化Simpson公式计算定积分；

（3）将计算结果与精确解比较，并比较两种算法的计算量。

1.复化梯形公式

程序：

程序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）

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=7.3891

等分数n=7019

已知值与计算值的误差R=2.8300e-008

2.Simpson公式

程序：

程序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）

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=7.3891

等分数n=24

已知值与计算值的误差R=2.7284e-008

用复化梯形公式计算的结果为：

7.3891，与精确解的误差为：

2.8300e-008。

等分数n=7019

用复化Simpson公式计算的结果为：

7.3891，与精确解的误差为：

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

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=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.1:

0.5）];%最多8个点，必须等距

%y1=NewtonCotes（fun1,0,0.5,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）

%由于科特斯系数最多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阶科特斯系数

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

functionf=intfun（t,n,k）

%科特斯系数中的积分表达式

f=1;

fori=[0:

k-1,k+1:

n]

f=f.*（t-i）;

end

输出结果：

fun=@（x）exp（x）;

a=-1;

b=1;

n=4;

NewtonCotes（fun,a,b,n）

ans=2.3505

二、三点数值微分

欢迎您的下载，

资料仅供参考！

致力为企业和个人提供合同协议，策划案计划书，学习资料等等

打造全网一站式需求