Solve 2nd order ODE with discrete time terms

I have a second order diferential equation:

k1 * d2a + k2 * da + a = e

where a and e are functions of t. I have e(t) in a matrix DATA where:

t = DATA(:,1)

e = DATA(:,2)

If I define the function to be used in ode45 solver as:

function dx = myFUN(t,x)
   dx = zeros(2,1);
   dx(1) = x(2);
   dx(2) = 1/k1*(-k2*x(2)-x(1)+ e(t));
end

How can I pass the value of e(t) on each time step?

Answer 1

[Paul]

Assuming some values for the constants, t1 is your t = DATA(:,1) and e1 is your e = DATA(:,2):

function dx = myFUN(t,x)
    t1 = 0:1/1000:1 ;
    e1 = 0:1000;  
    e_int=interp1(t1,e1,t);
    dx = zeros(2,1);
    k1=1;k2=2;
    dx(1) = x(2);
    dx(2) = 1/k1*(-k2*x(2)-x(1)+ e_int);
end

Because t during solving probably won't exactly be one of the t values in your data, you have to interpolate to estimate e at t (e_int).

[Fran]

Thanks Paul, I realized on this, but I didnt know how to do that. This solution works fine and that is what I needed, but just for curiosity, could it be possible to pass e(t) at each time step instead of to compute t1 and e1 on each iteration?

Thanks again for your answer

[Paul]

Yes, there are several solutions for this. You could parse t1 and e1 as extra arguments of the function file:

function dx = myFUN(t,x,t1,e1)
    e_int=interp1(t1,e1,t);
    dx = zeros(2,1);
    k1=1;k2=2;
    dx(1) = x(2);
    dx(2) = 1/k1*(-k2*x(2)-x(1)+ e_int);
end

t1 = 0:1/1000:1 ;
e1 = 0:1000;
[t,y]=ode45(@myFUN,[0;1],[0 0],[],t1,e1)

The [] are needed because that's normally where the options of ode45 (see odeset) would normally be. If you use the default options you put [] there to indicate this.

You could also for example use global variables, however this is frowned upon:

function dx = myFUN(t,x)
    global t1 e1
    e_int=interp1(t1,e1,t);
    dx = zeros(2,1);
    k1=1;k2=2;
    dx(1) = x(2);
    dx(2) = 1/k1*(-k2*x(2)-x(1)+ e_int);
end

global t1 e1
t1 = 0:1/1000:1 ;
e1 = 0:1000;
[t,y]=ode45(@myFUN,[0;1],[0 0])