Plot values changing with events in one continous graph

Hello,

I need help plotting a value, let's call it p. I have 3 different events set up and p should change depending on which event is active. I want a continous plot with only the value from the active event. I tried putting them inside the existing functions as nested functions, I tried codeing the values inside my main loop, I don't get a proper plot.

The functions should be as followed:

Events:

ZustandI: p=0

ZustandII: p=k2*(y(1)-y(3)-s)

ZustandIII: p=k2*(y(1)-y(3)+s)

Can you help me? My code is down below.

function mfc
clc
%%parameters
      m1=10; m2=10; m3=500; k3=300; k2=(k3*(m1+m2))/m3; my=m1/(m1+m2); FW=0; s=0.5; r=0.5; om=5; dv0s=1; dv0=dv0s-my*r*om;
%%starting parameters
    tspan = [0 100];
    tstart = tspan(1);
    tend = tspan(end);
    y0 = [0 dv0 0 0];

      % initiate starting function
      t = tstart;
      y = y0;
      fcn = @eqmi;
      opt = odeset('Events', @ZustandI);

      ateout = [];
      ayeout = [];
      aieout = [];
%%main loop
  while t(end) < tend
      % Integrate until event-function stops it
      [at,ay,ate,aye,aie] = ode45(fcn, [t(end) tend], y(end,:), opt);

      % new starting parameters
      t = cat(1, t, at(2:end));
      y = cat(1, y, ay(2:end,:));

      ateout = [ateout; ate];
      ayeout = [ayeout; aye];
      aieout = [aieout; aie];

      % new functions and events
      if (-s<y(end,1)-y(end,3) && y(end,1)-y(end,3)<s)
        fcn = @eqmi;
        opt = odeset('Events', @ZustandI);
      elseif y(end,1)-y(end,3)<=-s
        fcn = @eqmii;
        opt = odeset('Events', @ZustandII);
      elseif y(end,1)-y(end,3)>=s
        fcn = @eqmiii;
        opt = odeset('Events', @ZustandIII);
      end  
  end

figure(1);                           
plot(t,p)                            % this should detect events
hold on
plot(ateout,ayeout(:,1))
xlabel('time');
ylabel('test');
title('test');
legend('p');

%%functions
 function dy=eqmi(t,y)
     dy=zeros(5,1);
     dy(1)=y(2);  
     dy(2)=my*r*om^2*sin(om*t);                     
     dy(3)=y(4);        
     dy(4)=(FW*sin(w3*t)-k3*y(3))/m3;               
     dy(5)=0;                                        
  end

function dy=eqmii(t,y)
   dy=zeros(5,1);
   dy(1)=y(2);                                        
   dy(2)=my*r*om^2*sin(om*t)-(k2*(y(1)-y(3)-s))/(m1+m2); 
   dy(3)=y(4);                                             
   dy(4)=(FW*sin(w3*t)-k3*y(3)+k2*(y(1)-y(3)-s))/m3;   
   dy(5)=k2*(y(1)-y(3)-s);                               
end

    function dy=eqmiii(t,y)
        dy=zeros(5,1);
        dy(1)=y(2);                                            
        dy(2)=my*r*om^2*sin(om*t)-(k2*(y(1)-y(3)+s))/(m1+m2);   
        dy(3)=y(4);                                           
        dy(4)=(FW*sin(w3*t)-k3*y(3)+k2*(y(1)-y(3)+s))/m3;     
      dy(5)=k2*(y(1)-y(3)+s);                              
  end
    %%Events
    function [value,isterminal,direction] = ZustandI(t,y)
    value = [double(y(1)-y(3)>-s),double(y(1)-y(3)<s)];
    isterminal = [1,1];
    direction = [0,0];
    end
    function [value,isterminal,direction] = ZustandII(t,y)
    value = double(y(1)-y(3)<=-s);
    isterminal = 1;
    direction = 0;
    end
    function [value,isterminal,direction] = ZustandIII(t,y)
    value = double(y(1)-y(3)>=s);
    isterminal = 1;
    direction = 0;
    end
    end

[Jan]

I do not understand, what the question is exactly. "I want a continous plot with only the value from the active event" - this is not clear. In " I tried putting them", what is "them" and what does "nothing works" mean explicitly?

[Lennart]

Sorry, english is not my mother tongue but I try it as best as I can. I have set up 3 events with 3 corresponding functions. And I have a value p which I want to plot. p is equal to dy(5) from my 3 functions.

Now I want to plot p over the time t but as one graph. It should automatically switch to the corresponding function depending on which event is active. By the way my code works as intended, just the plotting of this one specific value is not.

Things I tried:

1. defining p as dy(5) inside my functions

2. defining p inside my if/else-loop

Those didnt work.

Answer 1

[Jan]

ateout = [];
ayeout = [];
aieout = [];
dy5    = [];

