Introduction to Mechanics for Engineersby Andy Ruina & Rudra Pratap2 Free-Body Diagrams 4 Units and estimation

Chapter 3 Computation: Algebra, ODEs, etc.

Download this chapter by section

3.1ODEs: Ordinary Differential Equationsstudy PDF + solutions
3.2Numerical solution of ODEsstudy PDF + solutions
3.3Introductionstudy PDF

Here is a quick review of some of the important math and computation for dynamics. These are things that engineers should just know, like they know that the derivative of $\sin(t)$ is $\cos(t)$ .

3.1 Basic summary of ODEs for mechanics

This section is not an aside. Know it well.

Here are some of the simplest useful ordinary differential equations (ODEs) and their general solutions. Think of $u$ as the distance an object has moved to the right of its ‘home’ at $u=0$ in time $t$ . The velocity and acceleration to the right are the derivatives of $u$ .

position	$\displaystyle u$
velocity	$\displaystyle du/dt=\dot{u}$
acceleration	$\displaystyle d^{2}u/dt^{2}=\ddot{u}.$	(3.1)

If $\dot{u}<0$ the particle is moving to the left. If $\ddot{u}<0$ the particle is accelerating to the left.

In all cases below $A$ and $B$ are constants and $\lambda$ is a positive constant. $C_{1}$ , $C_{2}$ , $C_{3}$ , and $C_{4}$ are arbitrary (undetermined) constants in the solutions that get pinned down (determined) by fixing the initial conditions.

This section is not about how to derive the solutions, it is about how to know them.

Constant position: ${\dot{u}=0}\ \Rightarrow u=C_{1}$

Filename:tfigure8-rel-ang-vel

Figure 3.1: Constant position.

$\dot{u}=0$ means that the velocity is zero. This equation would arise in dynamics if a particle has no initial velocity and no force is applied to it. The particle doesn’t move. Its position must be constant. But it could be anywhere, say at position $C_{1}$ . Hence the general solution $u=C_{1}$ , as can be found by direct integration.

Constant velocity: ${\dot{u}=A}\ \Rightarrow\ {u=At+C_{1}}$

Filename:tfigure8-ang-vel-ex

Figure 3.2: Constant velocity.

$\dot{u}=A$ means the object has constant speed. This equation describes the motion of a particle that starts with speed $v_{0}=A$ and because it has no force acting on it continues to move at constant speed. How far does it go in time $t$ ? It goes $v_{0}t$ . Where was it at time $t=0$ ? It could have been anywhere then, say $C_{1}$ . So where is it at time $t$ ? It’s at its original position plus how far it has moved, $u=v_{0}t+C_{1}$ , as can also be found by direct integration.

Zero acceleration: ${\ddot{u}=0}\ \Rightarrow\ u=C_{1}+C_{2}t$

Filename:tfigure8-ang-vel

Figure 3.3: Zero acceleration.

$\ddot{u}=0$ means the acceleration is zero. That is, the rate of change of velocity is zero. This constant-velocity motion is the general equation for a particle with no force acting on it. The velocity, if not changing, must be constant. What constant? It could be anything, say $C_{2}$ . Now we have the same situation as in case (b). So the position as a function of time is anything consistent with an object moving at constant velocity: $u=C_{1}+C_{2}t$ , where the constants $C_{1}$ and $C_{2}$ depend on the initial position and initial velocity. If you know that the position at $t=0$ is $u_{0}$ and the velocity at $t=0$ is $v_{0}$ , then the position is $u=u_{0}+v_{0}t$ .

Constant acceleration: ${\ddot{u}=A}\indent\Rightarrow\ \mbox{ }u=At^{2}/2+C_{1}t+C_{2}$

Filename:tfigure8-ang-accel

Figure 3.4: Constant acceleration.

This constant acceleration $A$ , constant rate of change of velocity, is the classic (all-too-often studied) case. This situation arises for vertical motion of an object in a constant gravitational field as well as in problems of constant acceleration or deceleration of vehicles. The velocity increases in proportion to the time that passes. The change in velocity in a given time is thus $A t$ and the velocity is $v={\dot{u}}=v_{0}+At$ (given that the velocity was $v_{0}$ at $t=0$ ). Because the velocity is increasing constantly over time, the average velocity in a trip of length $t$ occurs at $t/2$ and is $v_{0}+At/2$ . The distance traveled is the average velocity times the time of travel so the distance of travel is $t\cdot(v_{0}+At/2)=v_{0}t+At^{2}/2$ . The position is the position at $t=0$ , $u_{0}$ , plus the distance traveled since time zero. So $u=u_{0}+v_{0}t+At^{2}/2=C_{2}+C_{1}t+At^{2}/2$ . This solution can also be found by direct integration.

Exponential growth: ${\dot{u}=\lambda u}\ \Rightarrow\ {u=C_{1}e^{\lambda t}}$

Filename:tfigure8-ang-accel-ex

Figure 3.5: Exponential growth.

The displacement $u$ grows in proportion to its present size. This equation describes the initial falling of an inverted pendulum in a thick viscous fluid. The bigger the $u$ , the faster it moves. Such situations are called exponential growth (as in population growth or monetary inflation) for a good mathematical reason. The solution $u$ is an exponential function of time: $u(t)=C_{1}e^{\lambda t}$ , as can be found by separating variables or guessing.

Exponential decay: ${\dot{u}=-\lambda u}\ \Rightarrow\ {u=C_{1}e^{-\lambda t}}$

Filename:tfigure8-ang-accel1

Figure 3.6: Exponential decay.

The smaller $u$ is, the more slowly it gets smaller. $u$ gradually tapers towards nothing: $u$ decays exponentially. The solution to the equation is: $u(t)=C_{1}e^{-\lambda t}$ . This expression is essentially the same equation as in (e) above.

Acceleration proportional to position: ${\ddot{u}=\lambda^{2}u}$ $\Rightarrow{u=C_{1}e^{\lambda\ t}\ +\ C_{2}e^{-\lambda\ t}}\hfill\break\hbox{}% \hskip 28.90755pt\textrm{or}\Rightarrow{u=C_{3}\cosh(\lambda\ t)\ +\ C_{4}% \sinh(\lambda\ t)}$

Filename:sfig8-2-3

Figure 3.7: Acceleration proportional to position.

Note, $\sinh$ and $\cosh$ are just combinations of exponentials. For $\ddot{u}=\lambda^{2}u$ , the point accelerates more and more away from the origin in proportion to the distance from the origin. This equation describes the initial falling of a nearly vertical inverted pendulum when there is no friction. Most often, the solution of this equation gives roughly exponential growth. The pendulum accelerates away from being upright. The reason there is also an exponentially decaying solution to this equation is a little more subtle to understand intuitively: if a not quite upright pendulum is given just the right initial velocity it will slowly approach becoming just upright with an exponentially decaying displacement. This decaying solution is not easy to see experimentally because, without the perfect initial condition, the exponentially growing part of the solution eventually dominates and the pendulum accelerates away from being just upright.

Harmonic oscillator: $\ddot{u}=-\lambda^{2}u\mbox{~{}~{}~{}or~{}~{}~{}}\ddot{u}+\lambda^{2}u=0$ $\Rightarrow\ \mbox{ }u=C_{1}\sin(\lambda\ t)+C_{2}\cos(\lambda\ t)$ .

