第五章数值积分-席

时间：2025-07-07

第五章

数值积分

对定积分进行MATLAB符号计算可以使用命令： f=int(s,a,b) 输入量：s可以是被积函数，也可以是几个被积函数组成的矩阵；积分上下限a,b可以是数也可以是符号变量；积分变量默认为x,若积分变量为t,采用 f=int(s,t,a,b)

例:计算解：syms t f=5/((t-1)*(t-2)*(t-3));F=int(f,t,4,5),y=double(F) 例：

5 4 t 1 t 2 t 3 dt5

ax b dx n cx d

解： syms x a b c d n f=(a*x+b)/((c*x+d)^(1/n)); F=int(f,0,1),y=simple(F)

在求定积分 I f ( x)dx 时，a

问题：sin t dt 等； (1)积不出：如 a t (2)难积：当 f ( x)比较复杂时，不好使用N-L公式；b

(3)表格函数xi f ( xi ) x0 f ( x0 ) x1 xn f ( xn ) f ( x1 )

一、梯形法对积分

I f ( x)dxa

用线性插值函数 g ( x ) 近似代替被积函数 f ( x) , 其中： b x x a g ( x) f1 f2 b a b a f1 f (a) f 2 f (b) 则h I f ( x)dx g ( x)dx ( f1 f 2 ) a a 2 上式称为梯形积分公式。b b

h b a

也可写为：h I f ( x)dx ( f1 f 2 ) E a 2 其中 E 1 h3 f 为截断误差 12b

yf1y f ( x)

y f ( x)f3

y g ( x)f2

o x1 a

x2 b x o

复合梯形积分法：

若被积函数在等距的 n+1 个点上取值，可在每个区间段逐个利用梯形公式，整理后得： b h h h a f ( x)dx 2 ( f0 f1 ) 2 ( f1 f 2 ) 2 ( f n f n 1 ) E n h f 0 2 f k f n 1 +E 2 k 1

b a 其中： h , xi a (i 1)h, f i f ( xi ) i 1,2, , n 1 n b a 2 (b a)3 误差 E h f f 2 12 12n

f 为f "( x)在a x b上的平均值

Matlab里提供了复合梯形公式的计算程序— trapz.m,其调用格式为： Z=trapz(y) y为向量，则Z是y的积分，若y为矩阵，则行向量Z是y的每一列的积分； Z=trapz(x,y) 输入x、y为同长度的数组，计算y对x积分。例： >>X = 0:pi/100:pi; Y = sin(X); Z = trapz(X,Y) Z0 = pi/100*trapz(Y)

Z= 1.9998 Z0 = 1.9998 >>Y=[1 2 3 4]; trapz(Y) ans = 7.5000 根据上面的复合梯形积分公式： trapz(Y)=[1+2(2+3)+4]/2

n h f ( x)dx f 0 2 f k f n 1 +E 2 k 1

b a 其中： h , xi a (i 1)h, f i f ( xi ) i 1,2, , n 1 n

基于梯形法的积分函数文件有trapez_v,trapez_n和 trapez_g,分别如下： n b h a f ( x)dx 2 f0 2 fk f n 1 E k 1 function I = trapez_v(f, h) I = h*(sum(f) -(f(1) + f(length(f)))/2);%复合梯形公式 function I = trapez_n(f_name, a, b, n) h = (b-a)/n;

x = a+(0:n)*h;f = feval(f_name, x);%求f_name(x),f_name被积函数

%文件名:feval(' fun ',x)或feval(@fun,x)I = trapez_v(f,

function I = trapez_g(f_name, a, b, n) b h h = (b-a)/n; I f ( x)dx ( f1 f 2 ) E a x = a+(0:n)*h; f = feval(f_name, x); 2

I = h/2*(f(1) + f(n+1)); %n=1,积分区间没有进行分割 if n>1, I = I + h*sum(f(2:n));end n b h I; f ( x)dx f 0 2 f k f n 1 +E a 2 k 1 h2 = (b-a)/100; xc = a+(0:100)*h2; fc = feval(f_name, xc); plot(xc,fc,'r'); hold on title('Trapezoidal Rule');xlabel('x');ylabel('y');

plot(x,f, 'ko');plot(x,zeros(size(x))) %x轴上分点图形 for i=1:n; plot([x(i),x(i)], [0,f(i)]); end

例1 已知积分精确值 I 4.006994 ,分析用梯形法计算下述积分时，分割区间数 n 对误差的影响。

1 exp( x)dx

clear; Iexact = 4.006994;a = 0; b=2; n = 1; fprintf('\n n I Error\n');

for k=1:6n = 2*n; h = (b-a)/n; i = 1:n+1; %10.5f\n', n, I, Iexact - I); x = a + (i-1)*h; f = sqrt(1 + exp(x)); I = trapez_v(f,h); fprintf(' %3.0f %10.5f end

运行结果：n 2 4 8 16 32 64 I 4.08358 4.02619 4.01180 4.00819 4.00729 4.00707 Error -0.07659 -0.01919 -0.00480 -0.00120 -0.00030 -0.00008

结果表明：复合梯形公式的误差随着 n 的加倍减少为原来的1/4。

调用trapez_g函数：第一步建立被积函数文件： function y=f(x) y=sqrt(1 + exp(x)) 再调用： I = trapez_g(@f, 0, 2, 16)Trapezoidal Rule 3

2.5

1.5

0.5

0.2

0.4

0.6

0.8

1 x

1.2