数值积分与微分方程

2.3 数值积分

2.3.1 一元函数的数值积分

函数1 quad 、quadl 、quad8

功能 数值定积分，自适应Simpleson 积分法。

格式 q = quad(fun,a,b) %近似地从a 到b 计算函数fun 的数值积分，误差为10-6。

若给fun 输入向量x ，应返回向量y ，即fun 是一单值函数。

q = quad(fun,a,b,tol) %用指定的绝对误差tol 代替缺省误差。tol 越大，函数计

算的次数越少，速度越快，但结果精度变小。

q = quad(fun,a,b,tol,trace,p1,p2,…) %将可选参数p1,p2,…等传递给函数

fun(x,p1,p2,…)，再作数值积分。若tol=[]或trace=[]，则用缺省值进行计算。

[q,n] = quad(fun,a,b,…) %同时返回函数计算的次数n

… = quadl(fun,a,b,…) %用高精度进行计算，效率可能比quad 更好。 … = quad8(fun,a,b,…) %该命令是将废弃的命令，用quadl 代替。 例2-40

>>fun = inline(‘3*x.^2./(x.^3-2*x.^2+3)’); equivalent to: function y=funn(x)

y=3*x.^2./(x.^3-2*x.^2+3);

>>Q1 = quad(fun,0,2) >>Q2 = quadl(fun,0,2)

计算结果为：

Q1 =

3.7224 Q2 =

3.7224

补充：复化simpson 积分法程序

程序名称 Simpson.m

调用格式 I=Simpson('f_name',a,b,n)

程序功能 用复化Simpson 公式求定积分值

输入变量 f_name 为用户自己编写给定函数()y f x 的M 函数而命名的程序文件名 a 为积分下限

b 为积分上限

n 为积分区间[,]a b 划分成小区间的等份数 输出变量 I 为定积分值 程序

function I=simpson(f_name,a,b,n) h=(b-a)/n; x=a+(0:n)*h; f=feval(f_name,x); N=length(f)-1;

if N==1

fprintf('Data has only one interval') return; end if N==2

I=h/3*(f(1)+4*f(2)+f(3)); return; end if N==3

I=3/8*h*(f(1)+3*f(2)+3*f(3)+f(4)); return; end I=0;

if 2*floor(N/2)==N

I=h/3*(2*f(N-2)+2*f(N-1)+4*f(N)+f(N+1)); m=N-3; else m=N; end

I=I+(h/3)*(f(1)+4*sum(f(2:2:m))+2*f(m+1)); if m>2

I=I+(h/3)*2*sum(f(3:2:m)); end

例题 求0

sin I xdx π

=

⎰

。

解 先编制sin y x =的M 函数。程序文件命名为sin_x.m 。

function y=sin_x(x) y=sin(x)

将区间4等份，调用格式为

I=Simpson (’sin _x’,0,pi,4)

计算结果为

y =

0 0.7071 1.0000 0.7071 0.0000

I =

2.0046

将区间20等份，调用格式为

I=Simpson (’sin _x’,0,pi,20)

计算结果为

y =

0 0.1564 0.3090 0.4540 0.5878 0.7071 0.8090

0.8910 0.9511 0.9877 1.0000 0.9877 0.9511 0.8910 0.8090 0.7071 0.5878 0.4540 0.3090 0.1564 0.0000

I =

2.0000

重做上例2—40：simpson('funn',0,2,100)

函数2 trapz

功能 梯形法数值积分

格式 T = trapz(Y) %用等距梯形法近似计算Y 的积分。若Y 是一向量，则trapz(Y)

为Y 的积分；若Y 是一矩阵，则trapz(Y)为Y 的每一列的积分；若Y 是一多维阵列，则trapz(Y)沿着Y 的第一个非单元集的方向进行计算。

T = trapz(X,Y) %用梯形法计算Y 在X 点上的积分。若X 为一列向量，Y 为

矩阵，且size(Y,1) = length(X)，则trapz(X,Y)通过Y 的第一个非单元集方向进行计算。

T = trapz(…,dim) %沿着dim 指定的方向对Y 进行积分。若参量中包含X ，则

应有length(X)=size(Y ,dim)。

例2-41

>>X = -1:.1:1;

>>Y = 1./(1+25*X.^2); >>T = trapz(X,Y)

计算结果为：

T =

0.5492

补充： 复化梯形积分法程序

程序名称 Trapezd.m

调用格式 I=Trapezd('f_name',a,b,n) 程序功能 用复化梯形公式求定积分值

输入变量 f_name 为用户自己编写给定函数()y f x 的M 函数而命名的程序文件名 a 为积分下限

b 为积分上限

n 为积分区间[,]a b 划分成小区间的等份数 输出变量 I 为定积分值 程序

function I=Trapezd(f_name,a,b,n) h=(b-a)/n;