In this article, the simplest numeric method, namely the Euler method to solve 1st order ODEs will be demonstrated with examples (and then use it to solve 2nd order ODEs). Also, we shall see how to plot the phase lines (gradient fields) for an ODE and understand from examples how to qualitatively find a solution curve with the phaselines. All the problems are taken from the **edx Course: MITx — 18.03x: Introduction to Differential Equations**.

Let’s implement the above algorithm steps with the following R code snippet.

library(tidyverse) library(latex2exp) library(gridExtra) Euler <- function(f, xn, yn, h, N, n=0) { df <- NULL An <- f(xn, yn) while (n < N) { df <- rbind(df, data.frame(n=n, Xn=xn, Yn=yn, An=An, hAn=h*An)) An <- f(xn, yn) xn <- xn + h yn <- yn + h*An n <- n + 1 } return(df) }

Let’s define the function *f* next.

f <- function(x, y) { return(y*sin(x)) } df <- Euler(f, -3, 1, 0.1, 100) head(df)

The below table shows the solution curve computed by the Euler method for the first few iterations:

Let’s implement the following R functions to draw the isoclines and the phase lines for an ODE.

plot_isoclines <- function(f, rx, ry, rm, l, xgap, title) { lx <- rx[1] ux <- rx[2] x <- seq(lx, ux, xgap) idf <- NULL for (m in rm) { y <- f(x, m) dx <- l / 2 / sqrt(1 + m^2) idf %>% ggplot() + geom_segment(aes(x = x0, y = y0, xend = x1, yend = y1, col = m)) + xlim(rx[1], rx[2]) + ylim(ry[1], ry[2]) + scale_color_gradient2(low = "blue", mid = "green", high = "red") + ggtitle(title) + xlab('x') + ylab('y') + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) }

plot_phaselines <- function(f, rx, ry, dx, dy, l, ms=-1:1, title=NULL) { lx <- rx[1] ux <- rx[2] ly <- ry[1] uy <- ry[2] x <- seq(lx, ux, dx) y <- seq(ly, uy, dy) pdf <- expand.grid(x = x, y = y) pdf$m <- f(pdf$x, pdf$y) pdf$dx <- l / 2 / sqrt(1 + pdf$m^2) pdf$x0 <- pdf$x - pdf$dx pdf$x1 <- pdf$x + pdf$dx pdf$y0 <- pdf$y - pdf$m*pdf$dx pdf$y1 <- pdf$y + pdf$m*pdf$dx pdf %>% ggplot() + geom_segment(aes(x = x0, y = y0, xend = x1, yend = y1, col = m)) + xlim(rx[1], rx[2]) + ylim(ry[1], ry[2]) + ggtitle(title) + xlab('x') + ylab('y') + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) }

Now, let’s plot the 0 (the nullcline), 1, 2, 3 isoclines for the ODE: *y’ = y.sin(x)*.

ifunc <- function(x, m) m / sin(x) rx <- c(-4,4) ry <- c(-3,3) rm <- seq(-3,3,0.5) l <- 0.2 xgap <- 0.05 ygap <- 0.05 plot_isoclines(ifunc, rx, ry, rm, l, xgap, TeX('m-isoclines for $\\frac{dy}{dx} = y.sin(x)$'))

xgap <- 0.25 #05 ygap <- 0.25 #05 plot_phaselines(f, rx, ry, xgap, ygap, l, -1:1, TeX('phaselines and m-isoclines for $\\frac{dy}{dx} = y.sin(x)$'))

Now, let’s solve the ODE with Euler’s method, using step size *h = 0.1* with the initial value *y(-3) = 1*.

l <- 0.2 xgap <- 0.05 #05 ygap <- 0.05 #05 p <- plot_phaselines(f, rx, ry, xgap, ygap, l, -1:1, TeX('phaselines and m-isoclines for $\\frac{dy}{dx} = y.sin(x)$')) df <- Euler(f, -3, 1, 0.1, 100) p + geom_point(aes(Xn, Yn), data=df) + geom_path(aes(Xn, Yn), arrow = arrow(), data=df)

The black curve in the following figure shows the solution path with Euler.

Let’s plot the phase lines for the ODE *y’ = y ^{2} – x*, solve it with Euler and plot the solution curve again, using the following code.

f <- function(x, y) { return(y^2-x) } rx <- c(-4,4) ry <- c(-3,3) rm <- seq(-3,3,0.5) l <- 0.2 xgap <- 0.1 ygap <- 0.1 p <- plot_phaselines(f, rx, ry, xgap, ygap, l, -2:2, TeX('phaselines and m-isoclines for $\\frac{dy}{dx} = y^2-x$')) df <- Euler(f, -1, 0, 0.1, 100) p + geom_point(aes(Xn, Yn), data=df) + geom_path(aes(Xn, Yn), arrow = arrow(), data=df)

