Trapezoidal rule to approximate integrals

De MateWiki
Saltar a: navegación, buscar

In this article we focus on the implementation of the numerical approximations of integrals by the trapezoidal rule in one and two dimensions. We refer to Trapezoidal Rule for a deduction of the formula.

  1. One dimensional integrals

Let [math] [a,b] [/math] be an interval and [math] f:[a,b]\to \mathbb{R} [/math] a real function. We want to approximate the integral [math] \int_a^bf(u) \; du [/math].

Consider a partition of the interval [math] [a,b] [/math] in [math] N [/math] equal subintervals of length [math] h=\frac{b-a}{N} [/math], given by [math] u_n=a+nh, [/math] where [math] n=0,1,...,N [/math]. The trapezoidal rule is as follows:

[math] \int_a^b f(u) \; du \sim h\frac12 f(u_0)+h\sum_{i=1}^{N-1}f(u_i)+h\frac12 f(u_N) [/math]

that can be written as

[math] \int_a^b f(u) \; du \sim h\sum_{i=0}^{N}w_if(u_i)=hw \cdot f, [/math]

where [math] w_i [/math] are the components of the weight column vector [math] w=(1/2,1,1,...,1,1,1/2)^t [/math] and [math] f [/math] is the column vector [math] f=(f(u_0),f(u_1),...,f(u_N))^t [/math].

Example: In Matlab/Octave we approximate the integral of the function [math] e^{-x^2} [/math] in the interval [math] [-1,1] [/math] with [math] h=0.01 [/math]

1 N=200;                         %Number of points
2 a=-1; b=1;                     %Extremes of the interval
3 h=(b-a)/N;
4 u=a:h:b;                       %coordinates of the partition
5 f=exp(-u.^2)';                 %function
6 w=ones(N+1,1);                 %weights vector
7 w(1)=1/2; w(N+1)=1/2;
8 result=h*w'*f                  % result

  1. Two dimensional integrals

Let [math] [a,b]\times [c,d] [/math] be a rectangle and [math] f:[a,b]\times[c,d]\to \mathbb{R} [/math] a real function. To approximate the integral [math] \int_a^b\int_c^df(u,v) \; dv \; du [/math], we apply the trapeoidal rule iteratively. First, we consider a partition of the interval [math] [a,b] [/math] in [math] N_1 [/math] equal subintervals of length [math] h_1=\frac{b-a}{N_1}[/math]. Define [math] u_n=a+nh, [/math] where [math] n=0,1,...,N_1 [/math]. The trapezoidal rule gives us:

[math] \int_a^b \int_c^d f(u,v) \; du \; dv \sim h_1\sum_{n=0}^{N_1}w_{n}\int_c^df(u_n,v)dv. [/math]

Now, we use again the trapezoidal rule for the remaining integral with a partition of the interval [math] [c,d] [/math] in [math] N_2 [/math] equal subintervals of length [math] h_2=\frac{d-c}{N_2}[/math],

[math] \int_a^b \int_c^d f(u,v) \; du \; dv \sim h_1h_2\sum_{n=0}^{N_1}\sum_{m=0}^{N_2}w_{n}\hat w_mf(u_n,v_m)dv =h_1h_2w\cdot f\cdot \hat w. [/math]

where [math] f [/math] is the [math] (N_1+1)\times (N_2+1) [/math] matrix with components [math] f_{ij}=f(u_i,v_j)[/math] and [math] \hat w [/math] is similar to [math] w [/math] but with [math] N_2+1 [/math] rows.

Example: In Matlab/Octave we approximate the integral of the function [math] f(u,v)=e^{-u^2+v} [/math] in the interval [math] [-1,1]\times [0,1] [/math] with [math] h_1=h_2=0.01 [/math]

 1 N1=200; N2=100;                  %Number of points
 2 a=-1; b=1; c=0; d=1;             %Extremes of the interval
 3 h1=(b-a)/N1; h2=(d-c)/N2;
 4 u=a:h1:b; v=c:h2:d;              %coordinates of the partition
 5 [uu,vv]=meshgrid(u,v);           %coordinates of the rectangle
 6 f=exp(-uu.^2+vv);                %function
 7 w1=ones(N1+1,1);                 %weights vector
 8 w1(1)=1/2; w1(N1+1)=1/2;
 9 w2=ones(N2+1,1);                 %weights vector
10 w2(1)=1/2; w2(N2+1)=1/2;
11 result=h1*h2*w2'*f*w1            % result