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

问题描述：

我正在学习R解决一个二阶微分方程（可能使用deSolve包）。这是我写这两个一阶微分方程用Python编写的，并给出以下如何求解R中的二阶微分方程？

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中的ode（vode/bdf），但问题仍然存在。现在，我想查看如何与问题一起工作R。所以，我将非常感激，如果有人请给我一个例子，如何在[R解决这个问题，这样我可以尝试另一种在[R，也学会（基本上是一个代码翻译！）一些R（我知道这可能不是学习语言的理想方式）。

据我所知，这个问题可能没有什么建设性的价值，但我在[R我只是一个新手，所以请多多包涵！

[本文]（https://journal.r-project.org/archive/2010-2/RJournal_2010-2_Soetaert~et~al.pdf）给出了一个很好的介绍关于如何在'R'中使用'deSolve'，并且通过一个简单的例子来介绍。不过，我建议你考虑为什么Python会显示错误，以及为什么你认为使用'R'会解决你的问题。 –

我目前正试图了解为什么python显示错误，如果我找不到答案，我会在stackoverflow中寻求帮助。然而，寻找帮助，我碰到这个职位：https://stackoverflow.com/questions/16973036/odd-scipy-ode-integration-error 并试图实施ode（vode/bdf），但问题仍然存在。所以现在尝试R – Archimedes

我认为如果没有任何有关错误的信息，任何人都很难帮助您。您链接的帖子建议使用'scipy.integrate.ode'来指定是使用僵硬还是非僵硬的方法。我猜你已经做到了？ –

答

这应该是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中的二阶微分方程？

相关推荐