# DE solution in R (nonlinear oscillator)

**Dan Kelley Blog/R**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

# Introduction

The function `lsoda()`

from the `deSolve`

package is a handy function for

solving differential equations in R. This is illustrated here with a classic

example: the nonlinear oscillator.

# Theory

As explained in any introductory Physics textbook, the nondimensionalized

oscillator equation can be simplified to

provided that , i.e. in the “small

angle” limit.

The linear form has solution for initial conditions

and at .

The nonlinear form is harder to solve analytically, but numerical integration

can be applied to any situation of interest. This is made simpler in the

present statement of the problem in nondimensional form, since there is but a

single parameter, , that describes the system.

A few questions may come to mind before proceeding further. For example, it is

clear that if , the numerical solution should match the solution to

the linear equation. But just how small must be, for a given precision?

A second question is about the qualitative form of the numerical solution, e.g

is it still a sinusoid but a different frequency, or something of a different

character? Might there different classes of solutions in various ranges of

?

Showing that such questions are easily answered with R is the point of the

present exercise.

# Framing the DE for solution in R

As an exercise, `lsoda()`

from the `deSolve`

package can be used to

integrate the nonlinear DE without the small angle assumption.

The first step is to break the second-order DE down into two first-order DEs,

and , which are to be defined

with a function named `func`

that is used by the DE solver named `lsoda`

.

Of course, we also need initial conditions and a set of times at which to

report the solution. All of these things are set out in the R code given

below, which plots the solution for as the red dashed line, on top of

the theoretical solution as the blue line. Readers might wish to try this with

slightly larger and smaller values of , to get a feel for the solution.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
library(deSolve) de <- function(t, y, parms, ...) { theta <- y[1] phi <- y[2] list(c(dthetadt = phi, dphidt = -sin(theta))) } a <- 0.1 y0 <- c(0, a) t <- seq(0, 4 * pi, pi/100) sol <- lsoda(y = y0, times = t, func = de) ylim <- max(range(sol[, 2])) * c(-1, 1) par(mar = c(3.5, 3.5, 1, 1), mgp = c(2, 0.7, 0)) plot(sol[, 1], sol[, 2], type = "l", ylim = ylim, col = "blue", xlab = expression(t), ylab = expression(theta(t))) grid() lines(t, a * sin(t), col = "red", lty = "dashed") |

# Some test cases

For more exploration, it is convenient to define the above as a stand-alone

function that takes `a`

as a parameter.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
library(deSolve) oscillator <- function(a = 0.1) { de <- function(t, y, parms, ...) { theta <- y[1] phi <- y[2] list(c(dthetadt = phi, dphidt = -sin(theta))) } y0 <- c(0, a) t <- seq(0, 8 * pi, pi/100) sol <- lsoda(y = y0, times = t, func = de) ylim <- max(range(sol[, 2])) * c(-1, 1) par(mar = c(3.5, 3.5, 1, 1), mgp = c(2, 0.7, 0)) plot(sol[, 1], sol[, 2], type = "l", ylim = ylim, col = "blue", xlab = expression(t), ylab = expression(theta(t))) grid() lines(t, a * sin(t), col = "red") legend("bottomleft", col = c("blue", "red"), lwd = 1, legend = c("linear", "nonlinear"), bg = "white") } |

Now, a few examples are easy to construct.

Start with a somewhat nonlinear problem

1 |
oscillator(1) |

Readers should try increasing a bit at a time, e.g. the example below has

a distinctly non-sinusoidal character.

1 |
oscillator(1.999) |

# Exercises

Further explore the behaviour in the neighborhood of . Are changes

subtle or dramatic in that region? Are there other regions of interest?

Consult the literature if this problem interests you.

# Resources

- Source code: 2014-06-15-nonlinear-oscillator.R

**leave a comment**for the author, please follow the link and comment on their blog:

**Dan Kelley Blog/R**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.