Find integral upper limit knowing the result and having an array

Hello,

I am trying to find the upper limit of an integral, knowing the results, but I have an array instead of a function.

If I had a function, for example -2t, and an integral result of lets say -4, i could find the upper limit as shown below, using fzero.

fun = @(t)-2*t;

intfun = @(x)integral(fun,0,x) + 4;

xval = fzero(intfun,[0,10]);

and i find an xval=2 for this case.

Howover now I do not have a function, but an array, for exaple, y2 as shown below. I have different values for y2 as I find y2 solving a single degree of freedom system and it represents the velocity values for each time instance. I need to find the time tout (upper limit of integral) and lets say the results is 10. I tried to solve it as in the example above, however i get an error.

y2=0:0.1:15;

fun_1= y2;

vale1= @(tout)integral(fun_1,0,tout,'ArrayValued',true)+10;

ans2=fzero(vale1,[0,10]);

here is the error:

Error using fzero>localFirstFcnEval (line 729)

FZERO cannot continue because user-supplied function_handle ==> @(tout)integral(fun_1,0,tout,'ArrayValued',true)+10

failed with the error below.

First input argument must be a function handle.

Could you help me with it please?

Thank you!

Angie

[Torsten]

We need the values of t corresponding to the values of y2.

And for that the integral over y2 is 10, you need to subtract, not add 10 in your expression.

Or is y2 the negative of the function you want to integrate ?

[Angie]

Hi Torsen,

Regarding the 10, that's my mistake. It should be substracted indeed. The values of y2 i wrote are more of an example, but the corresponding values of t could be t=0:0.02:3.

Thank you!

Answer 1

[John D'Errico]

A basic rule is, you CANNOT use integral OR fzero with data, thus a vector of numbers.

You CAN approximate your data using a function, perhaps a spline interpolant, and then integrate that, and use fzero on the result. There is a fundamental difference here that you need to understand.

A second approach would be to form a cumulative integral, perhaps using cumtrapz. And then interpolate that.

For example, can you find the solution of int(sin(x)), on limits of [0,u], such that the integral is 1.75? So you would want to solve for u.

The problem is, I'll give you only data.

x = linspace(0,pi);
y = sin(x);

How would we solve it? A simple solution is as I said, to use cumtrapz.

yint = cumtrapz(x,y);
inttarget = 1.75;
plot(x,yint,'b-')
yline(inttarget,'r')

I''ve plotted the cumulative integral there, on the interval [0,pi].

The left endpoint of the integral is fixed at 0. The right endpoint is unknown. So all we need to do now is find the value of u, such that yint is exactly 1.5. Of course we will never be lucky enough that one of those data points will be exactly 1.75.

Here, since the function is monotone increasing on the interval, we can just swap x and y, and then use interp1.

format long g
u = interp1(yint,x,inttarget,'spline')

u = 
          2.41908043079086

How well did we do? We can check against the analytical result.

syms X U
Usol = solve(int(sin(X),X,[0,U]) == inttarget)

Usol = 

vpa(Usol)

ans = 

The second solution is the one we want. As you can see, the cumtrapz solution was not terrible. We got almost 4 correct digits in the solution.

The other approach may be a little more accurate, since it does not use a trapezoidal integration.

spl = spline(x,y); % Interpolate the curve first with a spline
splint = fnint(spl); % integrate
% Now solve for the upper limit of integration
[u,fval] = fzero(@(u) fnval(splint,u)-fnval(splint,0)-inttarget,[0,pi])

u = 
           2.4188584094742

fval = 
     0

This approach yields a few more significant digits. I would expect that to be the case as long as the data is not noisy. With noisy data, I would use other approaches that would probably involve smoothing the data first.

Answer 2

[Torsten]

t = 0:0.02:3;
y = 0:0.1:15;
fun = @(tt) interp1(t,y,tt,'spline');
tt0 = 1.0;
tt = fzero(@(x)integral(fun,0,x)-10,tt0)

tt = 2