Filename:sfig8-2-3a

Figure 3.8: Harmonic oscillator.

This equation describes a mass that is restrained by a spring which is relaxed when the mass is at $u=0$ . When $u$ is positive, $\ddot{u}$ is negative. That is, if the particle is on the right side of the origin it accelerates to the left. Similarly, if the particle is on the left it accelerates to the right. In the middle, where $u=0$ , it has no acceleration, so it neither speeds up nor slows down in its motion whether it is moving to the left or the right. So the particle goes back and forth: its position oscillates. A function that correctly describes this oscillation is $u=\sin(\lambda\,t)$ , that is, sinusoidal oscillations. The oscillations are faster if $\lambda$ is bigger. Another solution is $u=\cos(\lambda\,t)$ . The general solution is $u=C_{1}\sin(\lambda\,t)+C_{2}\cos(\lambda\,t)$ . A plot of this function reveals a sine wave shape for any value of $C_{1}$ or $C_{2}$ , although the phase depends on the relative values of $C_{1}$ and $C_{2}$ . The equation $\ddot{u}=-\lambda^{2}u$ or $\ddot{u}+\lambda^{2}u=0$ is called the ‘harmonic oscillator’ equation and is important in almost all branches of science. The solution may be found by guessing or other means (which are usually guessing in disguise). In the context of this equation, $\lambda$ is called the (angular) frequency of oscillation.

More ODEs

Besides the ODEs listed here there are a few others that are often solved by hand rather than with numerical simulation. Most famously there is the damped oscillator equation $A\ddot{u}+B\dot{u}+Cu=0$ discussed in Box 11.1 on page 11.1.

With the exception of the damped or forced oscillator, most engineers now-a-days will use numerical integration if they want to solve an ODE not in this appendix.

Learn, memorize, know, whatever

Conversely,

Anyone competent at dynamics knows all the equations and solutions on these few pages outside, inside, and inside-out.

Whether you have or have not learned these in a calculus course you should learn them now every way possible.

SAMPLE 3.1 A simple first order differential equation: Show that $y(t)=vt$ and $y(t)=C+vt$ are both solutions of the first order differential equation $\dot{y}=v$ where $v$ and $C$ are constants. Which of the two solutions is a general solution of the given differential equation and why?

Solution Any function $y(t)$ will be a solution of the given differential equation if it satisfies the differential equation. That is, if we plug the function into the differential equation, the left hand side must equal the right hand side. So, let us check if the given solutions satisfy the given differential equation:

(a)

$y(t)=vt$ : Differentiating $y(t)$ we get, $\dot{y}=\frac{d}{dt}(vt)=v$ . Substituting in the given differential equation, we get,

$\underbrace{v}_{\dot{y}}=v.$

Thus the left hand side is equal to the right hand side and hence the differential equation is satisfied. Therefore, $y=vt$ is a legitimate solution of the given differential equation.
(b)

$y(t)=C+vt$ : Differentiating $y(t)$ we get, $\dot{y}=\frac{d}{dt}(C+vt)=0+v=v$ . Substituting in the given differential equation, $\dot{y}=v$ , we get

$\underbrace{v}_{\dot{y}}=v.$

Thus again, the left hand side is equal to the right hand side and hence the differential equation is satisfied. Therefore, $y(t)=C+vt$ is a legitimate solution of the given differential equation.

General Solution: The general solution of a differential equation refers to a solution that has enough free parameters (unknown constants) to accommodate any given initial conditions (for Initial-Value-Problems) or boundary conditions (for Boundary-Value-Problems) and generate a unique solution for the given problem. In dynamics, we have mostly initial-value-problems where initial conditions (e.g., position, displacement, or velocity at $t=0$ ) are specified. We can then plug the given initial conditions in the general solution and set $t=0$ to find the unknown constants. The resulting solution then becomes a unique solution.

Here, $y(t)=C+vt$ is the general solution as it contains one free parameter $C$ that can accommodate any given initial condition $y(t=0)=y_{0}$ (say). By substituting this initial condition in the solution ( $y(t)=C+vt$ ), we get,

y_{0}=C+v\overbrace{t}^{0}\hbox to36.135pt{\hfil$\Rightarrow$\hfil}C=y_{0}.

Thus, the solutions to the differential equation, $\dot{y}=v$ , $y(0)=y_{0}$ , is

y(t)=y_{0}+vt.

Clearly, the first given solution $y=vt$ is a subset of the general solution, with $C=0$ , and thus it represents the solution of the given differential equation with the initial condition $y(0)=0$ . It is, therefore, a unique solution of the problem $\dot{y}=v$ , $y(0)=0$ , and not a general solution of the differential equation $\dot{y}=v$ .

Answer: $\mbox{General solution: }y(t)=C+vt$

SAMPLE 3.2 Initial conditions and constants of integration: Let $y(t)=C_{1}+C_{2}t$ be the solution of the differential equation $\ddot{y}=0$ . Find the solutions corresponding to the following initial conditions and sketch $y(t)$ vs $t$ in each case.

1.

$y(0)=1\,\mathrm{m}$ and $\dot{y}(0)=2\,\mathrm{m}/\mathrm{s}$ .
2.

$y(0)=0$ and $\dot{y}(0)=5\,\mathrm{m}/\mathrm{s}$ .

Solution

The general solution has two free constants $C_{1}$ and $C_{2}$ . We can find their values by substituting the given set of initial conditions and thus determine the corresponding solutions. From the general solution, we have,

	$\displaystyle y$	$\displaystyle=$	$\displaystyle C_{1}+C_{2}t$		(3.2)
	$\displaystyle\dot{y}$	$\displaystyle=$	$\displaystyle C_{2}.$		(3.3)

1.

$y(0)=1\,\mathrm{m};\dot{y}(0)=2\,\mathrm{m}/\mathrm{s}$ :

Substituting these values in eqn. (3.2) and eqn. (3.3), we get,

$\displaystyle 1\,\mathrm{m}$ $\displaystyle=$ $\displaystyle C_{1}$

$\displaystyle\hbox to36.135pt{and\hfil}2\,\mathrm{m}/\mathrm{s}$ $\displaystyle=$ $\displaystyle C_{2}.$

Filename:sfig8-2-3b

Figure 3.9: Plot of $y(t)=1\,\mathrm{m}+(2\,\mathrm{m}/\mathrm{s})t$ as the solution of $\ddot{y}=0$ ; $y(0)=1\,\mathrm{m}$ and $\dot{y}(0)=2\,\mathrm{m}/\mathrm{s}$ .

Thus, the corresponding solution is,

$y(t)=1\,\mathrm{m}+(2\,\mathrm{m}/\mathrm{s})t.$

We can now sketch $y(t)$ vs $t$ as shown in fig. 3.9.

Answer: $y(t)=1\,\mathrm{m}+(2\,\mathrm{m}/\mathrm{s})t$
2.

$y(0)=0;\dot{y}(0)=5\,\mathrm{m}/\mathrm{s}$ :

Again substituting the given values in eqn. (3.2) and eqn. (3.3), we get,

