Plotting heaviside unit step functions

I'm struggling to plot Z(t) (a function with respect to t) from a differential equation in terms of heaviside step functions. My code is as below:

%This function is used to integrate using "ode45", (beta, gamma, A = constants) i.e find "Z(t)" from the differential equation

function [ dydt ] = IntegratingFunction2(t,y,beta,gamma,A)

z = y(1);
zdot = y(2);

dydt = zeros(size(y));
dydt(1) = zdot;
dydt(2) = heaviside(A*t) + 4*heaviside(A*(t-1)) - 5*heaviside(A*(t-3)) - beta*gamma*zdot - beta*z;

return

Then in my main script I called the function:

[tLinear,y] = ode45(@(t,y)IntegratingFunction(t,y,beta,gamma,A(k)), tspan, xinit,options);

I then made:

%This statement means that the first column from "y" are all the values for "z" (a function in terms of t(time)

z = y(:,1);

% the output "y" is technically dydt, as written in my function earlier % I now plot "z" in terms of "t" (time)

plot(tLinear , z , 'r--');

Now here is the problem, nothing shows/plots. Why?

Answer 1

[Richard Brown]

If you replace the line with the heavisides in it with the one from my previous answer

dydt(2) = (A*t >= 0) + 4*(A*(t-1) >= 0) - 5*(A*(t-3) >= 0) - beta*gamma*zdot - beta*z;

your code works perfectly. The reason you can't see the dashed lines is because they are all stacked on top of each other - if you put a pause statement inside your loop, you can see the solutions build up.

The reason for this is because for A > 0, Heaviside(A*(t-1)) is exactly equal to Heaviside(t - 1), so the different values of A have no effect. With the tanh function the sigmoid gets closer and closer to the Heaviside function as A increases.

[Sarah Coleman]

Thanks Richard, I did replace it but it still didn't work. The only difference is that I get a black dashed line with the other 5 curves (tanh). The black "dash" curve intersects at the top of the black "solid line" curve. Does this mean there still stacking on top?

[Richard Brown]

Yes, the fact that you can see the black dashed line means your simulation is working. All of your other dashed lines are underneath it (all of the Heaviside simulations give the same result)

[Sarah Coleman]

Richard, I really appreciate your help. I thought that the code would produce 5 different curves from the Heaviside simulations, but like you said it'll produce one btw good to see another Kiwi on this forum!

[Richard Brown]

No problem. Whereabouts are you from?

[Sarah Coleman]

Auckland, doing electrical engineering at Auckland Uni,

[Andy Zelenak]

Thanks Richard. This was useful to me, too. I didn't know you could put inequalities in an equation like that.

Answer 2

[Richard Brown]

heaviside is a function from the symbolic toolbox - I'd be surprised your code doesn't complain when you try to call it.

Why not instead use

dydt(2) = (A*t >= 0) + 4*(A*(t-1) >= 0) - 5*(A*(t-3) >= 0) - beta*gamma*zdot - beta*z;

And see how that goes ...

[Sarah Coleman]

Tried it, didn't make a difference.

[Richard Brown]

see my new answer below

Answer 3

[Sarah Coleman]

I've given my complete coding as below, which I've annotated for better understanding. I really have no clue what's happening and have tried various things which have all been fruitless.

function Plotting_Of_Different_A_Values

tstart = 0;
tend = 10;
n = 10000;
tspan = linspace(tstart,tend,n);

%Initial conditions for differential equation i.e. Z(0) = 0, dZ/dt(0) = 0.
%NB Z and dZdt are with respect to t(time)
xinit = [0;0];

%constants
beta = 55;
gamma = 0.2;

%Different values of a co-efficient
A = [0.5,1,1.5,2,5];

options = odeset('RelTol', 1e-10, 'AbsTol', 1e-10);

%For loop is for plugging in different values of A into each of the
%functions and finding out how that affects Z(t) which is determined using
%ode45

for k = 1:length(A)

    [t,x] = ode45(@(t,x)IntegratingFunction(t,x,beta,gamma,A(k)), tspan, xinit,options);
    z = x(:,1);

    %This function is used to find Z(t) with a different d^2Z/dt^2 than
    %IntegratingFunction, look at function for details
    [t2,y] = ode45(@(t,y)IntegratingFunction2(t,y,beta,gamma,A(k)), tspan, xinit,options);
    z2 = y(:,1);

   %plotting Z for different A values. 
    if A(k) == 0.5
        figure(1);
        plot(t,z,'r-',t2,z2,'r--'); 
        text(2,0.0447, 'A=0.5');
        text(2,0.21, 'A=0.5');
    elseif A(k) == 1
        plot(t,z,'g-',t2,z2,'g--'); 
        text(2,0.07, 'A=1');
        text(2,0.345, 'A=1');
    elseif A(k) == 1.5      
        plot(t,z,'b-',t2,z2,'b--');
        text(1.5,0.07, 'A=1.5')
        text(2,0.4, 'A=1.5')
    elseif A(k) == 2
        plot(t,z,'m-',t2,z2,'m--');
        text(2.75,0.075, 'A=2')
        text(2,0.425, 'A=2')
    else    
        plot(t,z,'k-',t2,z2,'k--');
        text(2.65,0.4125, 'A=5')
        text(1.15,0.075, 'A=5')
    end    

    hold on
    xlabel('t(s)');
    ylabel('Z');

end
end

This is the problem I'm facing, I can plot all the different curves for %z = x(:,1) but cannot get any curves for z2. Why? It's almost as if %IntegratingFunction2 doesn't exist, maybe I'm mixing the variables between the 2 functions?

Below are the 2 Integrating functions:

IntegratingFunction:

function [ dxdt ] = IntegratingFunction( t,x,beta,gamma,A )

z = x(1);
zdot = x(2);

dxdt = zeros(size(x));
dxdt(1) = zdot;
dxdt(2) = 0.5*((tanh(A*t)) - tanh(A*(t-1)) + 5*tanh(A*(t-1)) - 5*tanh(A*(t-3))) - beta*gamma*zdot - beta*z;

return

IntegratingFunction2:

[Sarah Coleman]

Code for IntegratingFunction2:

function [ dydt ] = IntegratingFunction2(t,y,beta,gamma,A)

z = y(1);

zdot = y(2);

dydt = zeros(size(y));

dydt(1) = zdot;

dydt(2) = heaviside(A*t) + 4*heaviside(A*(t-1)) - 5*heaviside(A*(t-3)) - beta*gamma*zdot - beta*z;

return