Integration Limits on Trapz

Hi Everyone!

I have a very simple question. I used MATLAB to solve a set of coupled ODEs from zero to Inf. And now I have the vectors containing the solution for these functions and their first derivatives in this semi-infinite interval.

My question is: How can I ask trapz to integrate a functional (depending on these solutions and their derivatives) from xmin=1e-5 to xmax=100, only? Is it enough to insert extra inputs into trapz?

Thank you,

MdL.

Answer 1

[Walter Roberson]

x = sort(rand(1,20));
y = exp(sin(2*pi*x));
yint = cumtrapz(x, y);
plot(x, y, 'k', x, yint, 'b')

syms X B
Y = exp(sin(2*pi*X));
Yint = int(Y, X, 0, X);
fplot([Y, Yint], [0 1])

This illustrates that you can use trapz() with irregular spacing, and that if the values are close enough, it will approximate the true integral.

If you only want to integrate over a particular range, then find the endpoints of the range in the data and only do that portion:

mask = pi/10 <= x & x <= 2*pi/10;
px = x(mask);
py = y(mask);
pyint = cumtrapz(px, py);
plot(px, py, 'k', px, pyint, 'b')

The reason this does not look the same is that you are not bringing in the accumulated value from 0 to pi/10, and definite integration always starts at 0.

[MLPhysics]

Hi @Walter Roberson! I still have one doubt... In your last example (pretty similar to what I need to do) you have used cumtrapz, right? In my case, I just need to change it to trapz instead? I mean, I want to integrate in a finite interval a specific expression which comes from my numerical solution.

[MLPhysics]

@Walter Roberson And thank you for the answer!

[Walter Roberson]

The only difference between cumtrapz() and trapz() is that cumtrapz() gives you the intermediate results as well, such as you might want for plotting.

I could have just asked it to compute trapz() and int() at the endpoint, and given you numeric differences, but the meaning of the differences would not have been clear at all; it was better for illustration purposes to cumtrapz() so the progress towards the end could be plotted.

[MLPhysics]

All right, I get it. Thank you very much for your explanation.

[MLPhysics]

@Walter Roberson After I went back to my MATLAB code, I noticed there is a tiny detail I don't know how to deal with. In my case, numerical solutions are stored in the following arrays:

s.x1tau -> the whole set of points where functions are evaluated.

s.valx1tau(n,:) -> values of the functions for each point (n=1,2,3 stands for these 3 functions).

yderiv(n,:) -> values of the first derivative of the above functions for each point in the interval

I want to integrate something like this:

f=trapz(s.x1tau,(1/4)*(1/s.x1tau)*yderiv(3,:).^2. + (1/2)*s.x1tau*(yderiv(1,:).^2. + yderiv(2,:).^2.) + ...

(1/4)*(1/s.x1tau)*s.valx1tau(2,:).^2.*(1-s.valx1tau(3,:)).^2.)

So, I guess I should use a strategy similar to your last example in the answer above. However, I don't how to write that I only want to take a finite interval from the set of all points. I mean, should I substitute the ":" inside the arrays?

[Walter Roberson]

In the case where the subset can be found knowing only s.x1tau, and assuming that the subset is continguous:

LB = pi/10; UB = 2*pi/10;     %USE YOUR OWN VALUES
mask = LB <= s.x1tau & s.x1tau <= UB;
ssx1tau = s.x1tau(mask);
f = trapz(ssx1tau, (1/4) * (1./ssx1tau) .* yderiv(3,mask).^2 + (1/2) * ssx1tau .* (yderiv(1,mask).^2 + yderiv(2,mask).^2) + ...
	(1/4) * (1./ssx1tau) .* s.valx1tau(2,mask).^2 .* (1-s.valx1tau(3,mask)).^2);

[MLPhysics]

Thank you! Thank you! I was just thinking about redefining s.x1tau as you did, but I didn't think about substituting ':' by 'mask'. It helped me a lot. :D

[Walter Roberson]

Watch out for the places you had 1/s.sx1tau : s.sx1tau is a vector so 1/ it is matrix division . Notice I corrected them to 1./

The .* I put in are important too.

[MLPhysics]

Yes, I forgot about that. Thanks for remembering me. Now the calculation worked fine.