MATLAB数值积分法电子版

作者:凯鲁嘎吉 – 果壳网
http://www.cnblogs.com/kailugaji/

一、实验目的

诸多工程技术和数学钻探中要用到定积分,假诺不可能直接算不出精确值(如含在积分方程中的积分)或总结困难但可用近似值近似时,就用数值积分法方法加以解决。常用的算法有:复化梯形、辛甫生(辛普森(Simpson))、柯特斯(Cotes)求积法; 龙贝格(Romberg)算法;高斯(Gauss)算法。

二、实验原理

电子版 1

三、实验程序

下面给出复化Simpson求积法程序(梯形及柯特斯复化求积分程序可遵守编制):

电子版 2

四、实验内容

选一可精确算值的定积分,用复化的梯形法及复化Simpson求积法作近似总计,并相比较结实。

 

五、解答(按如下顺序提交电子版)

1.(程序)

xps.m:

function y=xps(x)
y=x^(3/2);

复化梯形公式:

   trap.m:

function [T,Y,esp]=trap(a,b,n)
h=(b-a)/n;
T=0;
for i=1:(n-1)
    x=a+h*i;
    T=T+xps(x);
end
T=h*(xps(a)+xps(b))/2+h*T;
syms x
Y=vpa(int(xps(x),x,a,b),8);
esp=abs(Y-T);

复化辛甫生(辛普森)公式:

   simpson.m:

function [SI,Y,esp]=simpson(a,b,m)
%a,b为区间左右端点,xps(x)为求积公式,m*2等分区间长度
h=(b-a)/(2*m);
SI0=xps(a)+xps(b);
SI1=0;
SI2=0;
for i=1:((2*m)-1)
    x=a+i*h;
    if mod(i,2)==0
        SI2=SI2+xps(x);
    else 
        SI1=SI1+xps(x);
    end
end
SI=h*(SI0+4*SI1+2*SI2)/3;
syms x
Y=vpa(int(xps(x),x,a,b),8);
esp=abs(Y-SI);

2.(运算结果)

>> [T,Y,esp]=trap(1,2,8)

T =

    1.8636


Y =

1.8627417


esp =

0.0008089288247354886607354274019599
>> [SI,Y,esp]=simpson(1,2,8)

SI =

    1.8627


Y =

1.8627417


esp =

0.000000020499792974248975951923057436943

从统计结果看:复化辛普森公式更纯粹。

3.(拓展(方法立异、体会等))

MATLAB中有部分平放函数,用于实施自适应求积分,都是基于Gander和Gautschi构造的算法编写的。

quad:使用辛普森(Simpson)求积,对于低精度或者不光滑函数功能更高

quadl:该函数使用了名叫洛巴托求积(Lobatto
Quadrature)的算法,对于高精度和光滑函数效率更高

拔取办法:

I=quad(func,a,b,tol);

func是被积函数,a,b是积分限,tol是期待的相对误差(即使不提供,默认为1e-6)

诸如对于函数f=xe^x在[0,3]上求积分,显明可以因此解析解知道结果是2e^3+1=41.171073846375336

先创立一个M文件xex.m

内容如下:

function f=xex(x)

f=x.*exp(x);

下一场调用:

>> format long

>> format compact

>> quad(@xex,0,3)

ans =

  41.171073850902332

可见有9位有效数字,精度如故蛮高的。

电子版,要是使用quadl:

>> quadl(@xex,0,3)

ans =

  41.171074668001779

反倒只有7位有效数字

相关文章