Solve 2nd order ODE with discrete time terms

Question

Fran 2014-2-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/118875-solve-2nd-order-ode-with-discrete-time-terms

编辑： Paul 2014-2-25

采纳的回答： Paul

在 MATLAB Online 中打开

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?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Paul 2014-2-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/118875-solve-2nd-order-ode-with-discrete-time-terms#answer_125953

编辑：Paul 2014-2-24

在 MATLAB Online 中打开

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 个评论
显示无隐藏无

Fran 2014-2-25

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 2014-2-25

编辑：Paul 2014-2-25

在 MATLAB Online 中打开

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])