$\displaystyle 0$ $\displaystyle=$ $\displaystyle C_{1}$

$\displaystyle\hbox to36.135pt{and\hfil}5\,\mathrm{m}/\mathrm{s}=C_{2}.$

Filename:sfig8-3-1

Figure 3.10: Plot of $y(t)=(5\,\mathrm{m}/\mathrm{s})t$ as the solution of $\ddot{y}=0$ ; $y(0)=0$ and $\dot{y}(0)=5\,\mathrm{m}/\mathrm{s}$ .

Thus the corresponding solution is

$y(t)=(5\,\mathrm{m}/\mathrm{s})t.$

This solution is shown in fig. 3.10.

Answer: $y(t)=(5\,\mathrm{m}/\mathrm{s})t$

Comments: The differential equation $\ddot{y}=0$ represents a zero acceleration, or alternatively, a constant speed motion in the $y$ direction. Therefore, the displacement $y$ is a linear function of $t$ as shown in the two graphs. The intercept, $y(t=0)$ , represents the initial displacement and the slope $\frac{dy}{dt}$ represents the constant velocity. In fig. 3.9, the initial displacement is 1 $\,\mathrm{m}$ which is shown as the intercept (at $t=0$ ) and the slope is the given velocity 2 $\,\mathrm{m}/\mathrm{s}$ . In fig. 3.10, the initial position is $y(0)=0$ and hence the graph passes through the origin. The slope is the given velocity $\dot{y}(0)=5\,\mathrm{m}/\mathrm{s}$ (since there is no acceleration, the initial velocity remains constant throughout the motion).

SAMPLE 3.3 Constant speed motion in 1-D again: A particle moves in the $x$ -direction with constant speed. What is the differential equation that represents this motion? Find the motion of the particle, i.e., find the expression for $x(t)$ , given the position of the particle at three instances: $x(0)=x_{0}$ , $x(t_{1})=x_{1}$ and $x(t_{2})=x_{2}$ . Do you need $x(t)$ at three instances to determine $x(t)$ completely? What could be alternate information that will help you specify the motion completely?

Solution

Let $x$ be the displacement. Then $\dot{x}$ represents speed and $\ddot{x}$ represents acceleration. Since it is motion in 1-D, we will not worry about the vector notation. The given information about constant speed can be represented by either of the two differential equations:

1.

$\dot{x}=C$ where $C$ is the given constant speed, or
2.

$\ddot{x}=0$ , i.e., the acceleration is zero because the speed is constant.

In the first case, the solution of the differential equation is,

x(t)=Ct+C_{1}

where $C_{1}$ is the constant of integration. We find this solution by direct integration as follows.

\dot{x}=C\hbox to36.135pt{\hfil$\Rightarrow$\hfil}dx=Cdt\hbox to36.135pt{\hfil% $\Rightarrow$\hfil}\int dx=\int Cdt\hbox to36.135pt{\hfil$\Rightarrow$\hfil}x=% Ct+C_{1}.

In the second case, we find the solution by integrating twice:

	$\displaystyle\ddot{x}$	$\displaystyle=0\hbox to36.135pt{\hfil$\Rightarrow$\hfil}\int d\dot{x}=\int 0\;% dt\hbox to36.135pt{\hfil$\Rightarrow$\hfil}\dot{x}=C_{1}$
	$\displaystyle\dot{x}$	$\displaystyle=C_{1}\hbox to36.135pt{\hfil$\Rightarrow$\hfil}\int dx=\int C_{1}% dt\hbox to36.135pt{\hfil$\Rightarrow$\hfil}x=C_{1}t+C_{2}.$

Filename:sfig8-3-1a

Figure 3.11: The graph of the general solution

x(t)=C_{1}t+C_{2}

Both solutions above represent the same motion, as they must. The displacement $x(t)$ varies linearly with time $t$ as we would expect in a constant speed motion. This solution is represented by a line in $t$ - $x$ coordinates as shown in fig. 3.11 with $C_{1}$ as the slope and $C_{2}$ as the intercept. Now, from the given information about values of $x$ at three time instances, it is easy to determine $C_{1}$ and $C_{2}$ . We can do it multiple ways:

Using the general solution:

Substituting the values of $x$ for the three given time instances, $x(t=0)=x_{0},x(t=t_{1})=x_{1}$ and $x(t=t_{2})=x_{2}$ , we get

	$\displaystyle x_{0}$	$\displaystyle=C_{1}\cdot 0+C_{2}\hbox to36.135pt{\hfil$\Rightarrow$\hfil}C_{2}% =x_{0}$
	$\displaystyle x_{1}$	$\displaystyle=C_{1}t_{1}+x_{0}$
	$\displaystyle x_{2}$	$\displaystyle=C_{1}t_{2}+x_{0}$

From the last two equations, we get $C_{1}=\frac{x_{2}-x_{1}}{t_{2}-t_{1}}$ and thus the solution that satisfies the given conditions is,

x(t)=x_{0}+\left(\frac{x_{2}-x_{1}}{t_{2}-t_{1}}\right)t.

From geometry of the function $x(t)$ :

Filename:sfig8-5-wiper

Figure 3.12: Given values of

x(t)

determine the slope and the intercept

We can mark the given values of $x(t)$ for the three specified time instants on the graph of $x(t)$ as shown in fig. 3.12. From the given values, it is obvious that $C_{2}=x_{0}$ and the slope of the line $C_{1}$ is $\Delta x/\Delta t=(x_{2}-x_{1})/(t_{2}-t_{1})$ , and hence the equation of the line becomes,

x(t)=x_{0}+\left(\frac{x_{2}-x_{1}}{t_{2}-t_{1}}\right)t

which is the solution of the differential equation that satisfies all given conditions.

Answer: $x(t)=x_{0}+\left(\frac{x_{2}-x_{1}}{t_{2}-t_{1}}\right)t$

SAMPLE 3.4 Find the solution of the differential equation

\ddot{x}+\alpha\dot{x}=0;\mbox{ with initial conditions }x(0)=0\mbox{ and }% \dot{x}(0)=u_{0}.

Solution

The given equation is a second order differential equation that involves derivative terms of $x$ but not $x$ itself. We can solve this equation in two steps. First we solve for $\dot{x}$ and then we solve for $x$ . Let us introduce a new variable $y$ by renaming $\dot{x}$ as $y$ , that is, let $y=\dot{x}$ . Then the given equation becomes a first order equation in $y$ :

\dot{y}+\alpha y=0,\mbox{ or }\dot{y}=-\alpha y.

We can solve this equation by integration as follows:

$\displaystyle\frac{dy}{dt}$	$\displaystyle=$	$\displaystyle-\alpha y$
$\displaystyle\hbox to36.135pt{\hfil$\Rightarrow$\hfil}\frac{dy}{y}$	$\displaystyle=$	$\displaystyle-\alpha dt$
$\displaystyle\hbox to36.135pt{\hfil$\Rightarrow$\hfil}\int\frac{1}{y}dy$	$\displaystyle=$	$\displaystyle-\alpha\int dt$
$\displaystyle\hbox to36.135pt{\hfil$\Rightarrow$\hfil}\ln y$	$\displaystyle=$	$\displaystyle-\alpha t+C$
$\displaystyle\hbox to36.135pt{\hfil$\Rightarrow$\hfil}y$	$\displaystyle=$	$\displaystyle e^{-\alpha t+C}$
$\displaystyle\mbox{or, }y$	$\displaystyle=$	$\displaystyle C_{1}e^{-\alpha t}$

