如何求解R中的二阶微分方程？

Question

我正在学习求解二阶微分方程（可能使用deSolve软件包）。我用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很容易解决这个问题，但是对于另一个类似但复杂的问题，Python显示出错误。我也尝试过python中的ode（vode/bdf），但问题仍然存在。现在，我想检查一下R如何处理这个问题。所以，如果有人能给我举个例子（基本上是一个代码翻译！）如何在R中解决这个问题，这样我就可以在R中尝试另一个问题，同时学习一些R（我知道这可能不是学习语言的理想方法）。在

我明白这个问题可能没有什么建设性的价值，但我只是一个新手，所以请您耐心等待！在

Answer 1

这应该是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()

有一种更快的方法可以使用以下代码直接计算R中的M：

^{pr2}$