ode 45 in a loop

I am trying to solve an ode in a loop, trying to get a single vx and tx value as ouput. I wish to use a loop because a parameter is to be changed everytime a new loop starts. I am getting length 0 output after calling the ode45 function. I wish to integrate the ODE over just one time interval that I have specifed (t(i:i+1)) to get a single value ouput in tx and vx Can you explain how to rectify the issue?

t = linspace(1,10);
vf(1) = 30 ; %some random constant value
vfx = vf(1); % needed in the function ode_research_x
x(1) = -3;% starting point for x
for i = 1:length(t)-1
    vfx = vf(1); % this vfx is then used by the ode_research_x function
    [tx,vx] = ode45('ode_research_x', t(i:i+1), vp(1)); % I am not sure but the time step seems like an issue
    x(i+1) =  x(i) + vx.*tx;
    vp(1) = vx;
    vf(1) = 30; % some new constant value calculated in the loop (I've shown a random value)
end

%function to define the ode
function dudt = ode_research_x(vx,vfx)
dudt = 18*mu_f*Cd*Rep*(vx - vfx)/(rho_p*d^2*24); % all the variables except vx and vfx are constant values in the same workspace so I haven't shown their values but they are defined 
end

Answer 1

[Bjorn Gustavsson]

You have a couple of problems with your code. Your ode_research does not follow the ordinary form (dydy = ode_f(t,y)) that ode45 expects. Then your use of a string for the first argument is a bit arcane. I would expect something like this would be better:

t = linspace(1,10);
vf(1) = 30 ; %some random constant value
vfx = vf(1); % needed in the function ode_research_x
x(1) = -3;% starting point for x
vp = v_initial; % or whatever you have as initial condition
for i = 1:length(t)-1
    vfx = vf(1); % this vfx is then used by the ode_research_x function
    [tx,vx] = ode45(@(t,vx) ode_research_x(t,vx,vfx), t(i:i+1), vp(1)); % I am not sure but the time step seems like an issue
    x(i+1) =  x(i) + vx(end).*tx;
    vp(1) = vx;
    vf(1) = 30; % some new constant value calculated in the loop (I've shown a random value)
end
%function to define the ode
function dudt = ode_research_x(t,vx,vfx)
% Some initialization of the constants, I assume.
dudt = 18*mu_f*Cd*Rep*(vx - vfx)/(rho_p*d^2*24); % all the variables except vx and vfx are constant values in the same workspace so I haven't shown their values but they are defined 
end

This looks like some attemt to integrate an equation of motion, if that is the case this is not a good way to go about it. In that case you should integrate the second-order ODE (i.e. the two coupled first order-ODEs) and handle the evolution of the force inside the ode_research_x - function.

HTH

[Sanjana Singh]

Thank you, the original problem got rectified but I don't understand your point about how to solve the Equation of Motion. I need to solve for vfx after every step in the ode because that gets updated after getting an updated x value. I am currently getting a really huge vector for vx (size is 2449 by 1 ) with zeros and NaN values. Can you explain in a bit more detail how I can solve this issue?

[Bjorn Gustavsson]

To integrate an equation of motion with ode45 you have to convert:

to

Then you integrate those coupled first-order ODEs from some initial condition [x0,v0] over your time-intervall of interest, since you seem to have a 1-D equation of motion this becomes 2 coupled ODEs. Something like this:

x0v0 = [x0;v0]; % your initial position and velocity as initial-condition-vector.
t_span = linspace(1,10); % Or whatever time-intervall you prefer.
[t,xv] = ode45(@(t,xv) ode_eq_motion_x(t,xv,other_pars),t_span,x0v0);

where ode_eq_motion_x now looks something like this:

function dxdxdvdt = ode_eq_motion_x(t,xv,other_pars)
% other_pars is just included as an example for how to feed the ODE with additional
% parameters,
dxdtdvdt = [xv(2); % That is just the second equation above
            F_over_m(t,xv,other_pars)];
    function F = F_over_m(t,xv,pars)
        mu_f = pars(1);
        Cd = pars(2);
        Rep = pars(3);
        d = pars(4);
        ... etc
        F = 18*mu_f*Cd*Rep*(vx - vfx)/(rho_p*d^2*24);
    end
end