Diferencia entre revisiones de «Trapezoidal rule to approximate integrals»

Revisión del 15:29 1 dic 2014

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]

N=200;                         %Number of points
a=-1; b=1;                     %Extremes of the interval
h=(b-a)/N;
u=a:h:b;                       %coordinates of the partition
f=exp(-u.^2)';                 %function
w=ones(N+1,1);                 %weights vector
w(1)=1/2; w(N+1)=1/2;
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\times N_2 [/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]

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