where the last expression is obtained by substituting $e^{C}$ with another constant $C_{1}$ (note that $e^{-\alpha t+C}=e^{C}e^{-\alpha t}$ ). So, we have found the solution for $y(t)$ . But $y(t)=\dot{x}$ , so we now have another differential equation to solve:

\dot{x}=C_{1}e^{-\alpha t}.

We can solve this equation by direct integration as follows:

$\displaystyle\frac{dx}{dt}$	$\displaystyle=$	$\displaystyle C_{1}e^{-\alpha t}$
$\displaystyle\hbox to36.135pt{\hfil$\Rightarrow$\hfil}\int dx$	$\displaystyle=$	$\displaystyle C_{1}\int e^{-\alpha t}dt$
$\displaystyle\hbox to36.135pt{\hfil$\Rightarrow$\hfil}x$	$\displaystyle=$	$\displaystyle C_{1}\frac{e^{-\alpha t}}{-\alpha}+C_{2}$
$\displaystyle\mbox{or, }x$	$\displaystyle=$	$\displaystyle C_{2}-\frac{C_{1}}{\alpha}e^{-\alpha t}.$

We now have expressions for both $x(t)$ and $\dot{x}(t)\equiv y(t)$ . We can substitute the given initial conditions and determine the constants $C_{1}$ and $C_{2}$ :

$\displaystyle\dot{x}(0)$	$\displaystyle\equiv$	$\displaystyle y(0)=u_{0}$
$\displaystyle\hbox to36.135pt{\hfil$\Rightarrow$\hfil}u_{0}$	$\displaystyle=$	$\displaystyle C_{1}e^{-\alpha 0}=C_{1}$
and
$\displaystyle x(0)$	$\displaystyle=$	$\displaystyle 0$
$\displaystyle\hbox to36.135pt{\hfil$\Rightarrow$\hfil}0$	$\displaystyle=$	$\displaystyle C_{2}-\frac{u_{0}}{\alpha}e^{-\alpha 0}$
$\displaystyle\hbox to36.135pt{\hfil$\Rightarrow$\hfil}C_{2}$	$\displaystyle=$	$\displaystyle\frac{u_{0}}{\alpha}.$

Filename:sfig8-5-wiper-a

Figure 3.13: Graph of

x(t)=\frac{u_{0}}{\alpha}\left(1-e^{-\alpha t}\right)

. The system reaches the equilibrium value of

x=u_{0}/\alpha

asymptotically as

t\rightarrow\infty

So, now substituting the values of $C_{1}$ and $C_{2}$ , we get the final solution:

x(t)=\frac{u_{0}}{\alpha}\left(1-e^{-\alpha t}\right).

The graph of $x(t)$ is shown in fig. 3.13.

Answer: $x(t)=\frac{u_{0}}{\alpha}\left(1-e^{-\alpha t}\right)$

SAMPLE 3.5 Units of coefficients: The equation of motion of a mechanical system is given as follows:

C_{1}\ddot{y}+C_{2}\dot{y}+C_{3}y=C_{4}

where $y$ has units of displacement (e.g., meters) and $C_{1}$ has units of mass (e.g., kg). What must be the units of $C_{2}$ , $C_{3}$ and $C_{4}$ in SI system?

Solution

The units of coefficients in any equation, differential or algebraic, have to be such that each term in the summation has the same units. So, in the given differential equation, $C_{1}\ddot{y}$ , $C_{2}\dot{y}$ , $C_{3}y$ , and $C_{4}$ must have the same units. Now, $y$ has the units of displacement (i.e., distance or length) and $C_{1}$ has the units of mass. So, in order to find the units of each term, let us first find the units of $\dot{y}$ and $\ddot{y}$ :

	$\displaystyle\dot{y}$	$\displaystyle=$	$\displaystyle\frac{dy}{dt}=\frac{\,\mathrm{m}}{\,\mathrm{s}}=\,\mathrm{m}/% \mathrm{s}$
	$\displaystyle\ddot{y}$	$\displaystyle=$	$\displaystyle\frac{d\dot{y}}{dt}=\frac{\,\mathrm{m}/\mathrm{s}}{\,\mathrm{s}}=% \,\mathrm{m}/\mathrm{s}^{2}.$

So, the units of the first term are:

C_{1}\ddot{y}=\,\mathrm{kg}\cdot\,\mathrm{m}/\mathrm{s}^{2}=\,{\mbox{N}},

that is, unit of force. For the equation to hold good, each term, therefore, must have units of force. So, we can now find the units of all coefficients as follows:

$\displaystyle C_{2}\dot{y}$	$\displaystyle=$	N
$\displaystyle\hbox to36.135pt{\hfil$\Rightarrow$\hfil}C_{2}$	$\displaystyle=$	$\displaystyle\frac{\,{\mbox{N}}}{\,\mathrm{m}/\mathrm{s}}=\frac{\,{\mbox{N}}% \cdot\,\mathrm{s}}{\,\mathrm{m}}=\frac{\,\mathrm{kg}}{\,\mathrm{s}},$
$\displaystyle C_{3}y$	$\displaystyle=$	N
$\displaystyle\hbox to36.135pt{\hfil$\Rightarrow$\hfil}C_{3}$	$\displaystyle=$	$\displaystyle\frac{\,{\mbox{N}}}{\,\mathrm{m}},\mbox{ and }$
$\displaystyle C_{4}$	$\displaystyle=$	$\displaystyle\,{\mbox{N}}.$

While units of $C_{4}$ are clearly that of force, the units of $C_{2}$ and $C_{3}$ are somewhat unusual and worth thinking about. $C_{3}$ has the units of stiffness, a quantity that measures the resistance of a body to deformation, and that is why its units are $\,{\mbox{N}}/\,\mathrm{m}$ . You may be familiar with this unit as the unit of spring stiffness. The units of $C_{2}$ are even more unusual. You can think of these units as $\,\mathrm{kg}/\,\mathrm{s}$ or $\,{\mbox{N}}/(\,\mathrm{m}/\mathrm{s})$ . $C_{2}$ has the units of damping, a quantity that measures the resistance to motion of a body in a viscous medium in terms of force per unit speed.

Answer: $C_{2}$ : $\,\mathrm{kg}$ / $\,\mathrm{s}$ , $C_{3}$ : N/ $\,\mathrm{m}$ , and $C_{4}$ : N

SAMPLE 3.6 Initial conditions dictate solutions: Consider the second order differential equation $\ddot{y}+\omega^{2}y=0$ . Show that $y(t)=A\sin\omega t$ , $y(t)=B\cos\omega t$ , and $y(t)=A\sin\omega t+B\cos\omega t$ are all valid solutions of the differential equation. Which of the three solutions will satisfy the following initial conditions?

1.

$y(0)=y_{0}$ and $\dot{y}(0)=0$ .
2.

