How to solve an ODE with parameters calculated in another function?

Hi

I want to solve an ODE with parameters calculated in another function.

For example, I have the following ODE:

dy/dx = -5*y + f

where f is obtained from another function.

Then, how can I import this f into the ode solver?

I would really appreciate if anyone can help me out.

Thank you.

[James Tursa]

What does this other function depend on? Can you call it from within the derivative function?

[Hancheol Cho]

f is obtained in another function with a very complicated procedure. f has its value in each time step. For example, f is the control force obtained using a PID controller. I want to see the trajectory to check if the control force works well.

[James Tursa]

So f depends on time? I will ask again, can you simply call it from within the derivative function since you have time available there?

[Hancheol Cho]

I'm sorry, but I don't exactly understand your question. f depends on time, but cannot express as a function of time. I have only the result for each time step. For example, at t=0, f=0, at t=0.1, f=0.2, at t=0,2, f=0.45, and so on.

[James Tursa]

"... f depends on time, but cannot express as a function of time. ..."

This statement doesn't make sense to me. How can you possibly use f if you don't know how to calculate it as a function of time, given that it depends on time? I.e., what do you really have to work with here? Do you have a vector of preset values? Or what? I don't understand how you are getting your f values calculated.

[Steven Lord]

James, I interpreted that description as f being a vector of data, but the poster not knowing a functional relationship f = someFunction(t). It's like having:

 f = [0 1 4 9 16 25];
 t = [0 1 2 3 4 5];

without knowing/recognizing/being able to take advantage of the fact that f is just t.^2.

Answer 1

[Jan]

It looks like your f is not continuous, but Matlab's ODE solvers fail for non smooth functions. See http://www.mathworks.com/matlabcentral/answers/59582#answer_72047

So the simple solution seems to run the integration in steps:

t  = [0, 0.1, 0.2];
y0 = ???;
for k = 2:length(t)
  fValue = f(t(k-1));
  [T, Y] = ode45(@(t, y)@yourFcn(t, y, fValue), t(k-1:k), y0);
  y0     = Y(end, :);
end

Then add some code to collect the T,Y in an array.

Answer 2

[Steven Lord]

Consider following the approach given in the third example on the documentation page for ODE45. That example passes two vectors, one containing time and the other values, for each of the two functions f(t) and g(t) and uses INTERP1 to interpolate the data to obtain a value at the time at which the ODE solver is evaluating the ODE function.

[Jan]

This is a good example of the problem, which occur when the function to be integrated is not smooth. What a pitty that it appears as example in the documentation.

When the tolerance is reduced to 1e-8, ODE45 rejects 58 steps at the locations, where INTERP1 creates non-smooth values:

 >> Result.stats
     nsteps: 203
    nfailed: 58
    nfevals: 1567

When the the interpolation is performed in 'pchip' mode, the integration runs with less rejected steps:

 >> Result.stats
      nsteps: 147
    nfailed: 4
    nfevals: 907

This seems to be not dramatic in this case. But when the trajectory is not stable, the differences can matter. An analysis of the sensitivity would suffer (a measurement of the variation of the results, when the initial values and/or parameters are varied).

So consider the specifications of the ODE integerators and do not integrate non-smooth functions for scientific purpusose.

Answer 3

[Marc]

I would try a brute force method having your ODE function calling the function for f(....) right before you set up your dy/dx equations.