Comment résoudre une équation différentielle du second ordre dans R?

Question

J'apprends R pour résoudre une équation différentielle de second ordre (probablement en utilisant le paquet deSolve). Ce que j'ai écrit en python en écrivant comme deux premier ordre des équations différentielles et est donné ci-dessous

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 résoudre ce problème facilement, mais pour un autre python problème similaire mais complexe indique l'erreur. J'ai aussi essayé ode (vode/bdf) en python mais le problème est toujours là. Maintenant, je voudrais vérifier comment R de travailler avec le problème. Donc, je serai bien obligé si quelqu'un s'il vous plaît me donner un exemple (essentiellement une traduction de code!) De la façon de résoudre ce problème dans R pour que je puisse essayer l'autre dans R et aussi apprendre certains R (Je sais que ce n'est peut-être pas le moyen idéal pour apprendre une langue).

Je comprends que cette question peut avoir peu de valeur constructive, mais je suis juste un novice dans R, donc s'il vous plaît garder avec moi!

Answer 1

Cela devrait être une traduction du code Python 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() 

Il y a un moyen plus rapide de calculer directement M dans R avec ce 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()