Usage of ode function
19 views (last 30 days)
Show older comments
So I need to find the values for dv/dt and the dr/dt, the first derivate of the functions. The problem is every else variable is a vector in a form nx1,
and i have different values for every variable depending of time. How do I use properly the ode function to get the solution?
0 Comments
Answers (2)
Jon
on 23 Jan 2024
Edited: Jon
on 23 Jan 2024
Here is one approach, see https://www.mathworks.com/help/optim/ug/passing-extra-parameters.html for further details on passing extra parameters.
This example is just intended to give you an idea of how to approach this. You will have to modify it based on your needs, and actual values
% Make up some values just to demonstrate program
m = 1;
I = 2;
a = 8;
b = 5;
u = 4;
nPts = 100;
tFinal = 23;
ty = linspace(0,tFinal,nPts);
Fyf = randn(nPts,1); % time varying force Fyf
Fyr = randn(nPts,1); % time varying force Fyr
v0 = 0; % initial velocity
r0 = 0; % initial position
[t,v,r] = mySolve(m,I,a,b,u,ty,Fyf,Fyr,v0,r0);
plot(t,v,t,r);
legend('v','r')
ylabel('position and velocity')
xlabel('time')
% Everything below this line could be put into a separate file
function [t,v,r] = mySolve(m,I,a,b,u,ty,Fyf,Fyr,v0,r0)
%motionSolve - solves my system of ode's
% m,I,a,b,u - System parameters (assumed to be scalar constants)
% ty - vector of time values where forces are defined
% Fyf,Fyr - vectors of force values occuring at times corresponding to
% elements of ty
% Find start and end time
tspan = [min(ty),max(ty)];
% Set initial condition
x0 = [v0,r0];
% Solve ode
[t,x] = ode45(@odefun,tspan,x0);
% Return desired values
v = x(:,1);
r = x(:,2);
% Use nested function so that external inputs, Fyf and Fyr and parameters
% are in scope when computing derivative values
function dxdt = odefun(t,x)
% Calculates derivatives with respect to time
% Interpolate to find force values at current time
ff = interp1(ty,Fyf,t);
fr = interp1(ty,Fyr,t);
% Assign local variable names for readability
v = x(1);
r = x(2);
% Calculate derivatives
dxdt(1,1) = 1/m*(ff + fr ) - u*r;
dxdt(2,1) = 1/I*(a*ff - b*fr);
end
end
1 Comment
Jon
on 23 Jan 2024
@Torsten exploits the special structure of your ode's in which dr/dt is only a function of the external inputs, so you can just integrate to obtain r(t). Since dv/dt depends only on r(t), and the external inputs (but not upon v(t)) you can then just integrate to find v(t).
The approach I have is more general, and would work for ode's that can not be directly integrated and require a solver such as ode45
Torsten
on 23 Jan 2024
Edited: Torsten
on 23 Jan 2024
I assume you know r0=r(t=0) and v0=v(t=0).
Then
r = r0 + cumtrapz(t,(a.*F_yf - b.*F_yr)./I)
v = v0 + cumtrapz(t,(F_yf+F_yr-m.*u.*r)./m)
where t is the nx1 time vector for which the other time-dependent variables (F_yf,...) are saved in nx1 vectors.
1 Comment
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!