ode 45 in a loop
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Sanjana Singh
am 4 Jun. 2020
Kommentiert: Bjorn Gustavsson
am 4 Jun. 2020
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
0 Kommentare
Akzeptierte Antwort
Bjorn Gustavsson
am 4 Jun. 2020
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
2 Kommentare
Bjorn Gustavsson
am 4 Jun. 2020
To integrate an equation of motion with ode45 you have to convert:
![d2x(t)/dt2 = F(t,x)](https://www.mathworks.com/matlabcentral/answers/uploaded_files/309845/dt2%20=%20F(t,x).png)
to
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/309848/image.png)
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
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Loops and Conditional Statements finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!