$y(0)=0$ and $\dot{y}(0)=v_{0}$ .
3.

$y(0)=y_{0}$ and $\dot{y}(0)=v_{0}$ .

Solution

In order to show that a given function is a possible solution of a differential equation, we just need to plug the function into the differential equation and see if it satisfies the equation. We can easily do this for the three given functions:

•

$y(t)=A\sin\omega t$ : We need to find $\ddot{y}$ and plug into the given differential equation.

$\dot{y}=-A\omega\cos\omega t\hbox to36.135pt{\hfil$\Rightarrow$\hfil}\ddot{y}=% -A\omega^{2}\sin\omega t\\$

Therefore,

$\ddot{y}+\omega^{2}y=\underbrace{-A\omega^{2}\sin\omega t}_{\ddot{y}}+\omega^{% 2}\underbrace{A\sin\omega t}_{y}=0.$

Thus $y=A\sin\omega t$ satisfies the given differential equation and hence it is a solution.

Filename:sfig8-5-2disks

Figure 3.14: Graph of the three solutions: (a) $y(t)=A\sin\omega t$ , (b) $y(t)=B\cos\omega t$ , and (c) $y(t)=A\sin\omega t+B\cos\omega t$
•

$y(t)=B\cos\omega t$ : In this case, $\ddot{y}=-B\omega^{2}\cos\omega t$ , and substituting $y$ and $\ddot{y}$ in the given differential equation, we again find that the equation is satisfied. Hence $y(t)=B\cos\omega t$ is also a solution.
•

$y(t)=A\sin\omega t+B\cos\omega t$ : Since each of the two terms here satisfies the equation individually, their sum will obviously also satisfy the equation ( $0+0=0$ ). Hence this is also a solution of the given differential equation. In fact, this is the most general solution.

The three solutions are shown in fig. 3.14. Now let us examine which ones can satisfy which initial conditions.

1.

$y(0)=y_{0}$ and $\dot{y}(0)=0$ : From fig. 3.14(a), it is clear that $y(t)=A\sin\omega t$ cannot support non-zero $y$ at $t=0$ , and hence cannot satisfy $y(0)=y_{0}$ . We also see that this function has a non-zero slope at $t=0$ , i.e., $\dot{y}(0)\neq 0$ and hence it cannot support the second initial condition either. On the other hand, fig. 3.14(b) can easily satisfy both initial conditions if we set $B=y_{0}$ . Similarly, the third solution shown graphically in fig. 3.14(c) can also satisfy the given initial conditions if we set $A=0$ and $B=y_{0}$ (then this solution reduces to the second solution).
2.

$y(0)=0$ and $\dot{y}(0)=v_{0}$ : Similar arguments as above lead us to conclude that $y(t)=A\sin\omega t$ satisfies the given initial conditions if we set $A=v_{0}/\omega$ , and $y(t)=A\sin\omega t+B\cos\omega t$ also satisfies the given initial conditions if we set $B=0$ and $A=v_{0}/\omega$ .
3.

$y(0)=y_{0}$ and $\dot{y}(0)=v_{0}$ : Now we have both initial displacement and initial velocity as non-zero quantities and the graphs in fig. 3.14(a) and (b) make it obvious that neither of these two functions can satisfy these initial conditions. However, the third solution $y(t)=A\sin\omega t+B\cos\omega t$ can certainly satisfy these conditions if we set $A=v_{0}/\omega$ and $B=y_{0}$ . Thus the solution will be

$y(t)=y_{0}\cos\omega t+\frac{v_{0}}{\omega}\sin\omega t.$

Thus the third solution can satisfy any given initial conditions ( $y_{0}$ can take any value including zero and so can $v_{0}$ ). This is why it is called the general solution.

Problems for 3.1 Basic computation: linear equations and plotting

3.1.1

Write the equations, and general solutions, for as many problems as you can.

3.2 Numerical solution of ODEs

Learn and master Euler’s method. Write it yourself. After you master the concepts, write it using a separate function for the Right-Hand-Side (RHS) file. And a separate function for the Euler integrator.

Once that works well for you. Then learn the syntax of a numerical method with your favorite software, for example MATLAB ’s ODE45.

SAMPLE 3.7 A simple first order ODE on the computer: Solve $\dot{x}-x=t$ , with initial condition $x(0)=0$ , numerically using Euler’s method. Plot the solution for $0\leq t\leq 2$ where $t$ is dimensionless time.

Solution As explained in the text, Euler’s method is the simplest, although not very efficient or stable, numerical method to solve ODEs. Here, we use Euler’s method to solve the given differential equation. The conceptual steps are as follows and the actual code is implemented in MATLAB.

Step-1:

Write the differential equation as a set of (if applicable) first order equations in the form $\dot{x}=f(x,t)$ :

\dot{x}=x+t

Step-2:

Specify initial conditions, time interval for solution, and time step for advancing the solution.

$\displaystyle x(0)$	$\displaystyle=0$	(given)
$\displaystyle t_{\rm span}$	$\displaystyle=[0\quad 2]$	(given)
$\displaystyle h(=\Delta t)$	$\displaystyle=0.002$	(our choice)

Step-3:

Implement Euler’s method for computing $x_{n+1}$ given $x_{n},t_{n}$ and $h$ , and advance the time instance to $t_{n+1}=t_{n}+h$ .

	$\displaystyle x_{n+1}$	$\displaystyle=$	$\displaystyle x_{n}+\dot{x}(x_{n},t_{n})h$
	$\displaystyle t_{n+1}$	$\displaystyle=$	$\displaystyle t_{n}+h$

Step-4:

Store solution at each time step in an array and plot the solution as $t$ vs $x$ .

Here is the MATLAB code with comments that solves the given ODE with the given initial condition.

% SCRIPT to solve xdot = x+t with x(0)=0 for t up to 2 sec
%----------------------------------------------------------
x0=0; t0=0;              % specify initial conditions
tf=2;                    % specify final time
x=x0;  t=t0;             % initialize arrays for x and t
h=0.002;                 % specify time step
nsteps = (tf-t0)/h;      % calculate number of time steps
for n = 1:nsteps         % compute solutions in a loop
  xdot = x(n) + t(n);    % calculate the derivate at (tn,xn)
  x(n+1) = x(n) + h*xdot;% calculate x(n+1)
  t(n+1) = t(n) + h;     % advance tn to t(n+1)
end
plot(t,x), xlabel(’t’), ylabel(’x’)

Filename:sfig8-4-1

Figure 3.15: Numerical solution of

\dot{x}=x+t

x(0)=0

for

t\leq 2

using Euler’s method in MATLAB.

The result of running this code in MATLAB is shown in fig. 3.15. Please note that this is a plain vanilla code that does not take advantage of several nice features of MATLAB. It is meant to illustrate how we can use Euler’s method in a straightforward manner to solve a differential equation.

SAMPLE 3.8 Numerical solution of a set of ODEs: Solve the following set of two first order differential equation using Euler’s method:

	$\displaystyle\dot{x}$	$\displaystyle=y$
	$\displaystyle\dot{y}$	$\displaystyle=-\alpha x$

