Preface

Euler's method is first order method. It is a straight-forward method that estimates the next point based on the rate of change at the current point and it is easy to code. It is a single step method. Notably, Forward Euler's method is unconditionally unstable for un-damped oscillating systems (such as a spring-mass system or wave equations) in space desctretization. For complex problems and/or boundary condition it may fail. It can be used for basic numerical analysis. This method is not commonly used for spatial discretization but some times used in time discretization. This scheme is not recommended for hyperbolic differential equation because this is more diffusive. Order of convergence of this scheme with grid refinement is very poor. Extending Euler method to higher order method is easy and straight forward.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the first course APMA0330
Return to the main page for the second course APMA0340
Return to Part IV of the course APMA0340
Introduction to Linear Algebra with Mathematica

Glossary

Euler Methods

A common example of a physics problem that requires the solution of a differential equation is the motion of a particle acted on by a force. For simplicity, we first discuss one-dimensional motion so that only a single vector component of position, velocity, and acceleration are needed. These variables are related using the language of differential calculus

\[ \frac{{\text d}y(t)}{{\text d}t} = v(t) , \qquad \frac{{\text d}v(t)}{{\text d}t} = a(t) , \]
where y is the displacement of a particle, v is its velocity, and a is the acceleration. These quantities are known as kinematical quantities, because they describe the motion without regard to its cause.

Why do we need the concept of acceleration in kinematics? The answer can be found only aposteriori. Thanks to Newton, we know that the net force acting on a particle determines its acceleration. Newton’s second law of motion tells us that

\[ F(y, v, t) = m\, a(t) . \]
where F is the net force and m is the inertial mass. In general, the force depends on position, velocity, and time. Note that Newton’s law implies that the motion of a particle does not depend on \( \ddot{v} = {\text d}^2 v / {\text d}t^2 \) or on any higher derivative of the velocity. It is a property of nature, not of mathematics, that we can find simple explanations for motion.

The two first-order equations can be combined to obtain the familiar second- order differential equation for the position:

\[ m\,\frac{{\text d}^2 y}{{\text d}t^2} = F(t,y,\dot{y}) . \]
However, we will find it easier to describe the motion of a particle by two coupled first-order differential equations \( \dot{y} =v, \ \dot{v} = a . \)

Because many types of systems can be modeled by the above system of differential equations such, it is important to know how to solve such equations. In general, analytical solutions of differential equations, that is, solutions in terms of well-known functions, do not exist. We are therefore motivated to find numerical solutions of differential equations. However, analytical solutions are very important and often exist in special or limiting cases. We often use them to test our numerical solutions.

The power of mathematics when applied to physics comes in part from the fact that frequently seemingly unrelated problems can have the same mathematical formulation. Hence, if we can solve one problem, we can solve other problems that might appear to be unrelated.

We start demonstrating the Euler methods to Newton’s equations of motion. We write Newton’s second as a system of coupled first-order differential equations and apply the Euler algorithm to each variable. The variables of interest are

\[ \begin{split} t_{k+1} &= t_k + \Delta t , \\ y_{k+1} &= y_k + v_k \Delta t , \\ v_{k+1} &= v_k + a_k \Delta t . \end{split} \]

In this section, we extend Euler methods to systems of differential equations. It is a custom to denote by t the independent variable that we associate with time. Then denoting the derivatives by dots, we illustrate the concepts by considering the initial value problem

\[ \begin{cases} \dot{x} &= f(t, x, y) , \\ \dot{y} &= g(t,x,y) , \end{cases} \qquad\mbox{subject to }\qquad \begin{cases} x(t_0 ) = x_0 , \\ y(t_0 ) = y_0 . \end{cases} \]
Here f and g are given functions and t0, x0, and y0 are specified points. Using vectors, we can rewrite the initial value problem in compact form:
\[ \dot{\bf x} = {\bf f} (t,{\bf x}) , \qquad {\bf x} (t_0 ) = {\bf x}_0 , \]
where
\[ {\bf x}(t) = \begin{bmatrix} x(t) \\ y(t) \end{bmatrix} , \qquad {\bf f}(t, {\bf x}) = \begin{bmatrix} f(t,x,y) \\ g(t,x,y) \end{bmatrix} , \qquad {\bf x}_0 = \begin{bmatrix} x_0 \\ y_0 \end{bmatrix} . \]
A solution to the above problem is a pair of differentiable functions x(t) and y(t) with the property that when x(t) and y(t) are substituted in f(t,x,y) and g(t,x,y), the result is equal to the derivative \( \dot{x} \) and \( \dot{y} , \) respectively.

