Plotting heaviside unit step functions

6 views (last 30 days)
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?

Accepted Answer

Richard Brown
Richard Brown on 19 Apr 2012
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.
  6 Comments
Sarah  Coleman
Sarah Coleman on 19 Apr 2012
Auckland, doing electrical engineering at Auckland Uni,
Andy Zelenak
Andy Zelenak on 7 Dec 2014
Edited: Andy Zelenak on 7 Dec 2014
Thanks Richard. This was useful to me, too. I didn't know you could put inequalities in an equation like that.

Sign in to comment.

More Answers (2)

Richard Brown
Richard Brown on 19 Apr 2012
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
Sarah Coleman on 19 Apr 2012
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:
  1 Comment
Sarah  Coleman
Sarah Coleman on 19 Apr 2012
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

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!