Trapezoidal rule to approximate integrals
In this article we approximate integrals by the trapezoidal rule.
- 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^t \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
- 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. We want to approximate the integral [math] \int_a^b\int_c^df(u) \; dv \; du [/math].
We apply the trapeoidal rule to both integrals iteratively. First 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}w_{n}\int_c^df(u_n,v)dv. [/math]
Now, we use the formula for the remaining integral,
- [math] \int_a^b \int_c^d f(u,v) \; du \; dv \sim h_1h_2\sum_{n=0}^{N}\sum_{m=0}^{N}w_{n}\hat w_mf(u_n,v_m)dv =h_1h_2w^t\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*w1'*f*w2 % result