A numerical solution to the given system over the interval \( a \le t \le b \) is found by considering the differentials

\[ {\text d}x = f(t,x,y)\,{\text d}t \qquad\mbox{and}\qquad {\text d}y = g(t,x,y)\,{\text d}t \]
Euler's method for solving the system is obtained by substituting approximations to the differential dt = tk+1 - tk,    dx = xk+1 - xk,   , and   dy = yk+1 - yk into definition of the differential to get
\[ \begin{split} x_{k+1} - x_k &= f(t_k , x_k , y_k ) , \\ y_{k+1} - y_k &= g(t_k , x_k , y_k ) . \end{split} \]
The interval is divided into m subintervals of width \( h = (b-a)/m , \) and the mesh points become \( t_{k+1} = t_k +h . \) This is used to get the recursive formulas for Euler's rule:
\[ \begin{split} t_{k+1} &= t_k +h , \\ x_{k+1} &= x_k + h \,f(t_k , x_k , y_k ) , \\ y_{k+1} &= y_k + h \,g(t_k , x_k , y_k ) , \qquad \mbox{for}\quad k=0,1,2,\ldots , m-1. \end{split} \]

Euler method for systems of differential equations

euler[F_, a_, Y0_, b_, Steps_] :=
Module[{t, Y, h = (b - a)/Steps//N, i},
t[0] = a; Y[0] = Y0;
Do[
t[i] = a + i h;
Y[i] = Y[i - 1] + h F[t[i - 1], Y[i - 1]],
{i, 1, n}
];
Table[{t[i], Y[i]}, {i, 0, n}]
]

Using for simplicity the initial point to be the origin. we present another way to implement the Euler rule:

h = .1;                (* step size *)
TimeMax = 3.5;
maxN = TimeMax/h;
plotmax = 5;
plotmin = -5;
fx[x_, y_] = y + x
Out[8]= x+y
fy[x_, y_] = -2 x + 3 y
Out[9]= -2 x + 3 y
xxo = -1.5;
yyo = -0.5;
xx[0] := xxo
yy[0] := yyo
xx[n_] := xx[n] = xx[n - 1] + h*fx[xx[n - 1], yy[n - 1]]
yy[n_] := yy[n] = yy[n - 1] + h*fy[xx[n - 1], yy[n - 1]]
MatrixForm[Table[{h*n, xx[n], yy[n]}, {n, 0, maxN}]]
Out[15]//MatrixForm=

We can make several different types of plots. Here we plot on the phase plane using dots and then using line segments

eulerdotplot =
ListPlot[Table[{xx[n], yy[n]}, {n, 0, maxN}],
PlotStyle -> {PointSize[Large], Red}]
eulerlineplot =
ListLinePlot[Table[{xx[n], yy[n]}, {n, 0, maxN}],
PlotStyle -> {Thick, Green}]

Here we create the associated streamplot

streamplot =
StreamPlot[{fx[x, y], fy[x, y]}, {x, plotmin, plotmax}, {y, plotmin,
plotmax}]

 

We plot x-component versus t

xPlot = ListLinePlot[Table[{h*n, xx[n]}, {n, 0, maxN}],
PlotStyle -> {Thick, Blue}]

We plot y-component versus t

yPlot = ListLinePlot[Table[{h*n, yy[n]}, {n, 0, maxN}],
PlotStyle -> {Thick, Green}]

Finally, we mash them all together

Show[streamplot, eulerlineplot]


 

We know from the previous section that Euler's rule is not recommended for practical usage. Recalling the backward Euler method, we arrive at the Euler--Cromer algorithm:

\[ \begin{split} t_{k+1} &= t_k + \Delta t , \\ v_{k+1} &= v_k + a_k \Delta t , \\ y_{k+1} &= y_k + v_{k+1} \Delta t . \end{split} \]
It might occur to you that it would be better to compute the velocity at the middle of the interval rather than at the beginning or at the end of the interval. The Euler--Richardson algorithm is based on this idea. This algorithm is particularly useful for velocity-dependent forces, but does as well as other simple algorithms for forces that do not depend on the velocity. The algorithm consists of using the Euler algorithm to find the intermediate position ymid and velocity vmid at a time tmid = t +∆t/2. Then we compute the force, F(tmid, ymid, vmid) and the acceleration amid at tmid. The new position yk+1 and velocity vk+1 at time tk+1 is found using vmid and amid. We summarize the Euler--Richardson algorithm as follows:
\[ \begin{split} a_k &= F(t_k , y_k , v_k )/m , \\ v_{\small mid} &= v_k + \frac{1}[2}\, a_k \Delta t , \\ y_{\small mid} &= y_k + \frac{1}[2}\, v_k \Delta t , \\ y_{\small mid} &= y_k + \frac{1}[2}\, v_k \Delta t , \\ a_{\small mid} &= F \left( t + \frac{1}[2}\, \Delta t ,y_k, v_k \right) \\ \end{split} \]
and
\[ \begin{split} t_{k+1} &= t_k + \Delta t , \\ v_{k+1} &= v_k + a_{\small mid} \Delta t , \\ y_{k+1} &= y_k + v_{\small mid} \Delta t \qquad\mbox{Euler--Richardson algorithm}. \end{split} \]
Although we need to do twice as many computations per time step, the Euler--Richardson algorithm is much faster than the Euler algorithm because we can make the time step larger and still obtain better accuracy than with either the Euler or Euler--Cromer algorithms.

Example: Consider the initial https://www.researchgate.net/publication/237621801_Motion_in_One_Dimension ■

 

Example: Consider the initial ■

 

Stanley, R.W., Numerical metods in mechnaics, American Journal of Physics, 1984, Vol. 52, pp. 499--507.
  1. William R. Bennett, Scientific and Engineering Problem-Solving with the Computer, Prentice Hall (1976)
  2. Byron L. Coulter and Carl G. Adler, “Can a body pass a body falling through the air?,” American Journal of Physics, 47, 841 (1979).
  3. A. P. French, Newtonian Mechanics, W. W. Norton & Company (1971). Chapter 7 has an excellent discussion of air resistance and a detailed analysis of motion in the presence of drag resistance
  4. Ian R. Gatland, “Numerical integration of Newton’s equations including velocity-dependent forces,” American Journal of Physics, 62, 259 (1994). The author discusses the Euler-Richardson algorithm
  5. Margaret Greenwood, Charles Hanna, and John Milton, “Air resistance acting on a sphere: numerical analysis, strobe photographs, and videotapes,” The Physics Teacher, 1986, Vol. 24, 153 (1986); https://doi.org/10.1119/1.2341969. More experimental data and theoretical analysis are given for the fall of ping-pong and styrofoam balls. Also see Mark Peastrel, Rosemary Lynch, and Angelo Armenti, “Terminal velocity of ashuttlecock in vertical fall,” American Journal of Physics, 48, 511 (1980).
  6. K. S. Krane, “The falling raindrop: variations on a theme of Newton,” Am. J. Phys. 49, 113 (1981). The author discusses the problem of mass accretion by a drop falling through a cloud of droplets
  7. Mehta, R.D., “Aerodynamics of sports balls,” Annual Review of Fluid Mechanics, 1985, Vol. 17, pp. 151--189. https://doi.org/10.1146/annurev.fl.17.010185.001055
  8. Neville de Mestre, The Mathematics of Projectiles in Sport, Cambridge University Press (1990). The emphasis of this text is on solving many problems in projectile motion, for example, baseball, basketball, and golf, in the context of mathematical modeling. Many references to the relevant literature are given.
  9. Forman S. Acton, Numerical Methods that work, Harper & Row (1970); corrected edition, Mathematical Association of America (1990). A somewhat advanced, but clearly written text in the 1960s; it is more about the 1950s world of computation; it makes only brief mention of the QR algorithm for eigenvalues and does not cover the singular value decomposition or variable step size ODE solvers. Moreover, the author has an aversion to library routines and to rigorous error bounds.

 

Return to Mathematica page

Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations
Return to the Part 7 Special Functions