如何求解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我只是一個新手，所以請多多包涵！

來源

2017-09-06 Archimedes

[本文]（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()

來源

2017-09-06 13:33:43

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

回答

相關問題