Getting a Local Variable out of an ODE function

I have the following code to perform an ODE integration:

PLOT_ELEMENT=[1];

Mg=100;

options = odeset('RelTol',INT_TOL,'AbsTol',INT_TOL*ones(1,MATRIX_SIZE(1)*2),'OutputSel',PLOT_ELEMENT,'OutputFcn',@odeplot);  

ODEFCN= @(time,X) mass_spring_damp(time,X,Mg);

[time,X] = ode15s(ODEFCN,[TMIN TMAX],X_DOT_INIT,options);

Inside the mass_spring_damp.m file, I have a number of calculations that define local variables and ultimately the XDOT equation. the integration is working properly, but I want to store one of the local variables to a workspace at each time step for future plotting. How do I do that? I've been trying to make my own outputfcn, but have not had any luck transferring the local variable from my ODE function to the outputfcn.

I have also been able to set up a for loop after the fact and call my mass_spring_damp function for each time step and X, but that is very slow to process. I'd rather write something each step as I'm going through the integration.

Answer 1

[Walter Roberson]

You would use the same basic techniques as described in http://matlab.wikia.com/wiki/FAQ#How_can_I_share_data_between_callback_functions_in_my_GUI.28s.29.3F

In your situation my preference would be to create mass_spring_damp as a nested function with a shared variable that mass_spring_damp added to.

Caution:

The t, y outputs of ode15s are not necessarily the same as the t, y inputs. The ode* routines probe near values and project what the output values would be at the times they wish to report on. The ode* routines are, generally speaking, permitted to go backwards and forwards in time, such as to locate potential discontinuities or as needed to meet integration tolerances. My experiments with ode45() suggest that it might possibly not change direction in time but that when length(y)>1 that it will definitely use the same t multiple times, to probe the behavior of changing the other values.

Because of this, to make sense of your intermediate values, you really need to record the t and y inputs along with the intermediate values.

Furthermore, if your expectation was that recording the intermediate values would save you from having to recompute them with the final t, y output, then because the t, y inputs are not necessarily in the same location as the t, y outputs, recording the intermediate values might not help.

[Doug Scott]

Thank you, Walter. Do you know of a way to take advantage of the "flag" variable that is used in the outputfcn call? I know that one tests to see if the integrator is done before writing data. If I could get access to that in my mass_spring_damp function, I could test to see if it's done before writing using the assignin or similar function.

[Walter Roberson]

https://www.mathworks.com/help/matlab/ref/odeset.html#input_argument_namevalue_OutputFcn

"The solver calls myOutputFcn([],[],'done') once integration is complete to allow the output function to perform cleanup tasks."

Notice that no data is passed, just the flag value.

If you are holding data inside your ode function in a persistent variable then you could call your function with some parameters that signal it is time to do the saving.

Note: if the data you are holding on to in the ode function is used as part of the computations of the ode, then you are at risk that the ode function is discontinuous, since a different path to the same boundary condition might not give the same result due to the history of computation that you are saving.

Answer 2

[Doug Scott]

Would it be possible to test for flag==[] to determine when the integrator is writing the data for the current converged time step and somehow gain access to the local variables in the mass_spring_damp ODE function?

[Walter Roberson]

You can set an outputFcn and in that function you can test isempty() on the flag. However, the outputFcn is called after your ode function has finished executing for the particular set of conditions it was given, so your ode function has already exited. Any data that the ode function did not save somewhere else (using one of the techniques I linked to) is gone, with the exception of any data the ode function saved in a persistent variable. The only way to access that persistent variable after the ode function has finished executing is to call the ode function again, passing in something that it can tell means you want it to do something with the values it saved.

At the time the ode function is executed normally, the flag does not yet exist. The ode15s() routine has not yet decided it is going to terminate -- it does not decide that until it sees the result of the ode function call and decides that the results of the integration are within tolerances (or, alternately, that the results indicate that a singularity has probably been encountered.)