How do I solve an ODE with time-dependent parameters in MATLAB?

Consider the following ODE with time-dependent parameters:

y'(t) + f(t)y(t) = g(t)

and the given initial condition:

y(0) = 1

This is an example of an ODE with time-dependent terms. Suppose the time-dependent terms are defined only through the set of data points given in two vectors. Which of the MATLAB ODE solvers should I use, and how do I set up this problem?

Answer 1

[MathWorks Support Team]

To solve this ODE, you must pass the data sets to the derivative function as additional parameters. When the ODE solver calls the derivative function, it will pass a specified time as the first input argument. You must then interpolate the datasets to obtain the value of the time-dependent terms at the specified time. This is performed within the following function, called myODE.

function dydt = myODE(t, y, ft, f, gt, g)
f = interp1(ft, f, t); % Interpolate the data set (ft, f) at times t
g = interp1(gt, g, t); % Interpolate the data set (gt, g) at times t
dydt = -f.*y + g; % Evalute ODE at times t

The function f is defined through the n-by-1 vectors tf and f, and the function g is defined through the m-by-1 vectors tg and g.

Now, you can refer to myODE within a call to a MATLAB ODE solver. Assume that the time-dependent parameters ft and gt are defined within the data sets generated by the following code.

ft = linspace(0, 5, 25); % Generate t for f 
f = ft.^2 - ft - 3; % Generate f(t)
gt = linspace(1, 6, 25); % Generate t for g
g = 3*sin(gt - 0.25); % Generate g(t)

The following code uses the ODE45 function to solve this time-dependent ODE.

TSPAN = [1 5]; % Solve from t=1 to t=5
IC = 1; % y(t=0) = 1
[T Y] = ode45(@(t,y) myODE(t, y, ft, f, gt, g), TSPAN, IC); % Solve ODE

Note that if you are using a version of MATLAB prior to MATLAB 7.0 (R14), you will need to pass the four additional parameters ft, f, gt, and g, into the ODE solver as follows.

TSPAN = [1 5]; % Solve from t=1 to t=5
IC = 1; % y(t=0) = 1
[T Y] = ode45(@myODE, TSPAN, IC, [], ft, f, gt, g); % Solve ODE

Now you can plot the solution y(t) as a function of time.

plot(T, Y);
title('Plot of y as a function of time');
xlabel('Time'); ylabel('Y(t)');

[Jan]

No, I do not agree. The ODE integrators are designed to operate on differentiable functions. The linear interpolation does not fullfill this requirement. In consequence the step-size control of the integrator cannot determine a proper value anymore and the results can be dominated by accumulated round-off error (in case of too small steps) or discretization errors (in case of too large steps). In both cases an analysis of the sensitivity (changes of the result to a variation of the initial positions or paramters) can yield to random results.

See https://www.mathworks.com/matlabcentral/answers/59582-ode45-errer-input-argument-y-is-undefined

I've seen many scientific projects which simply ignore the specifications of the integrators. But the calculated trajectories are not stable and this is not state of the art.

The only trustworthy solution is to run the integration in steps over the intervals with smooth parameters. Everything else is closing the eyes, cheating the step-size control and accepting a final value without considering the numerical tolerances.

[Steven Lord]

   function dydt = myODE(t, y, ft, f, gt, g)

   %snip most of the body of myODE

   [T Y] = ode45(@(t,y) myODE(t, y, ft, f, gt, g), TSPAN, IC); % Solve ODE

Don't call ode45 from within myODE to solve the system of equations whose right-hand side is defined by myODE itself. At best you receive the error you quoted. At worst you call myODE with six inputs which eventually results in ode45 calling myODE which calls ode45 which calls myODE which calls ode45 which calls myODE .... Eventually either MATLAB will complain about the recursion or (if you've set your recursion limit too high) MATLAB and/or your computer could crash!

Move the ode45 call (and whatever code is a prerequisite for that call, like the code that defines ft, f, gt, g, TSPAN, and IC) to a separate script or function or evaluate that call in the Command Window.

[Steven Lord]

Please show your updated scripts or functions.

As a guess, if you have a script file (your call to ode45) that contains a function (myODE) put the script content before the function rather than after.

[Steven Lord]

As I suspected, you put the function before the script content. Move the function to after the script content. But you need to do the interpolation inside myODE.

ft = linspace(0, 5, 25); % Generate t for f 
f = ft.^2 - ft - 3; % Generate f(t)
gt = linspace(1, 6, 25); % Generate t for g
g = 3*sin(gt - 0.25); % Generate g(t)

TSPAN = [1 5]; % Solve from t=1 to t=5
IC = 1; % y(t=0) = 1
[T Y] = ode45(@(t,y) myODE(t, y, ft, f, gt, g), TSPAN, IC); % Solve ODE

plot(T, Y);
title('Plot of y as a function of time');
xlabel('Time'); ylabel('Y(t)')

function dydt = myODE(t, y, ft, f, gt, g)
% You probably want to change one of the uses of the variables f and g
% to a different name
f = interp1(ft, f, t); % Interpolate the data set (ft, f) at times t
g = interp1(gt, g, t); % Interpolate the data set (gt, g) at times t

dydt = -f.*y + g; % Evalute ODE at times t
t = [1 5]; % this isn't necessary
end