# How to simulate a system output with a sine wave input?

I wish to simulate the output of a certain gear system I have. How the gear system looks isn't particularly important to the problem, I managed to get the differential equation needed from the mechanical system. Here is the code I have (big thanks to the person who helped me with this)

% parameters
N2  = 90;
N1  = 36;
Jn1 = 0.5;
Jn2 = 0.8;
J2  = 2;
D   = 8;
K   = 5;
J   = (N2/N1)^2 * Jn1 + Jn2 + J2;

% define the system
sys = ss([0 1; -K/J -D/J], [0; N2/(N1*J)], [1 0], 0);

x0 = [0; 0];

% define the time span
t = linspace(0, 15, 10000)';

% define the input step
T1 = zeros(length(t), 1);
T1(t>=0) = 1;

% compute the system step response at once
theta1 = lsim(sys, T1, t, x0);

% compute the system response as aggregate of the forced and unforced
% temporal evolutions
theta2 = lsim(sys, T1, t, [0; 0]) + initial(sys, x0, t);

% plot results
figure('color', 'white');
hold on;
yyaxis left;
plot(t, T1, '-.', 'linewidth', 2);
ylabel('[N]');
yyaxis right;
plot(t, theta1, 'linewidth', 3);
plot(t, theta2, 'k--');
xlabel('t [s]');
grid minor;
legend({'$$T_1$$', '$$\theta_1$$', '$$\theta_2$$'}, 'Interpreter', 'latex',...
'location', 'southeast');
hold off;


This should work in generating a graph that shows the positions, my outputs, for a Heaviside/step input. My question is, how would I do this for a sine wave input. I figure I should have sin(w*t) instead of (t>=0), where w is my pulse frequency. Still, I can't seem to make this work. Any help would be really appreciated! :)

• Could you elaborate on what you have tried for the sine wave? Did you use T1=sin(w*t)? Feb 3 at 5:29
• Hello there, yes I did and it didn't work, an error appears. What I suspect the problem may be are the negative values in the arrays which are generated by the sine wave.
– Marx
Feb 3 at 6:30
• Could you show in more detail what you tried and what errors you got in return? Feb 3 at 6:48

Marx, you already asked me this and I've provided the answer to that in the bottom line of https://robotics.stackexchange.com/a/21740/6941.

Let me recall the function here as well:

function x = integrate(t, u, x0)
% parameters
N2  = 90;
N1  = 36;
Jn1 = 0.5;
Jn2 = 0.8;
J2  = 2;
D   = 8;
K   = 5;
J   = (N2/N1)^2 * Jn1 + Jn2 + J2;

% integrate the differential equation
[t, x] = ode23(@fun, t, x0);

% plot results
figure('color', 'white');

% plot position
yyaxis left;
plot(t, x(:, 1));
ylabel('$$x$$ [rad]', 'Interpreter', 'latex');

% plot velocity
yyaxis right;
plot(t, x(:, 2));
ylabel('$$\dot{x}$$ [rad/s]', 'Interpreter', 'latex');

grid minor;
xlabel('$$t$$ [s]', 'Interpreter', 'latex');

function g = fun(t, x)
g = zeros(2, 1);
g(1) = x(2);
g(2) = (-K/J)*x(1) + (-D/J)*x(2) + (N2/(N1*J)*u(t));
end
end


Example call:

t = linspace(0, 15, 10000)';
x0 = [0.1; 0];
x = integrate(t, @(t)(sin(0.5*t)), x0);

• I've just seen your answer, thank you so very much! :)
– Marx
Feb 3 at 13:47
• For linear time invariant systems lsim() might be more accurate than ode23() (depending on the time step size for lsim and the tolerances for ode23). Feb 3 at 16:30
• I know, but the OP specifically asked for ode23 in the other Q that I pointed at. Feb 3 at 17:23