# 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)? Commented Feb 3, 2021 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. Commented Feb 3, 2021 at 6:30
• Could you show in more detail what you tried and what errors you got in return? Commented Feb 3, 2021 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! :) Commented Feb 3, 2021 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). Commented Feb 3, 2021 at 16:30
• I know, but the OP specifically asked for ode23 in the other Q that I pointed at. Commented Feb 3, 2021 at 17:23