Solve ODE with array forcing function

I want to solve a simple ODE of the form like below ( current in a RL circuit):

dIdt= k1*I+k2*V, where k1 and k2 are constants.

V is a vector of recorded voltage data that can be thought of as V(t) at a specified sampling frequency fs. All examples I find have the forcing term with a dependency on t like sin(t).

How can I formulate this so that my forcing function V gets used without an explicit dependency on t?

Something like:

[t,I]=ode45( @rhs, t, initial_I, k1,k2 V);

function dIdt=rhs(t,x)

dIdt = k1*I + k2*V;

end

Answer 1

[Swastik Sarkar]

Hi @MR,

The ode45 function is typically used with continuous functions. When working with custom discrete functions, interpolation is necessary so that the solver can handle the data. Here is an example demonstrating this approach:

k1 = -0.5; 
k2 = 1.0; 

V = [0, 1, 0.5, -0.5, -1, 0, 0.8, 0.2, -0.2, -0.8, 0];
T = numel(V) - 1;
t_V = 0:T;

V_interp = @(t) interp1(t_V, V, t, 'linear', 'extrap');

rhs = @(t, I) k1 * I + k2 * V_interp(t);

t_span = [0, T]; % Going beyond this value may bring in undesired results

[t, I] = ode45(rhs, t_span, 0);

figure;
plot(t, I, 'b-', 'LineWidth', 2);
xlabel('Time (s)');
ylabel('Current (I)');
title('Current in RL Circuit');
grid on;

figure;
stairs(t_V, V, 'r-', 'LineWidth', 1.5);
xlabel('Time (s)');
ylabel('Voltage (V)');
title('Voltage Input');
grid on

To learn more about ode45, please refer to the following documentation:

https://www.mathworks.com/help/matlab/ref/ode45.html

Hope this is helpful.

[MR]

Well then is there a built-in way to ''numerically'' solve these type of equations? By this I mean keeping everything discrete, no interpolation? I don't ''absolutely'' want to use ODE45.

Thanks.

[Torsten]

Well then is there a built-in way to ''numerically'' solve these type of equations?

This is the way to solve "these type of equations".

Look at the example

"ODE with Time-Dependent Terms"

under

https://uk.mathworks.com/help/matlab/ref/ode45.html

If you want to keep everything discrete, you will have to code "Euler forward" with the step size dt chosen as the difference between subsequent times at which V(t) is recorded.

V = ...;
T = ...;
I(1) = I0;
for i = 1:numel(T)-1
    I(i+1) = I(i) + (T(i+1)-T(i))*(k1*I(i)+k2*V(i));
end