with the initial conditions $x(0)=1$ and $y(0)=0$ . Write separate codes for implementing the solution algorithm (e.g., Euler’s method) and the given differential equation, so that you can use the method of solution with any set of differential equations.

Solution Unlike the previous sample problem, we now have a set of two first order equations. We can put them in an array, say $\dot{z}$ , and solve for $z$ where $x$ and $y$ are treated as two separate elements of $z$ , i.e.,

z=\left(\begin{array}[]{c}x\\ y\end{array}\right)\hbox to36.135pt{\hfil$\Rightarrow$\hfil}\dot{z}\equiv\left% (\begin{array}[]{c}\dot{x}\\ \dot{y}\end{array}\right)=\left(\begin{array}[]{c}\dot{y}\\ -\alpha\dot{x}\end{array}\right),

and the initial conditions for $z$ can be written as,

z(0)=\left(\begin{array}[]{c}x(0)\\ y(0)\end{array}\right)=\left(\begin{array}[]{c}1\\ 0\end{array}\right).

Let us write a MATLAB code where the equation is coded in a separate function and Euler’s method is coded in another function assuming that we are solving for a set of first order equations.

•

Coding the differential equation: Let us write a function that takes an array variable $z$ , time $t$ , and additional problem parameter $p$ (e.g., $\alpha$ in our equation), as the input, computes $\dot{z}$ inside the function and puts it out as the output (see schematic diagram in fig. 3.16). Here is the MATLAB code:
```
function [ zdot ] = MyDIffEqn(z,t)
% MyDiffEqn takes the value of array z at time t
%  and outputs the derivative zdot

alpha = 1;          % specify problem parameters
x = z(1); y = z(2); % unpack z
xdot = y            % code the given equation
ydot = -alpha*x
zdot = [xdot; ydot] % pack the derivatives in zdot
end                 % end of function
```
Filename:sfig8-4-1a

Figure 3.16: Schematic of coding a differential equation as a MATLAB function
•

