% SOLVEEXAMPLE--Example of mechanical system solved using the
% symbolic toolbox of Matlab.
echo on % to display the progress of the program in the command window
% Define all variables and parameters that will appear in equations
syms Z1 Z2 F z2_0
syms M k s
% The equations will be written in a variable 'e'. This program is most
% easily modified if any previous version of 'e' is deleted from Matlab's
% workspace.
clear('e'); % removes 'e' from workspace if previously defined
% These are the equations of the system in the Laplace domain:
e{1} = F == k*(Z1-Z2);
e{2} = k*(Z1-Z2) == M*s^2*Z2 - M*s*z2_0 + k*Z2;
u = {Z1, Z2}; % these are the unknowns
% It is assumed that z2_0, denoting z2(0), and F are known.
% The constants k and M are also assumed to be known.
% The following solves the system of equations for Z1 and Z2
x = [e, u];
a = solve(x{:});
% This instruction prints the calculated Z1
pretty(a.Z1)
% Another way to print Z1
a.Z1
% If f = f0, we can substitute F by f0/s:
syms f0
expr = subs(a.Z1, F, f0/s);
expr = simplify(expr) % rewrite expr to a form easier to read
% Let's tell Matlab that k and M are positive and real; by default they
% are assumed to be complex
syms k M positive
% Now, let's apply the inverse Laplace Transform
z1t = ilaplace(expr)
% Let's now assume some numerical values: M = 1kg, k = 55 N/m, z2_0 =
% 5E-3 m and f0 = 2 N. The command below substitutes the numerical values:
z1t = subs(z1t, {M, k, z2_0, f0}, [1,55,5E-3,2]);
% Generally, z1t may not be a simple enough expression. It could be
% simplified in several steps, such as the following:
z1tn = vpa(z1t); % evaluate numerically all coefficients
% Uncomment the following line if using an old version of the symbolic toolbox (2013 or earlier)
% [z1tn, h] = simple(z1tn);
z1tn = vpa(simplify(z1tn),3) % display the result with 3 significant digits
% Now we can plot the response of the system
t = 0:0.01:6; % time between 0 and 6 with interval 0.01 between points
z1 = subs(z1tn,'t',t); % finds the numerical values of z1(t)
% You may have trouble with the line above if you use an
% old version of Matlab. If you have Maple, you could use the following
% line instead:
% z1 = subs(maple('evalf',z1t),'t',t);
plot(t,z1) % plots the graph
title('z_1(t) for t > 0.'); % writes the title of the graph
ylabel('z_1(t) [m]'); % lables the ordinate
xlabel('t [s]'); % labels the abscissa
echo off
% He is wise in heart, and mighty in strength: who hath hardened himself
% against him, and hath prospered? Job 9:4