Solve 2nd order ODE with discrete time terms

5 views (last 30 days)
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?

Accepted Answer

Paul
Paul on 24 Feb 2014
Edited: Paul on 24 Feb 2014
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).
  2 Comments
Fran
Fran on 25 Feb 2014
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
Paul on 25 Feb 2014
Edited: Paul on 25 Feb 2014
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])

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!