Coding Euler’s method: We already coded Euler’s method in the previous sample. We just need to convert it into an independent function which takes an array $z_{n}$ at time $t_{n}$ , calls the coded differential equation function (here, MyDiffEqn to compute $\dot{z}_{n}$ and advances the solution to $t_{n+1}$ by computing $z_{n+1}$ . Here is the MATLAB code:
```
function [zn1, tn1] = EulerForArray(zn, tn, h);
% EulerForArray computes the soltution zn1 at time
% tn1 of a set of differential equations  given the values
% of zn at time tn where tn1 = tn+h. The time step h must
% be specified by the user

zdot = MyDiffEqn(zn,tn);
zn1 = zn + h*zdot;
tn1 = tn + h;
end
```

•

Putting it all together: We now have these two functions. EulerForArray already calls the function MyDiffEqn. We only need to write a master code that sets up the entire solution by specifying initial conditions $z(0)$ , time step $h$ , and time span $t_{\rm span}$ as we did in the previous sample.

% ODESolver: Script file that solves the differential equation
% coded in MyDiffEqn.m using Euler’s method coded
% in EulerForArray.m with specified initial conditions
% and over the specified time interval

t0 = 0; tf = 10;         %specify time span [t0, tf]
z0 = [1; 0];             %specify initial conditions
h = 0.01;                %specify time step
nsteps = (tf-t0)/h;      %calculate number of time steps

t = t0; z = z0’;         %initialize output arrays t and z
tn = t; zn = z0;         %initialize tn and zn
for k=1:nsteps           %start a loop for computing z
    [zn1,tn1] = EulerForArray(zn,tn,h); %compute zn1 and tn1
            z = [z; zn1’];              %append zn1 in z
            t = [t; tn1];               %append tn1 in t
            zn = zn1; tn = tn1;         % reset zn and tn
end
x = z(:,1);              %unpack x from z (1st column)
y = z(:,2);              %unpack y from z (2nd column)
subplot(2,1,1), plot(t,x), xlabel(’t’), ylabel(’x’)
subplot(2,1,2), plot(t,y), xlabel(’t’), ylabel(’y’)

Filename:sfig8-6-2

Figure 3.17: Plot of

x

t

and

y

t

as obtained from the numerical solution

The two plots obtained from executing the script file are shown in fig. 3.17. The given differential equation represents a harmonic oscillator $\ddot{x}+\alpha x=0$ and hence the solution for $x$ and $\dot{x}(=y)$ are sinusoidal functions as expected.

Comments:

•

The most important thing to realize here is that $z$ is an array of independent variables and $\dot{z}$ is an array of first derivative of the variables in $z$ . If we had a set of $n$ first order equations, we could pack them all in the same $\dot{z}$ . Having all equations in an array is very convenient for solving them using numerical solution algorithms. Note that Euler’s method works the same way on an array as it does on a single equation.
•
The codes presented here are written for clarity, not for minimizing code length or maximizing efficiency. We could, for example, condense all lines of code in MyDiffEqns to just two lines:
```
alpha = 1;
zdot = [z(2); -alpha*z(1)];
```
Similarly, we could reduce the number of lines of code in the final script file ODESolver by combining different lines. The code then, admittedly, becomes a bit cryptic.
•

Specifying direct numerical values to variables inside a function is usually not a good idea (it is called hard coding). Such values should be passed through input variables. For example, we specified the value of alpha in MyDiffEqns and values of time step h in EulerForArray. It would be much nicer to pass their values as input to the function so that if we want to solve the differential equations with various values of such parameters, we do not have to edit these functions every time.

SAMPLE 3.9 Using canned programs for solving ODEs: Find the solution of the following nonlinear pendulum equation over the time interval $0\leq t\leq 20$ using any numerical ODE solver:

\ddot{\theta}+\sin\theta=0;\quad\theta(0)=\frac{19\pi}{20},\quad\dot{\theta}(0% )=0.

Solution Most mathematical software packages offer canned programs for solving ODEs. Here we will use MATLAB’s most popular and generally robust ODE solver, ode45, which is an implementation of fourth and fifth order Runge-Kutta algorithm. Since we are going to use a canned program (MATLAB calls them functions), we will have to follow the syntax specified by the program for input and output variables. The built-in function ode45 requires the differential equations to be coded in the following format:

Filename:sfig8-6-2a

Figure 3.18: The ODEs need to be coded in a function with a specified syntax for the input and output so that the MATLAB’s built-in solvers can use them properly.

Here the function, myODEs, contains the ODEs of interest as a set of first order ODES cast in the form of an array where the output array zdot represents the time derivative of the input array z and the input instant t. In this sample, we show how to form and code these arrays.

However, before we proceed to write the function, we need to cast the given second order ODE as a set of first order ODEs. Following the discussion in the text, we do it by introducing another variable, say, $\alpha$ , where $\alpha=\dot{\theta}$ . Then, we can write the given pendulum equation as a set of first order equations in $\theta$ and $\alpha$ as follows:

	$\displaystyle\dot{\theta}$	$\displaystyle=$	$\displaystyle\alpha$
	$\displaystyle\dot{\alpha}$	$\displaystyle=$	$\displaystyle-\sin\theta.$

The initial conditions, now in terms of the new variables, are $\theta(0)=\frac{19\pi}{20}$ and $\alpha(0)=0$ .

For coding these ODEs, we need to put them in a convenient array form. Here the two first order ODEs are in variables $\theta$ and $\alpha$ . So, let $z_{1}\equiv\theta$ and $z_{2}\equiv\alpha$ , where $z_{1}$ and $z_{2}$ are the elements of input array z. Thus we have,

{\bf z}\equiv\left(\begin{array}[]{c}\dot{z}_{1}\\ \dot{z}_{2}\end{array}\right)=\left(\begin{array}[]{c}z_{2}\\ -\sin(z_{1})\end{array}\right)

and the initial conditions are

{\bf z(0)}\equiv\left(\begin{array}[]{c}{z}_{1}(0)\\ {z}_{2}(0)\end{array}\right)=\left(\begin{array}[]{c}\frac{19\pi}{20}\\ 0\end{array}\right).

Now we are ready to code the given differential equations in a function with the required input and output syntax:

   function [zdot] = myODE(t,z);
   % Call syntax:  [zdot] = myODE(t, z);
   %  Inputs are:  t = time
   %               z = [z(1); z(2)] = [theta; alpha]
   %   Output is:  zdot = [z(2); -sin(z(1))]

   % unpack z for clarity
   z1 = z(1);    z2 = z(2);
   % compute derivatives
   z1dot = z2;
   z2dot = -sin(z1);
   % now pack derivatives in the output array zdot
   zdot = [z1dot;  z2dot];

The next step is to use MATLAB’s built-in ODE solver, ode45, to get the solution of the given differential equation. We need to specify the time span and the initial conditions for ode45 to work. Time span is already given in the problem: $0\leq t\leq 20$ ; so tspan = [0 20], and the initial condition is z0 = [19*pi/20; 0].

Here is the script file containing the commands that run ode45 and show the solution in the form of a plot of $\theta$ vs $t$ over the specified time domain.

  tspan = [0 20];       % specify time span
  z0 = [19*pi/20; 0];   % specify initial conditions
  [t, z] = ode45(@myODE, tspan, z0);     % run ode45

  theta = z(:,1);       % unpack solution array z;
                        % theta is column 1 of z
  thetadot = z(:,2);    % thetadot is column 2 of z

  figure(1)             % open a figure window
  plot(t,theta)         % plot theta as a function of t
  xlabel(’Time (t)’), ylabel(’theta’)    % label axes

  figure(2)             % open another figure window
  plot(t,thetadot)      % plot thetadot as a function of t
  xlabel(’Time (t)’), ylabel(’thetadot’) % label axes

Filename:efig1-2-28

Figure 3.19: The solution obtained from MATLAB’s ODE solver is plotted here: (a) plot of

\theta

as a function of time

t

, and (b) plot of

\dot{\theta}

as a function of time

t

Problems for 3.2 ODEs, especially the common simple ODEs

3.2.1

Write as many different ODEs, and their general solutions, as you can.

3.3 More things to do with a computer, besides solving linear equations and ODEs

Mechanics is a physical subject. The concepts in mechanics do not depend on computers. But mechanics is also a quantitative subject; relevant amounts (of length, mass, force, moment, time, etc) are described with numbers, and relations are described using equations and formulas. Computers are very good with numbers and formulas. Thus, the modern practice of engineering mechanics uses computers. The most-needed computer skills for mechanics are:

•

solution of simultaneous linear algebraic equations,
•

plotting, and
•

numerical solution of ODEs (Ordinary Differential Equations).

More basically, an engineer also needs the ability to routinely evaluate standard functions ( $x^{3}$ , $\cos^{-1}\theta$ , etc.), to enter and manipulate lists and arrays of numbers, and to write short programs.

Classical languages, applied packages, and simulators

Programming in standard languages such as Fortran, Basic, C++, Java or Python probably takes too much time to use in solving simple mechanics problems. Thus an engineer needs to learn to use one or another widely available computational package (e.g., MATLAB, SCI-LAB, OCTAVE, MAPLE, MATHEMATICA, MATHCAD, LABVIEW, etc). We assume that students have learned, or are learning such a package. Although none of the homework here depends on such, we also encourage you to play with packaged mechanics simulators (e.g., INVENTOR, ADAMS, DADS, ODE, ALGODOO, etc) for testing and building your intuition.

How we explain computation

Solving a mechanics problem involves

1.

Reducing a physical problem to a well posed mathematical problem;
2.

Solving the math problem using some combination of pencil and paper and numerical computation; and
3.

Giving physical interpretation of the mathematical solution.

This book is primarily about setup (a) and interpretation (c), neither of which particularly depends on what method is used to solve the equations. If a problem requires computation, the exact computer commands vary from package to package. Because we don’t know which computer package you are using we show computer calculations using an informal pseudo-computer language. For reference, typical commands are summarized on page 3.3.

Required computer skills

Here, in a little more detail, are the primary computer skills you need.

•
Linear algebraic equations. Many mechanics problems are statics or ‘instantaneous mechanics’ problems. These problems involve trying to find some forces or accelerations at a given configuration of a system. These problems can generally be reduced to the solution of linear algebraic equations of this general type: solve

$\begin{array}[]{rlcrlcl}3&x&+&4&y&=&8\\ -7&x&+&\sqrt{2}&y&=&3.5\end{array}$

for $x$ and $y$ . In practice the number of variables and equations can be quite large. Some computer packages will let you enter equations almost as written above. In our pseudo language we would write:
```
           eqset =  {  3*x +       4*y = 8
                      -7*x + sqrt(2)*y = 3.5  }
           solve eqset for x and y
```
Other packages may require you to set up your equations in matrix form

$\underbrace{\left[\begin{array}[]{cc}\phantom{-}3&4\\ -7&\sqrt{2}\\ \end{array}\right]}_{A}\underbrace{\left[\begin{array}[]{c}x\\ y\\ \end{array}\right]}_{z}=\underbrace{\left[\begin{array}[]{c}8\\ 3.5\\ \end{array}\right]}_{b}\hbox to36.135pt{\hfil or\hfil}Az=b$

which in computer-speak might look something like this:
```
                  A =  [ 3      4
                        -7  sqrt(2) ]
                  b =  [ 8 3.5 ]’
                  solve  A*z=b for z
```
where A is a $2\times 2$ matrix, b is a column of 2 numbers (the ’ indicates that the row of numbers b should be transposed into a column), and the two elements of z are $x$ and $y$ . For systems of two equations, like above, a computer is hardly needed. But for systems of three equations, pencil and paper work is sometimes error-prone. Given the tedium, the propensity for error, and the availability of electronic alternatives, pencil and paper solution of four or more equations is an anachronism.
•

Plotting. In order to see how a result depends on a parameter, or to see how a quantity varies with position or time, it is useful to see a plot. Any plot based on more than a few data points or a complex formula is far more easily drawn using a computer than by hand. Most often you can organize your data into a set of $(x,y)$ pairs stored in an x list and a corresponding y list. A simple computer command will then plot $x$ vs $y$ . The pseudo-code below, for example, plots a circle using 100 points
```
        npoints = [0 1 2 3 ... 100]
        theta   = npoints * 2 * pi / 100
        x       = cos(theta)
        y       = sin(theta)
        plot y vs x
```
where npoints is the list of numbers from 0 to 100, theta is a list of 100 numbers evenly spaced between 0 and $2\pi$ , and x and y are lists of 100 corresponding $x, y$ coordinate points on a circle.
•
ODEs The result of using the laws of dynamics is often a set of ordinary differential equations which need to be solved. A simple example would be:

Find $x$ at $t=5$ given that $\displaystyle\frac{dx}{dt}=x$ and that at $t=0$ , $x=1$ .

The solution to this problem can be found easily enough by hand to be $x(5)=e^{5}$ . But often the differential equations are just too hard for pencil and paper solution. Fortunately the numerical solution of Ordinary Differential Equations (ODEs) is already programmed into scientific and engineering computer packages. The simple problem above is solved with computer code equivalent to these informal commands:
```
       ODES =  { xdot  = x }
       ICS  =  { xzero = 1 }
       solve ODES with ICS until t=5
```
which will yield a list of values for paired values for $t$ and $x$ the last of which will be $t=5$ and $x$ close to $e^{5}\approx 148.4$ .

Examples of informal computer commands

We use informal computer commands that are not as strict as any real computer package. You will translate informal commands, like those below, into commands your package understands. This reference table uses mathematical ideas which you may or may not know before you read this book, but these will be introduced in the text when needed.

x=7 Set the variable $x$ to 7.

omega=13 Set $\omega$ to 13.

u=[1 0 -1 0]
v=[2 3 4 pi] Define $u$ and $v$ to be the lists shown.

t= [.1 .2 .3 ... 5] Set $t$ to the list of 50 numbers implied by the expression.

y=v(3) sets $y$ to the third value of $v$ (in this case 4).

A=[1 2 3 6.9
5 0 1 12 ] Set $A$ to the array shown.

z= A(2,3) Set z to the element of $A$ in the second row and third column.

w=[3
    4
    2
    5] Define $w$ to be a column vector.

w = [3 4 2 5]’ Same as above. ’ means transpose.

u+v Vector addition. In this case the result is $[3\ 3\ 3\ \pi]$ .

u*v Element by element multiplication, in this case $[2\ 0\ -4\ 0]$ .

sum(w) Add the elements of $w$ , in this case 14.

cos(w) Make a new list, each element of which is the cosine of the corresponding element of $[w]$ .

mag(u) The square root of the sum of the squares of the elements in $[u]$ , in this case 1.41421…

u dot v The vector dot product of component lists $[u]$ and $[v]$ , (we could also write sum(u*v).

C cross D The vector cross product of ${\accentset{\rightharpoonup}{\boldsymbol{C}}}$ and ${\accentset{\rightharpoonup}{\boldsymbol{D}}}$ , assuming the three element component lists for $[C]$ and $[D]$ have been defined.

A matmult w Use the rules of matrix multiplication to multiply $[A]$ and $[w]$ .

eqset = {3x + 2y = 6
6x + 7y = 8} Define ‘eqset’ to stand for the set of 2 equations in braces.

solve eqset
for x and y Solve the equations in ‘eqset’ for $x$ and $y$ .

solve Ax=b for x Solve the matrix equation $[A][x]=[b]$ for the list of numbers $x$ . This assumes $A$ and $b$ have already been defined.

for i = 1 to N
such and such
end Execute the commands ‘such and such’ $N$ times, the first time with $i=1$ , the second with $i=2$ , etc

plot y vs x Assuming $x$ and $y$ are two lists of numbers of the same length, plot the $y$ values vs the $x$ values.

solve ODEs
with ICs
until t=5 Assuming a set of ODEs and ICs have been defined, use numerical integration to solve them and evaluate the result at $t=5$ .

With an informality consistent with what is written above, other commands will be introduced as needed.

SAMPLE 3.10

Backslash is a great command: x = A $\backslash$ b. Incredible!! :). Learn it. Use it. Your abilities are expanded. OK, just kidding. We need some simple samples here of algebra, linear equations, and plotting.

	$\displaystyle 1\,\mathrm{m}$	$\displaystyle=$	$\displaystyle C_{1}$
	$\displaystyle\hbox to36.135pt{and\hfil}2\,\mathrm{m}/\mathrm{s}$	$\displaystyle=$	$\displaystyle C_{2}.$

Chapter 3 Computation: Algebra, ODEs, etc.

Download this chapter by section

3.1 Basic summary of ODEs for mechanics

Constant position: u˙=0⇒u=C1

Constant velocity: u˙=A⇒u=A⁢t+C1

Zero acceleration: u¨=0⇒u=C1+C2⁢t

Constant acceleration: u¨=A⇒ ⁢u=A⁢t2/2+C1⁢t+C2

Exponential growth: u˙=λ⁢u⇒u=C1⁢eλ⁢t

Exponential decay: u˙=−λ⁢u⇒u=C1⁢e−λ⁢t

Acceleration proportional to position: u¨=λ2⁢u ⇒u=C1⁢eλ⁢t+C2⁢e−λ⁢t⁢ or⇒u=C3⁢cosh⁡(λ⁢t)+C4⁢sinh⁡(λ⁢t)

Harmonic oscillator: u¨=−λ2⁢u⁢ or ⁢u¨+λ2⁢u=0 ⇒ ⁢u=C1⁢sin⁡(λ⁢t)+C2⁢cos⁡(λ⁢t).

More ODEs

Learn, memorize, know, whatever

Problems for 3.1 Basic computation: linear equations and plotting

3.2 Numerical solution of ODEs

Problems for 3.2 ODEs, especially the common simple ODEs

3.3 More things to do with a computer, besides solving linear equations and ODEs

Classical languages, applied packages, and simulators

How we explain computation

Required computer skills

Constant position: ${\dot{u}=0}\ \Rightarrow u=C_{1}$

Constant velocity: ${\dot{u}=A}\ \Rightarrow\ {u=At+C_{1}}$

Zero acceleration: ${\ddot{u}=0}\ \Rightarrow\ u=C_{1}+C_{2}t$

Constant acceleration: ${\ddot{u}=A}\indent\Rightarrow\ \mbox{ }u=At^{2}/2+C_{1}t+C_{2}$

Exponential growth: ${\dot{u}=\lambda u}\ \Rightarrow\ {u=C_{1}e^{\lambda t}}$

Exponential decay: ${\dot{u}=-\lambda u}\ \Rightarrow\ {u=C_{1}e^{-\lambda t}}$

Acceleration proportional to position: ${\ddot{u}=\lambda^{2}u}$ $\Rightarrow{u=C_{1}e^{\lambda\ t}\ +\ C_{2}e^{-\lambda\ t}}\hfill\break\hbox{}% \hskip 28.90755pt\textrm{or}\Rightarrow{u=C_{3}\cosh(\lambda\ t)\ +\ C_{4}% \sinh(\lambda\ t)}$

Harmonic oscillator: $\ddot{u}=-\lambda^{2}u\mbox{~{}~{}~{}or~{}~{}~{}}\ddot{u}+\lambda^{2}u=0$ $\Rightarrow\ \mbox{ }u=C_{1}\sin(\lambda\ t)+C_{2}\cos(\lambda\ t)$ .