The black solution curve shows the trajectory of the solution obtained with Euler.

Finally, let’ plot the isoclines for the ODE *y’ = y ^{2} – x^{2}* and solve it with Euler. Let’s try starting with different initial values and observe the solution curves in each case.

f <- function(x, y) { return(y^2-x^2) } rx <- c(-4,4) ry <- c(-4,4) l <- 0.2 xgap <- 0.1 ygap <- 0.1 p <- plot_phaselines(f, rx, ry, xgap, ygap, l, -2:2, TeX('phaselines and m-isoclines for $\\frac{dy}{dx} = y^2-x^2$, solving with Euler with diff. initial values')) df <- Euler(f, -3.25, 3, 0.1, 100) p <- p + geom_point(aes(Xn, Yn), data=df) + geom_path(aes(Xn, Yn), arrow = arrow(), data=df) df <- Euler(f, -1, 1, 0.1, 50)

As can be seen, the solution curve starting with different initial values lead to entirely different solution paths, it depends on the phase lines and isoclines it gets trapped into.

Finally, let’s notice the impact of the step size on the solution obtained by the Euler method, in terms of the error and the number of iterations taken to obtain the solution curve. As can be seen from the following animation, with bigger step size h=0.5, the algorithm converges to the solution path faster, but with much higher error values, particularly for the initial iterations. On the contrary, smaller step size h=0.1 taken many more iterations to converge, but the error rate is much smaller, the solution curve is pretty smooth and follows the isoclines without any zigzag movement.

**Spring-mass-dashpot system** – the Damped Harmonic Oscillator

A spring is attached to a wall on its left and a cart on its right. A dashpot is attached to the cart on its left and the wall on its right. The cart has wheels and can roll either left or right.

The spring constant is labeled *k*, the mass is labeled *m*, and the damping coefficient of the dashpot is labeled b. The center of the cart is at the rest at the location *x* equals 0, and we take rightwards to be the positive *x* direction.

Consider the spring-mass-dashpot system in which mass is* m *(kg), and spring constant is *k* (Newton per meter). The position *x(t)* (meters) of the mass at (seconds) is governed by the DE

where m, b, k > 0. Let’s investigate the effect of the damping constant, on the solutions to the DE.

The next animation shows how the solution changes for different values of 𝑏, given 𝑚=1 and 𝑘=16, with couple of different initial values (c).

Now, let’s observe how the amplitude of the system response increases with the increasing frequency of the input sinusoidal force, achieving the peak at resonance frequency (for *m=1*, *k=16* and *b=√14*, the resonance occurs at *w _{r} = 3*)

## Solving the 2nd order ODE for the damped Harmonic Oscillator with Euler method

Now, let’s see how the 2nd order ODE corresponding to the damped Harmonic oscillator system can be numerically solved using the Euler method. In order to be able to use Euler, we need to represent the 2nd order ODE as a collection of 1st order ODEs and then solve those 1st order ODEs simultaneously using Euler, as shown in the following figure:

Let’s implement a numerical solver with Euler for a 2nd order ODE using the following R code snippet:

Euler2 <- function(f, tn, xn, un, h, N, n=0) { df <- NULL while (n < N) { df <- rbind(df, data.frame(n=n, tn=tn, xn=xn, un=un)) tn <- tn + h xn <- xn + h*un un <- un + h*f(tn, xn, un) n <- n + 1 } return(df) }

For the damped Harmonic Oscillator, we have the 2nd order ODE

from which we get the following two 1st order ODEs:

Let’s solve the damped harmonic oscillator with Euler method:

f <- function(t, x, u) { return(-b*u/m - k*x/m) } m <- 1 b <- 1 k <- 16 df <- Euler2(f, 0, 1, 0, 0.1, 100) head(df)

Now, let’s plot the solution curves using the following code block:

par(mfrow=c(1,2)) plot(df$tn, df$xn, main="x(t) vs. t", xlab='t', ylab='x(t)', col='red', pch=19, cex=0.7, xaxs="i", yaxs="i") lines(df$tn, df$xn, col='red') grid(10, 10) plot(df$tn, df$un, main=TeX('$\\dot{x}(t)$ vs. t'), xlab='t', ylab=TeX('$\\dot{x}(t)$'), col='blue', pch=19, cex=0.7, xaxs="i", yaxs="i") lines(df$tn, df$un, col='blue') grid(10, 10)

The following animation shows the solution curves obtained by joining the points computed with the Euler method for the damped Harmonic Oscillator: