How do I solve a second order differential equation in R?

Question

I am learning R to solve a second order differential equation(probably using deSolve package). Which I have written in python by writing it as two first order differential equations and is given below

import  numpy  as np
import  matplotlib.pyplot  as plt
from  scipy.integrate  import  odeint

def  fun(X, t):
    y , dy , z = X
    M = np.sqrt (1./3. * (1/2. * dy **2 + 1./2.  * y **2))
    dz = (M*z) # dz/dt
    ddy =  -3.* M * dy -  y # ddy/dt

    return [dy ,ddy ,dz]

y0 = 1

dy0 =  -0.1

z0 = 1.

X0 = [y0, dy0, z0]

M0 = np.sqrt (1./3. * (1./2. *  dy0 **2. + 1./2.*  y0 **2)) 

t = np.linspace(0., 100., 10001.) # time spacing

sol = odeint(fun, X0, t)

y = sol[:, 0]

dy = sol[:, 1]

z = sol[:, 2]

M = np.sqrt (1./3. * (1./2. *  dy**2. + 1./2.*  y **2)) 

#Graph plotting
plt.figure()
plt.plot(t, y)
plt.plot(t, z)
plt.plot(t, M)
plt.grid()
plt.show()

Python solve this problem easily, however for another similar but complicated problem python shows error. I have also tried ode(vode/bdf) in python but the problem is still there. Now, I would like to check how R work with the problem. So, I will be much obliged if someone please provide me an example(basically a code translation!) of how to solve this problem in R so that I can try the other one in R and also learn some R(I know that this may not be the ideal way to learn a language).

I understand that this question may have little constructive value, but I am just a novice in R, so please bear with me!

[This paper](https://journal.r-project.org/archive/2010-2/RJournal_2010-2_Soetaert~et~al.pdf) gives a good introduction on how to use `deSolve` in `R` with a minimal example about half way through. However I'd suggest you think about why Python shows an error and why you think using `R` will solve your problem. — Michael Bird, Sep 06 '17 at 12:04
I am currently trying to understand why python is showing the error, if I don't find an answer I will ask for help in stackoverflow. However, searching for help, I came across this post: https://stackoverflow.com/questions/16973036/odd-scipy-ode-integration-error And tried to implement ode(vode/bdf) but the problem persists. So now trying R — Archimedes, Sep 06 '17 at 12:25
I think it will be hard for anyone to help you further without information about the error you have. The post you have linked suggests using `scipy.integrate.ode` to specify whether to use stiff or non-stiff methods. I guess you've done this? — Michael Bird, Sep 06 '17 at 12:32
Yes, I have implemented `scipy.integrate.ode` . For the stiff problem in my case, I first tried the `vode/bdf` and then eventually exhausted all the possible combinations and different orders! But, the problem is still there. — Archimedes, Sep 06 '17 at 12:45
I think @MichaelBird is right: translating this into R probably won't solve the problem. I'd venture that both R and Python probably ultimately rely on the same old school ODE code that has just been repackaged. — Dan, Sep 06 '17 at 12:48
Yeah, I would be very happy if both of them have the same problem. Then, I will try to reformulate the problem (In fact, I am trying to check the physics of the issue). But, at least I would like to check it with `R` as I have already done with `Mathematica` and `Matlab`. But, learning to code in `R` from scratch is now a problem for me, hence I asked for help. — Archimedes, Sep 06 '17 at 12:56
You're on the right track. If "this works but another similar but complicated problem python shows error", you have to show us the problem that *fails* (and, ideally, the error; the error alone *might* give someone a hint). Is it the same problem as in the linked question? — Ben Bolker, Sep 06 '17 at 12:59
FYI: In addition to the paper suggested above, they also wrote an accompanying [book](https://www.amazon.ca/Practical-Guide-Ecological-Modelling-Simulation/dp/1402086237) that is helpful for learning modelling with R. — Dan, Sep 06 '17 at 13:01
@BenBolker, I am investigating the reason of the error in python. As this equation came in study of physics, I am trying to check my physics. If all fails, I will ask for help here. Now, searching in stack I came across the linked problem and the sentence there "This also uses lsoda, but seems to be going off without a hitch" gave me motivation to check with `R` — Archimedes, Sep 06 '17 at 13:12
@Lyngbakr thanks for the reference, I will try to get the book. — Archimedes, Sep 06 '17 at 13:21

J_F · Accepted Answer · 2017-09-06T14:36:15.640

This should be a translation of the Python code to R

library(deSolve)

deriv <- function(t, state, parameters){
  with(as.list(c(state, parameters)),{

    M <- sqrt(1/3 * (1/2 * dy^2 + 1/2 * y^2))
    dz <- M*z # dz/dt
    ddy <-  -3* M * dy -  y # ddy/dt

    list(c(dy, ddy, dz))

  })
}

state <- c(y = 1,
           dy = -0.1,
           z = 1)

times <- seq(0, 100, length.out = 10001)

sol <- ode(func = deriv, y = state, times = times, parms = NULL)

y <- sol[, "y"]

dy <- sol[, "dy"]

z <- sol[, "z"]

M <- sqrt(1/3 * (1/2 *  dy^2 + 1/2*  y^2)) 

plot(times, z, col = "red", ylim = c(-1, 18), type = "l")
lines(times, y, col = "blue")
lines(times, M, col = "green")
grid()

There is a faster way to directly calculate M in R with this code:

library(deSolve)

deriv <- function(t, state, parameters){
  with(as.list(c(state, parameters)),{

    M <- sqrt(1/3 * (1/2 * dy^2 + 1/2 * y^2))
    dz <- M*z # dz/dt
    ddy <-  -3* M * dy -  y # ddy/dt

    list(c(dy, ddy, dz), M = M)

  })
}

state <- c(y = 1,
           dy = -0.1,
       z = 1)

times <- seq(0, 100, length.out = 10001)

sol <- ode(func = deriv, y = state, times = times, parms = NULL)

## save to file

write.csv2(sol,file = "path_to_folder/R_ODE.csv")

## plot

matplot(sol[,"time"], sol[,c("y", "z", "M")], type = "l")
grid()

How do I solve a second order differential equation in R?

1 Answers1

Linked