MATLAB: What are some ways to use trapz,integral or quad in reverse

differential equationsnumerical integration

I have used ode45 to solve two simultaneous differential equations and plotted their solutions, such that : Area=trapz(x) or Area=integral(x), where 'x' is a function of 'phi'
My intention is to define the 'Area' first hand and find the associated 'phi'.
Example:
    Let, Area of plot=trapz(x).
          x=func(phi)
           I give the Area and I want to find the associated 'phi'.
Another way to put it is:
If, V=trapz(x) for V=1, stop/print trapz(x)

Best Answer

It will essentially never happen that area is EXACTLY the value you want to find.
You could of course trivially use cumtrapz, to compute a cumulative trapezoidal integral. Then use find to find the first x where the area exceeds the target value.
x = linspace(0,2*pi,10);
y = cos(x);
Of course, we know that the integral is sin(x). But can we find the point where the integral is exactly 0.85?
Iy = cumtrapz(x,y)
Iy =
          0      0.61647      0.94448      0.83056      0.32801     -0.32801     -0.83056     -0.94448     -0.61647  -4.4409e-16
It seems to lie between the second and third points, although there are several locations where it would happen. Find could find at least the first location where that value is exceeded.
ind = find(Iy > 0.85)
ind =
     3
x(ind)
ans =
       1.3963
Pretty crappy, really, in terms of the true value, which should arise at:
asin(0.85)
ans =
        1.016
You point out that you understand tools like ODE45. Note that trapz and cumtrapz are equivalent to tools like Euler's method, so relatively poor ways to solve an ODE.
Could you write a cumulative trapezoidal rule that solves for the exact point where the trapezoidal integral hits some target value? Well, yes. It would still be terribly poor in terms of accuracy. Were I to want to do that (god knows why) I would probably write it as a piecewise linear spline. Then I would integrate the spline, and finally solve for the location(s) of interest.
But, if I was going to do it using a spline, just use a real cubic spline in the first place.
spl = spline(x,y);
Ispl = fnint(spl);
fnint comes with the curve fitting toolbox as I recall. If you lack that TB, it is easily enough written though.
Next, solve for the location of interest. In my case, I would just be lazy, and use the tools I've already written. Here, slmsolve, part of my SLM toolbox.
slmsolve(Ispl,0.85)
ans =
     1.0126       2.1286
It found both locations where that event takes place.
But it is quite easy to solve using just fzero too.
target = 0.85;
fzero(@(x) ppval(Ispl,x) - target,1)
ans =
       1.0126
fzero(@(x) ppval(Ispl,x) - target,2)
ans =
       2.1286
Finally, could I have used an ODE solver, like ODE45 to solve this? Well, yes, but WHY would I? But yes, you could formulate an ODE to solve this problem using events.
Related Question