## Ordinary Differential Equations: Introduction and Examples

Ordinary Differential Equations (ODEs) are equations that relate the derivatives of one or more smooth functions with an independent variable . The term “ordinary” indicates that there is only one independent variable appearing in the equations.

### Examples of Ordinary Differential Equations

ODEs appear naturally in almost all engineering applications. Here are some examples:

#### Newton’s Second Law of Motion

One of the most ubiquitously used ordinary differential equations is Newton’s second law of motion, which relates the second derivative of the position of a particle (i.e., the acceleration) to the applied force on the particle. The relationship has the form:

where is the mass of the particle. The time is the independent variable, while the dependent variable is the position . Newton’s equation adopts different forms depending on the application. For example, to describe the free vibration of a spring-mass system with spring constant (Figure 1), then the equation has the form:

where is the position of the mass at the tip of the spring from the spring’s static equilibrium position.

To clarify, assume a spring-mass system hung vertically. The forces acting on the mass are the spring force and the weight of the mass. The weight of the mass equals equals where is the gravitational acceleration of the Earth. If denotes the (arbitrary) position of the mass measured from the unstretched length, , of the spring (Figure 2), the force by the spring equals and the Newton’s second law of motion for the hung mass can be written as,

The negative sign is due to the fact that the weight acts downward and the spring force acts upward when the spring is stretched as shown by the free body diagram in Figure 2.

At the static equilibrium (no vibration), the weight of the mass pulls the spring vertically down and initially gives the spring a constant stretch of . Therefore, can be written as , where is the position of the mass with respect to the static equilibrium of the mass. Therefore, the above equation can be written as,

which leads to,

The explicit solution to this equation has the following form:

where and the initial position and initial velocity of the mass, and is the natural frequency of the system (denoted as or ).

The displacement response of the system has the form shown in Figure 3.

If the external forces include a damping force (a force that is a function of the velocity and that is always opposite to the direction of motion), then, the equation has the form:

where is the damping coefficient.

The explicit solution to this equation has the following form:

where is the damping ratio and is the natural frequency of the system. Again and are the initial displacement and initial velocity of the mass, respectively.

The displacement response for the free vibration of a damped spring-mass system is shown in Figure 4.

The following code utilizes the built-in DSolve function in Mathematica to solve the above simple equations. Also, the Manipulate built-in function is used to visualize the effect of varying and on the resulting vibration of the system.

