Numerical Solution of Differential Equations

In engineering, differential equations are everywhere. Their most significant use is as a mathematical description of a dynamic system, such as a solid mechanics system, electronics system, or a fluid mechanics system.

We can solve a few special types of differential equations analytically, but many we cannot. In these cases, especially, we will turn to the numerical solution of differential equations. Usually, analytic solutions are superior to numerical solutions because they are more general (they are usually a family of solutions) and provide more insight. When analytic solutions are unknown or inconvenient to find, numerical solutions can provide insight.

Writing a higher-order differential equation as a first-order system

MATLAB can only solve systems of ordinary differential equations that are first-order. Fortunately, we can rewrite any $n$th-order ode

as a system of first-order odes. Let the the dependent variable be $y$ and the independent variable $t$ (usually “time”). We can define $n$ new variables $y_1(t)$, $y_2(t)$, …, $y_n(t)$ as follows \begin{align} y_1(t) &= y(t) \\ y_2(t) &= \frac{dy}{dt} \\ y_3(t) &= \frac{d^2y}{dt^2} \\ \vdots \\ y_n(t) &= \frac{d^{n-1}y}{dt^{n-1}}. \end{align} Define a vector function from these:

We can now write our original differential equation as the system

where $g$ is essentially just $f$ from above. Note that this technique is valid for both linear and nonlinear ordinary differential equations.

Using MATLAB to solve differential equations

Once we have our system of first-order ordinary differential equations, we have to set them up to be solved in MATLAB. There are several different numerical solvers in MATLAB. Which is best in a given situation is a subject for a course on numerical methods, but generally we begin with a Runge-Kutta solver like ode45.

Solver inputs

All solvers have three basic inputs: odefun, the differential equation function; tspan, a vector with the range of time over which to solve; and y0, a vector of initial conditions. The first of these is the most challenging.

The odefun can be named anything, but must have as its arguments (at least) t (time) and y (the dependent variable array) and must have as its output the time-derivative of y.

Solver outputs

The solver outputs a time-array t and solution-array y. These can easily be plotted using plot(t,y), for instance.

Example

Derive a dynamical model of the horizontal motion of a motor of mass $m$ mounted on a workbench with legs that have combined horizontal stiffness $k$ and combined damping coefficient $B$. Model the motor’s running as a sinusoidal force input horizontally. Simulate several cases in MATLAB.

Some notes on this example.

Another example

Some notes on the electrical example.

State space

State space models can be solved using the solution method described above, but MATLAB provides an easier method. A special state space object is created using the ss command. Several functions can be used on this object, some of the most common of which are step, lsim, eig, and bodeplot. We will explore these through solving the following example.

Example

Consider a flywheel of moment of inertia $J = 1\ kg\, m^2$ connected to a flexible shaft with rotational spring constant $k = 10^4 N\, m/rad$. Let the bearings have viscous damping coefficient $10 N\, m\, s/rad$. Let the shaft have a torque source. Develop a state space model using the linear graph methods of ME 345, enter it into MATLAB using ss, and simulate the system’s response to a shaft torque input of $T_s = 5\ N\, m$. Assume zeroed initial conditions. Moreover, simulate the system’s response to a sinusoidal input $T_s = 5 \sin{(\omega_n/2\,t)}$, where $\omega_n$ is the natural frequency. Now, double the frequency. Interpret the results in the light of the bode plot. Finally, simulate the system’s response to a sawtooth input at the same amplitude and two frequencies as the sinusoid. Bask in your godlike powers.

Tunable models

In engineering design, it often occurs that we can tune one or more parameter in a system in order to make its dynamic response behave in a desired manner. MATLAB allows us to create state space model objects with tunable parameters using the following syntax (example taken from the MATLAB documentation).

a = realp('a',3);   % tunable parameter
b = realp('b',-5);  % tunable parameter

A = [-1 a;0 a*b];
B = [-3.0;1.5];
C = [0.3 0];
D = 0;

sys = ss(A,B,C,D);  % build state space object

a_vals = [1;10];
b_vals = [-1 -1.5 -2];
[a_grid,b_grid] = ndgrid(a_vals,b_vals);
sys_sample = replaceBlock(sys,'a',a_grid,'b',b_grid);     % sample system
sys_sample.SamplingGrid = struct('a',a_grid,'b',b_grid);  % just labeling

stepplot(sys_sample)  % see what you've done

Example

Let’s say you’re an engineer at a manufacturing plant with a two-shelf table bolted to the floor. The top shelf of the table has a piece of relatively sensitive equipment that needs to be isolated from the vibration. However, currently, the velocity of the sensitive equipment has twice the allowable amplitude. You notice that the table is near heavy machinery that excites the floor with a sinusoidal velocity of $V_s(t) = A \sin{(\omega t)}$, where $A = 1\ m/s$ and $\omega = 50\ rad/s$. It is desirable to try to reduce the velocity of the top shelf of the table by adding mass to it, which is easy to do. Determine by how much its mass $m_2$ should increase. Model the bottom and top shelves as masses $m_1 = m_2 = 100\ kg$. Model the legs as springs with spring constants $k_1 = k_2 = 10^5\ N/m$ and damping $B_1 = B_2 = 100\ N\,s/m$. Can we reduce the vibration speed by half without more than doubling the original mass $m_2$? If not, explore another possibility.