Approximating derivative of numerical solution within event function

The issue I have is having to compute the derivative (in real time) of the solution produced by ode45 within the events function.

Some pseudo-code to explain what I'm mean is,

function dx = myfunc(t,x,xevent)
persistent xevent
% xevent is the solution at the event
dx(1) = x(2);
dx(2) = complicated function of x(1),x(2), and xevent;
end

function [value,isterminal,direction] = myeventfcn(t,x)
position = function of x(1), x(2), and dx(2); 
isterminal = 1;
direction = 0;
end

I know that if I didn't need to use the solution at the event within `myfunc` I could just compute dx=myfunc(t,x) within my event function and use `dx(2)`, yet since `xevent` is used within `myfunc` I can't input `xevent`.

I know there is a way of inputting constant parameters within the event function, but since it's the solution at the event location, which also changes, I'm not sure how to go about doing that.

My work around to this is to approximate `dx(2)` using the solution `x`. What I wanted to know is if it will be satisfactory to just use a finite difference approximation here, using a small fixed step size relative to the step size od45 takes, which is a variable step size.

Thanks for any help!

[Torsten]

To compute dx(2), you need "xevent" even before the event occured (i.e. starting from t=0) ?

In this case, you will have to solve your ODE several times in order to iteratively determine tevent and xevent.

Best wishes

Torsten.

[Shant Danielian]

I have xevent initialized at zero before the first nonzero event is detected, same with tevent. I already use the ODE event locator to find when tevent and xevent happen. I then stop the solver, gather the solution, reinput the new values/initial conditions into ode45 and continue until the next event. My problem is that I want a rough approximation of dx(2) within the event function, without having to call myfunc.

[Walter Roberson]

If you have duplicate calculations then put them into a function that is called by both myfunc and the event function.

[Shant Danielian]

What do you mean by duplicate calculations?

[Walter Roberson]

Calculations that are the same in both places.

Instead of trying to have myfunc estimate a value and push it out somewhere that the event function can read it, put the calculations into a routine and have both of them call that routine.

If you have R2017a or later, you could try memoize() to see if that helps -- it avoids recomputing values that have already been computed. However, the event function is not certain to be called with the same parameters as the ode function, and when ode45 is trying to find a place that is inside the boundary it can call the event function a number of times in a row looking for a place that is not rejected, which could be a problem for you if the calculation depends upon the derivative of places you have not yet evaluated.

ode45 is, by the way, permitted to move forwards and backwards in time or in boundary conditions. It mostly moves forward in time but when x is not a scalar then at times ode45 will call more than once with the same time and different boundary conditions.

[Shant Danielian]

I'm a little confused on what you mean by "put the calculations into a routine and have both of them call that routine". Do you mean have the myeventfcn and myfunc be nested in the same file?

I looked into memoize() and it's not entirely clear how to use it, but I'll mess with it to see where it goes.

I was wondering if I can use the options function to input the new values of xevent into the event function (and in turn use dx=myfunc(t,x,xevent) within the event function) after the solver detects the event and stops to gather the solution.

[Walter Roberson]

"Do you mean have the myeventfcn and myfunc be nested in the same file?"

That would be one approach, but I am not talking about sharing variables (like would be possible with nested functions). You need to get rid of the reliance of the event function on a value calculated by the ode function, so you need to be able to perform the same calculation in both routines: I am saying to create a common routine (could be in a different file) to call from both to do the common work.

With regards to using memoize:

   memoized_calc_dx2 = memoize(@calc_dx2);
   memoized_calc_dx2.CacheSize = 25;   %a guess
   memoized_calc_dx2.CachePolicy = 'LFU';  %release least frequently used
   memoized_calc_dx2.clearCache;

   ...

   options = ... 'event', @(t, x) event_function(t, x, xevent, memoized_calc_dx2) )
   [T, X] ode45(@(t, x) myfunc(t, x, xevent, memoized_calc_dx2), .... options)

   %did the cache help?

   memoized_calc_dx2.stats.Cache

   function dx = myfunc(t, x, xevent, calc_dx2)
      dx(1) = x(2);
      dx(2) = calc_dx2(x, xevent);

    function [value, isterminal, direction] = event_function(t, x, xevent, calc_dx2)
      ...
      dx2 = calc_dx2(x, xevent);
      ...

[Shant Danielian]

Hi Walter, I finally got a copy of R2017a installed and attempted to use your pseudo-code as a guide to update my own code. I attempted to run the code and got the following error message:

No public property CachePolicy exists for class matlab.lang.MemoizedFunction

I did a google search on it and nothing helpful came up. I also searched what you wrote in lines 2-4 to figure out what you're doing and nothing comes up as well. Could you please explain what they mean or point me in the right direction?

[Walter Roberson]

Change CachePolicy to CachingPolicy

See toolbox/matlab/lang/+matlab/+lang/MemoizedFunction.m near line 138

[Shant Danielian]

Thank you Walter, that fixed the issue and the code runs.

Answer 1

[Walter Roberson]

http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

parameterize both the ode function and the event function to pass in xevent

Do not make xevent persistent within the ode function: when you do that, the variable's identity as persistent overrides the value passed in as a parameter.

If xevent is something being calculated inside the ode function rather than something being passed in, then do not try to remember it for later use in the event function by using persistent or global or a shared variable: There are circumstances under which ode45 might not call the event function for a while, and there are circumstances under which ode45 might call the event function several times in a row with different parameters without calling the ode function between. Therefore any value that you calculate in the ode function that you would like to also use in the event function needs to be recalculated in the event function.

[Shant Danielian]

Thank you for your answer.

The reason I have xevent as persistent is because it's value needs to stay the same until the solver detects the next event, and then an new xevent will be passed into myfunc. It isn't something explicitly calculated within myfunc.

I'm aware that the events function is called multiple times and I too want to avoid using persistent variables within the event function. So would using a finite difference to get a rough approximation to dx(2) be worthless?

I will look into trying to parameterize the myfunc and myeventsfcn and see where I get, but the solution to my issue isn't super clear as of yet.

[Walter Roberson]

If xevent is changing during any one call to ode45() then you run the risk that you are introducing a discontinuity in the calculation or in one of of the next two derivatives of the calculation.

Discontinuities in the next one derivative will typically be detected and ode45 will narrow in on the time it happens at, trying to find a smooth transition and being unable to find it.

Discontinuities in the second derivative after what is directly calculated will not necessarily be detected by ode45 but will typically lead to errors in calculations.

Therefore if you have a parameter that is changing as it appears you are saying it might, then you should use an event function to detect the circumstances under which the transition should be made and terminate the ode45, and then call ode45 again with the rest of the timespan and with the boundary conditions set to the output of the previous call, but with the changed value of the parameter. Which you would pass in by parameterizing the call like I linked to.

[Shant Danielian]

That's exactly what I did. I stopped the solver when the discontinuity occurs, inputted a new xevent and restarted the solver at the updated initial conditions. I'm just a little confused on how to use the parameterization to input the xevent and the rest of the parameters that calculate dx(2) into the event function. The event I'm looking for is a function of dx(2) as well.