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

Question

我正在学习 R 来求解二阶微分方程（可能使用 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（我知道这可能不是学习语言的理想方式）。

我明白这个问题可能没有什么建设性价值，但我只是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有一个更快的方法，代码如下：

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

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

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

python

r

ode

differential-equations