# Trapezoidal rule to approximate integrals

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 $[a,b]$ be an interval and $f:[a,b]\to \mathbb{R}$ a real function. We want to approximate the integral $\int_a^bf(u) \; du$.

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

$\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)$

that can be written as

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

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

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

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

$\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.$

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

$\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.$

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

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

 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