Problem using OutputFcn to get parameters out of ODE45
8 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Hi,
I'm using ode45 to solve a three degree of freedom mass, spring, damper problem. It's been running well until I've modified the code to extract parameters calculated during the solving process using odeset with OutputFcn. I've used the following previous questions for help:
I am now receiving an error saying that an index exceeds the number of array elements (0) and the line in question is in myOutputFcn. It is presumably x that is returning the error:
tanadr = tan(deg2rad(adr_dbz.*(x(1)+lr.*x(3))+adr));
Where my code may differ from the examples I've used is that the calculations in myOutputFcn use the current values of the displacement and velocity variables that are being solved for in the ode. x(1) is the body vertical displacement and x(3) is the body pitch displacement.
The solver setup and calling code is:
% Set up ODE function to pass to solver
odeFcn = @(t, x) Three_DOF_Half_Car_ODE(t, x, ix, iz, ip, tvec, mb, Iby, lf, lr, h, ...
kfz, cfz, kfx, cfx, adf, adf_dbz, krz, crz, krx, crx, adr, adr_dbz, Fxr_max);
% Set up options for output function
options = odeset('OutputFcn',@(t, x, flag) myOutputFcn(t, x, flag, ix, iz, ip, tvec, mb, Iby, lf, lr, h, ...
kfz, cfz, kfx, cfx, adf, adf_dbz, krz, crz, krx, crx, adr, adr_dbz, Fxr_max));
% Run solver
[T,x] = ode45(odeFcn, tvec, x0, options);
The ode contains, for example:
function dx = Three_DOF_Half_Car_ODE(t, x, ix, iz, ip, tvec, mb, Iby, lf, lr, h, ...
kfz, cfz, kfx, cfx, adf, adf_dbz, krz, crz, krx, crx, adr, adr_dbz, Fxr_max)
ix = interp1(tvec, ix, t); % this is interpolating ix at t within tvec
iz = interp1(tvec, iz, t);
ip = interp1(tvec, ip, t);
Fxr_max = interp1(tvec, Fxr_max, t);
% Calculation of anti-dive tangent wrt body displacement
tanadf = tan(deg2rad(adf_dbz.*(x(1)-lf.*x(3))+adf));
tanadr = tan(deg2rad(adr_dbz.*(x(1)+lr.*x(3))+adr));
% Contact patch longitudinal displacement and velocity equations
xif = -(x(1)-lf.*x(3)).*tanadf;
xir = (x(1)+lr.*x(3)).*tanadr;
xdotif = -(x(4)-lf.*x(6)).*tanadf;
xdotir = (x(4)+lr.*x(6)).*tanadr;
% Calculate rear axle longitudinal force and limit it according to decel
% Only valid during braking, when fxr is negative
Fxr = krx.*(-xir-x(2))+crx.*(-xdotir-x(5));
if Fxr < Fxr_max
Fxr = Fxr_max;
end
dx(1) = x(4);
dx(2) = x(5);
dx(3) = x(6);
dx(4) = 1./mb.*(kfz.*(-x(1)+lf.*x(3))+cfz.*(-x(4)+lf.*x(6))+krz.*(-x(1)-lr.*x(3))+crz.*(-x(4)-lr.*x(6))-...
(kfx.*(-xif-x(2))+cfx.*(-xdotif-x(5))).*tanadf+Fxr.*tanadr+mb.*iz);
dx(5) = 1./mb.*(kfx.*(-xif-x(2))+cfx.*(-xdotif-x(5))+Fxr+mb.*ix);
dx(6) = 1./Iby.*(-kfz.*lf.*(-x(1)+lf.*x(3))-cfz.*lf.*(-x(4)+lf.*x(6))+krz.*lr.*(-x(1)-lr.*x(3))+crz.*lr.*(-x(4)-lr.*x(6))-...
kfx.*h.*(-xif-x(2))-cfx.*h.*(-xdotif-x(5))-Fxr.*h+...
lf.*(kfx.*(-xif-x(2))+cfx.*(-xdotif-x(5))).*tanadf-lr.*Fxr.*tanadr+Iby.*ip);
It is the variable Fxr that I am trying to extract, so the following is in myOutputFcn (including the line beginning tanadr at which the code fails):
function status = myOutputFcn(t, x, flag, ix, iz, ip, tvec, mb, Iby, lf, lr, h, ...
kfz, cfz, kfx, cfx, adf, adf_dbz, krz, crz, krx, crx, adr, adr_dbz, Fxr_max)
persistent Fxr
Fxr_max = interp1(tvec, Fxr_max, t);
tanadr = tan(deg2rad(adr_dbz.*(x(1)+lr.*x(3))+adr));
xir = (x(1)+lr.*x(3)).*tanadr;
xdotir = (x(4)+lr.*x(6)).*tanadr;
Fxr_tmp = krx.*(-xir-x(2))+crx.*(-xdotir-x(5));
if Fxr_tmp < Fxr_max
Fxr_tmp = Fxr_max;
end
switch flag
case 'init'
Fxr = Fxr_tmp;
case ''
Fxr = [Fxr, Fxr_tmp];
case 'done'
assignin('base','Fxr',Fxr);
end
status = 0;
I'd be grateful for any guidance.
Regards, Simon.
2 comentarios
Jan
el 7 de En. de 2021
By the way:
if Fxr < Fxr_max
Fxr = Fxr_max;
end
This looks like your function to be integrated is not smooth. Then the computed trajectory can be dominated by rounding errors or the integration can stop with an error. Using ODE45 for non-smooth functions replies a trajectory, but this is not a reliable result from a scientifiv point of view.
Use an event function to stop and restart the integration instead.
The same problem concerns the linear interpolation by interp1.
It is still not clear to me, what you want to achieve. You calculate the trajectory already. As far as I can see, you try export one component a second time, but why?
Respuestas (0)
Ver también
Categorías
Más información sobre Ordinary Differential Equations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!