while t(end) < tend
  % Integrate until event-function stops it
  [at,ay,ate,aye,aie] = ode45(fcn, [t(end) tend], y(end,:), opt);

    % new starting parameters
    t = cat(1, t, at(2:end));
    y = cat(1, y, ay(2:end,:));

    dy  = fcn(at.', ay.');
    dy5 = cat(1, dy5, dy(:, 5));  % Or:  cat(1, dy5, dy(5, :).') ???

Now modify your functions to be integrated such, that they handle a vector input also:

function dy = eqmi(t, y)
dy       = zeros(5, numel(t));
dy(1, :) = y(2, :);  
dy(2, :) = my * r * om^2 * sin(om * t);                     
dy(3, :) = y(4, :);        
dy(4, :) = (FW * sin(w3 * t) - k3 * y(3, :)) / m3;               
dy(5, :) = 0;                                        
end

!UNTESTED! Care for transposing the inputs and/or outputs t and y.

[Lennart]

Hey, your help is greatly appreciated. First off, it doesn't affect the other calculations, so this is a good thing. Also the values for dy5 seem to be correct, I calculated them by hand and they coincide but only if I rename this:

dy5 = cat(1, dy5, dy(5,:));  % Or:  cat(1, dy5, dy(5,:).')

in something else like this:

p = cat(1, dy5, dy(5,:)); or cat(1, dy5, dy(5,:).')

Otherwise I get some weird starting entrys which aren't zero.

But how do I plot them correctly?

plot(t,p) or plot(t,p(:))

doesn't work because the vectors aren't the same length.

[Jan]

p = cat(1, dy5, dy(5,:)); or cat(1, dy5, dy(5,:).')

The second part of this is useless and should drop an error. I've added this comment only, because I'm not sure, if the variables have to be transposed or not.

t and p must have the same size after this code ran. Please post your code, if there are still problems.

[Lennart]

Yeah, I only included both parts here because I don't know if I should use

p = cat(1, dy5, dy(5,:)) or p=cat(1, dy5, dy(5,:).')

of course I only included one of this in my code. I tried it with both versions. I don't know why but the vectors of t and p aren't the same length.

My code is basically what I already posted in my initial questions but I added this like you suggested:

      ateout = [];
      ayeout = [];
      aieout = [];
      dy5    = [];

and

% new starting parameters
    t = cat(1, t, at(2:end));
    y = cat(1, y, ay(2:end,:));
    dy = fcn(at.', ay.');
    p=cat(1, dy5, dy(5,:).');

and I changed my functions to these:

% functions
function dy=eqmi(t,y)
    dy=zeros(5,numel(t));
    dy(1,:)=y(2,:);
    dy(2,:)=my*r*om^2*sin(om*t);
    dy(3,:)=y(4,:);
    dy(4,:)=(FW*sin(w3*t)-k3*y(3,:))/m3;
    dy(5,:)=0; 
end 
function dy=eqmii(t,y)
    dy=zeros(5,numel(t));
    dy(1,:)=y(2,:); 
    dy(2,:)=my*r*om^2*sin(om*t)-(k2*(y(1,:)-y(3,:)-s))/(m1+m2);
    dy(3,:)=y(4,:);
    dy(4,:)=(FW*sin(w3*t)-k3*y(3,:)+k2*(y(1,:)-y(3,:)-s))/m3;
    dy(5,:)=k2*(y(1,:)-y(3,:)-s);
end
function dy=eqmiii(t,y)
    dy=zeros(5,numel(t));
    dy(1,:)=y(2,:);
    dy(2,:)=my*r*om^2*sin(om*t)-(k2*(y(1,:)-y(3,:)+s))/(m1+m2);
    dy(3,:)=y(4,:);
    dy(4,:)=(FW*sin(w3*t)-k3*y(3,:)+k2*(y(1,:)-y(3,:)+s))/m3;
    dy(5,:)=k2*(y(1,:)-y(3,:)+s);
end

but my plot command plot(t,p) displays an error message because the vectors aren't the same length.

[Jan]

Ah, yes: If you have

t = cat(1, t, at(2:end));

you need crop the first elemet of dy also.

dy5 = cat(1, dy5, dy(5, 2:end).');

But

p = cat(1, dy5, dy(5,:).');

is not useful, because you do not accumulate dy5 anymore. Equivalently to concatenating t cumulatively with at, you collect the dy in dy5.

[Lennart]

Hey, ok I probably made a mistake earlier. This does seem to yield the correct values:

dy5 = cat(1, dy5, dy(5, 2:end).');

But it still hasn't the same length as my vector t:

t = cat(1, t, at(2:end));

But it is only a difference of one entrance. E.g. if I use the debugger and t is a 97x1 double vector at the same time dy5 is only a 96x1 double vector. I'm a little lost here on how to fix this.

[Lennart]

Ok, I found the problem. I just have to change

this dy5 = []; to this dy5 = 0;

Now it works just fine. Thank you for your help.