0

I have the following ode function (below), I am trying to adjust the function to allow the variable 'k_mb' to be dependent on the value of z(1) as seen in the function torsional_stiffness_changer, but also be equal to 5 when z(1) = 0.

The current state of the code doesn't compute as I'm not sure how to properly use the z(1) variable in this manner.

Any help is much appreciated, thanks!

    function zdot = odenesttest(t,z)
global r_p r_g I_p I_g m_p m_g V
persistent k_mb;
k_mb = 5;
A_or_B = [5,10];
V = [0,5,15,20,30,35,45,50,60];
%% state space equations are defined below:
zdot(2) = ((r_p*k_mb)/I_p)*(-r_p*z(1)+r_g*z(3)+z(5)-z(7));  % Eq 1ss
zdot(6) = (k_mb/m_p)*(r_p*z(1)-r_g*z(3)-z(5)+z(7));         % Eq 2ss
zdot(4) = ((r_g*k_mb)/I_g)*(-r_g*z(3)+r_p*z(1)+z(7)-z(5));  % Eq 3ss
zdot(8) = (k_mb/m_g)*(-r_p*z(1)+r_g*z(3)-z(7)+z(5));        % Eq 4ss

zdot(1) = z(2); % Eq 5ss
zdot(3) = z(4); % Eq 6ss
zdot(5) = z(6); % Eq 7ss
zdot(7) = z(8); % Eq 8ss

torsional_stiffness_changer
   function torsional_stiffness_changer
       if 0 <= z(1)
           lower_index = find(V < z(1), 1, 'last');
           even_or_odd = mod(lower_index,2);
           k_mb = A_or_B(2 - even_or_odd);
       end
   end
% Outputs
zdot = zdot';
end
Sirsh
  • 1
  • ODEs and ODE solvers by their very nature assume a certain degree of continuity. Inserting a discontinuity will result in a [stiff system](http://en.wikipedia.org/wiki/Stiff_equation) that can be difficult to integrate numerically and even lead to errors. Instead, you should simulate this as a piecewise continuous system and use [event detection](http://mathworks.com/help/matlab/math/ode-event-location.html) to switch between them. See [this](http://stackoverflow.com/q/13755352/2278029) or [this](http://stackoverflow.com/a/22819748/2278029). – horchler Jan 20 '17 at 05:24
  • Hi @horchler thanks for your reply. I've looked at the mathworks documentation, but cannot seem to wrap my head around how to implement what I'm trying to achieve. I've attempted it by using: [code] options = ('Events', @events'); [/code] Where events is: [code] function [k_mb,isterminal,direction] = events(t,z) lower_index = find(V < zdot(9), 1, 'last'); even_or_odd = mod(lower_index,2); k_mb = A_or_B(2 - even_or_odd); isterminal = 0; direction = 0; end [/code] However, I can't seem to make it work as k_mb needs [1] – Sirsh Jan 22 '17 at 02:34
  • to be defined before the ode can solve.. [2] – Sirsh Jan 22 '17 at 02:38

0 Answers0