Get integrated vector values

Question

I am trying to integrate the sine-function. My goal is to get not just the value of the area inbetween a certain distance but the specific values of the integrated course.

One way to achieve this is by using cumtrapz. I want to get the same result using integral or quad. So I am wondering if there is something like cumquad?

I tried to write something for myself but it works very slow and seems to be even worse than cumtrapz. Later on I want to integrate measured data. So it won't be as simple as a sine.

Here is my current code:

a = 0; b = 10;

x = a:0.1:b; 
y = 2*sin(3*x);
pp = spline(x,y);

y2=zeros(1,length(y));
y3=zeros(1,length(y));

y2(1)=integral(@(x)ppval(pp,x),x(1),x(2));
y3(1)=integral(@(x)ppval(pp,x),x(1),x(2));
for a=2:(length(y)-1)

   y2(a) = y2(a-1)+integral(@(x)ppval(pp,x),x(a-1),x(a));
   y3(a) = y3(a-1)+quad(@(x)ppval(pp,x),x(a-1),x(a));

end
y4=cumtrapz(x,y);
% y5=cumsum(y);

plot(x,y)
hold on
plot(x,y2,'-ro')
plot(x,y3,'-kx')
plot(x,y4,'g')
syms x % compare with analytical result
ya=2*sin(3*x);
ya5=int(ya)+(2/3);
ezplot(x,ya5)

Answer 1

Using integral

I don't think there is a way to have MATLAB return the integral along the path, so you are correct in performing the integration one Δx at a time. The slowness comes from the loop and subsequent restart of every integral call. You can avoid the loop by posing the integral over every interval as a vector-valued function.

The Math

Suppose we divide x into N-1 intervals with N total boundaries and denote an interval boundary as x_n where n ∈{1,2,3...,N} such that x₁ ≤ x₂ ≤ x₃ ... ≤ x_N. Then any integral over the interval would be

Using the u-substitution:

The integral becomes:

where Δx_n = x_n - x_n-1

The Code

So now, we can pose the interval integration of any function by specifying the lower bound x_n-1, specifying the interval width Δx, and integrating from 0 to 1. The best part is that if the lower bound and interval widths are vectors, we can create a vector-valued function in terms of u and have integral integrate with the option 'ArrayValued' = true.

x    = a:0.1:b; 
xnm1 = x(1:end-1);
dx   = x(2:end) - xnm1;
fx   = @(x) 2*sin(3*x);
f    = @(u) dx .* fx(dx*u+xnm1);
y    = cumsum([0,integral(@(u)f(u),0,1,'ArrayValued',true)]);

The cumsum accounts for the fact that each integral over a given interval needs to have the value of the previous interval added to it.

On my machine, this is at least order of magnitude faster than the loop version and gets better as the interval count increases.

Using ode45

Use can also use ode45 to perform the integration. It is not nearly as efficient as the integral method, but it may be easier conceptually and look cleaner. In fact, ode45 is about 10 times slower than the integral method above when required to return an absolute error on par with that of integral.

a    = 0;
b    = 10;

% These options are necessary to approach the accuracy of integral
opt = odeset('RelTol',100*eps(),'AbsTol',eps());
sol = ode45(@(x,y) 2*sin(3*x),[a,b],0,opt);

x    = a:0.01:b; 
yint = deval(sol,x);