Problem using OutputFcn to get parameters out of ODE45
8 次查看(过去 30 天)
显示 更早的评论
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 个评论
Jan
2021-1-7
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?
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!