Plotting heaviside unit step functions

9 次查看（过去 30 天）
Sarah Coleman
Sarah Coleman 2012-4-19
编辑： Andy Zelenak 2014-12-7
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?

请先登录，再回答此问题。

采纳的回答

Richard Brown
Richard Brown 2012-4-19
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 个评论
显示 4更早的评论
Sarah Coleman
Sarah Coleman 2012-4-19
Auckland, doing electrical engineering at Auckland Uni,
Andy Zelenak
Andy Zelenak 2014-12-7
编辑：Andy Zelenak 2014-12-7
Thanks Richard. This was useful to me, too. I didn't know you could put inequalities in an equation like that.

请先登录，再进行评论。

更多回答（2 个）

Richard Brown
Richard Brown 2012-4-19
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 2012-4-19
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 个评论
Sarah Coleman
Sarah Coleman 2012-4-19
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 MATLAB 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by