DE solution in R (nonlinear oscillator)

June 15, 2014
By

(This article was first published on Dan Kelley Blog/R, and kindly contributed to R-bloggers)

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")

center

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)

center

Readers should try increasing a bit at a time, e.g. the example below has
a distinctly non-sinusoidal character.

1
oscillator(1.999)

center

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

To 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 on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Sponsors

Mango solutions



RStudio homepage



Zero Inflated Models and Generalized Linear Mixed Models with R

Quantide: statistical consulting and training

datasociety

http://www.eoda.de





ODSC

ODSC

CRC R books series





Six Sigma Online Training









Contact us if you wish to help support R-bloggers, and place your banner here.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)