Rの2階微分方程式を解くにはどうすればよいですか？

Question

私は二次微分方程式（おそらくdeSolveパッケージを使用）を解くためにRを学習しています。私は2つの1次微分方程式としてそれを書くことによって、Pythonで書かれている。しかし、他の類似しているが、複雑な問題のpythonのために、簡単にこの問題を解決するためのPython

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でode（vode/bdf）を試しましたが、問題はまだあります。今、私はどのように問題を解決するにはRを確認したいと思います。誰かが私の例（基本的にコード変換！）私はRに他のものを試してみても学ぶことができるようにRでこの問題を解決する方法のを提供してくださいあれば、私はずっと義務がありますいくつかR（これは言語を学ぶのに理想的な方法ではないことがわかります）。そう私と一緒に負担してください、

私はこの質問は少し建設的価値を有することができることを理解し、私はRにちょうど初心者です！

Answer 1

これは、直接このコードでMでRを計算するより高速な方法はあり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() 

にPythonコードの翻訳する必要があります：

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