View Mathematica CodeClear[x, g, k, m] (*General Solution*) a = DSolve[m*x''[t] == - k*x[t] - beta*x'[t], x, t]; Print["General Solution (damped system)="] x = x[t] /. a[[1]] (*Solution*) Clear[x, g, k, m, beta] a = DSolve[{m*x''[t] == - k*x[t] - beta*x'[t], x[0] == 0, x'[0] == 0}, x, t]; Print["Solution with initial velocity and displacement =0"] x = x[t] /. a[[1]]; x = Simplify[x] Manipulate[g = 1;Graphics3D[{Cuboid[{-1, -1, -5}, {0, 0, -((E^(-(((beta + Sqrt[beta^2 - 4 k m]) t)/(2 m)))m (beta (-1 + E^((Sqrt[beta^2 - 4 k m] t)/m)) + (1 + E^((Sqrt[beta^2 - 4 k m] t)/m)-2 E^(((beta + Sqrt[beta^2 - 4 k m]) t)/(2 m))) Sqrt[beta^2 - 4 k m]))/(2 k Sqrt[beta^2 - 4 k m]))}]}, Axes -> True, PlotRange -> {{-1, 0}, {-1, 0}, {-5, 10}}, AxesOrigin -> {0, 0, 0}, ViewVertical -> {0, 0, -1}], {t, 0, 2*4*Pi*Sqrt[5]}, {m, 1, 5}, {beta, 0, 5}, {k, 1, 5}]

# UPGRADE: need Sympy 1.2 or later, upgrade by running: "!pip install sympy --upgrade" in a code cell # !pip install sympy --upgrade import cmath import numpy as np import sympy as sp import matplotlib.pyplot as plt from IPython.display import HTML from ipywidgets.widgets import interact, Play sp.init_printing(use_latex=True) display(HTML('<link rel="stylesheet" href="//stackpath.bootstrapcdn.com/font-awesome/4.7.0/css/font-awesome.min.css"/>')) x = sp.Function('x') m, k, t, beta = sp.symbols('m k t beta') # General Solution sol = sp.dsolve(m*x(t).diff(t,2) + k*x(t) + beta*x(t).diff(t)) display("General Solution (damped system):", sol) # Solution sol = sp.dsolve(m*x(t).diff(t,2) + k*x(t) + beta*x(t).diff(t), ics={x(0): 0, x(t).diff(t).subs(t, 0): 0}) display("Solution with initial velocity and displacement = 0:", sp.simplify(sol)) @interact(t=Play(value=0, min=0, max=60, step=1, interval=500), m=(1,5,0.01), k=(1,5,0.01), beta=(0,5,0.01)) def update(t, m=1, k=1, beta=0): g = 1 s = cmath.sqrt(beta**2 - 4*k*m) x = -((np.exp(-(((beta + s)*t)/(2*m)))*m*(beta*(-1 + np.exp((s*t)/m)) + (1 + np.exp((s*t)/m)-2*np.exp(((beta + s)*t)/(2*m)))*s))/(2*k*s)) plt.xlim(0,10); plt.ylim(-0.1,0.1) plt.yticks([]) plt.title('x = '+str(np.round(x,3))) plt.plot([0,x],[0,0],lw=50)

View Python Code

# Spring Animation # UPGRADE: need Sympy 1.2 or later, upgrade by running: "!pip install sympy --upgrade" in a code cell # !pip install sympy --upgrade import numpy as np import sympy as sp import matplotlib.pyplot as plt from IPython.display import HTML from scipy.integrate import odeint from matplotlib.patches import Rectangle from ipywidgets.widgets import interact display(HTML('<link rel="stylesheet" href="//stackpath.bootstrapcdn.com/font-awesome/4.7.0/css/font-awesome.min.css"/>')) r = sp.Function('r') t = sp.symbols('t') x_val = np.arange(0,5,0.01) @interact(displacement=(-0.1,0.1,0.05), velocity=(-1,1,0.05), mass=(5,10,0.05), damping_ratio=(0,1,0.01), stiffness=(200,750,1)) def update(displacement=0, velocity=0, mass=5, damping_ratio=0, stiffness=200): run = 1 cd = 1.5*(2*np.sqrt(stiffness*mass)); w = np.sqrt(stiffness/(mass)) w1 = w/(2*np.pi) Tn = 1/w1 U = (displacement + velocity)*Tn nt = run/Tn da = damping_ratio*cd sol = sp.dsolve(mass*r(t).diff(t,2) + da*r(t).diff(t) + stiffness*r(t), ics={r(0): displacement, r(t).diff(t).subs(t,0): velocity}) l = sol.rhs.subs(t,run) plt.xlim(-1, 20) plt.ylim(-6, 6) plt.plot([-0.5,0,0],[0,0,2], 'k') plt.plot([3 + l, 3 + l, 4.5 + l],[-2, 0, 0], 'k') plt.plot([3*nx/11 + l*nx/11 for nx in range(12)], [2*np.cos(np.pi*nx) for nx in range(12)], 'k') ax = plt.gca() ax.add_patch(Rectangle((-1,-6),0.5,12,color='black')) ax.add_patch(Rectangle((4.5 + l, -2),4,4,color='rosybrown')) plt.text(4,3,"f_n(Hz) = {}".format(round(w1,5))) plt.text(4,4,"T_n(s) = {}".format(round(Tn,5))) plt.show() plt.plot(x_val, [sol.rhs.subs(t,i) for i in x_val]) plt.grid(); plt.show()

#### Beam Deflection

The deflection of an Euler-Bernoulli beam follows the following equation:

where is the beam’s Young’s modulus, is the cross sectional moment of inertia, is the beam’s deflection (the dependent variable), is the applied distributed load on the beam, and is the position along the beam’s length (the independent variable). For a fixed-ends beam; if the deflection and rotation at both ends of a beam of length is equal to zero, i.e. , and if is constant, then, the deflection of the beam has the form:

#### Growth and Decay

If is the number of organisms at any time , and if the growth of the number depends on the number that is already there, then, the ODE that would describe this uninhibited growth is given by:

Where is a positive constant. The time is the independent variable, while the dependent variable is the number of organisms . The general solution of this equation has the form:

where is a constant depending on the boundary conditions. If the initial condition is such that then, the solution has the form:

Similarly, if rate of decrease (decay) of a quantity is a function of the amount available, then, the uninhibited decay can be described by the equation:

where is a positive constant. The time is the independent variable, while the dependent variable is the quantity . Radioactive decay follows this ODE. The general solution of this equation has the form:

where is a constant depending on the boundary conditions. If the initial condition is such that then, the solution has the form:

The following tool plots the quantities or when for different values of the exponent:

View Mathematica CodeClear[x] a = DSolve[x'[t] == k*x[t], x, t]; x = x[t] /. a[[1]] Clear[x] a = DSolve[{x'[t] == k*x[t], x[0] == A}, x, t] x = x[t] /. a[[1]] Manipulate[A = 1; r = Grid[{{"k=", k}}, Spacings -> 0]; Grid[{{Plot[A*E^(k*t), {t, 0, 20}, AxesLabel -> {"t", "x(t) or y(t)"}, ImageSize -> Medium]}, {r}}], {k, -0.5, 0.5}]

# UPGRADE: need Sympy 1.2 or later, upgrade by running: "!pip install sympy --upgrade" in a code cell # !pip install sympy --upgrade import numpy as np import sympy as sp import matplotlib.pyplot as plt from ipywidgets.widgets import interact sp.init_printing(use_latex=True) x = sp.Function('x') t, k, A = sp.symbols('t k A') sol = sp.dsolve(x(t).diff(t) - k*x(t)) display(sol) sol = sp.dsolve(x(t).diff(t) - k*x(t), ics={x(0): A}) display(sol) @interact(k=(-0.5,0.5,0.01)) def update(k=-0.5): x = np.arange(0,20,0.1) y = np.exp(k*x) plt.plot(x,y) plt.xlabel('t'); plt.ylabel('x(t) or y(t)') plt.grid(); plt.show()

#### Electrical Capacitors and Inductors

Simple electric circuits are composed of three components: resistances, capacitors, and inductors. The relationship between the voltage drop across a resistance, and the current passing through it is given by:

where is the constant of proportionality, namely the value of the electric resistance. A capacitor on the other hand is an element that has the following relationship between the current passing through it and the rate of change of the voltage drop across it :

where is the constant of proportionality, namely the capacitance. An inductor is an element that has the following relationship between the voltage drop across it and the rate of change of the current that passes through it :

Where is the constant of proportionality, namely the inductance. When the elements are combined, ODEs are generated with being the independent variable, and and being the dependent variables.