如何解决R中的二阶微分方程? [英] How do I solve a second order differential equation in R?

问题描述

我正在学习 R 来求解二阶微分方程(可能使用deSolve软件包).我用python写成两个一阶微分方程写的，下面给出了

I am learning R to solve a second order differential equation(probably using deSolve package). Which I have written in python by writing it as two first order differential equations and is given below

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 (我知道这可能不是学习A的理想方法语言).

Python solve this problem easily, however for another similar but complicated problem python shows error. I have also tried ode(vode/bdf) in python but the problem is still there. Now, I would like to check how R work with the problem. So, I will be much obliged if someone please provide me an example(basically a code translation!) of how to solve this problem in R so that I can try the other one in R and also learn some R(I know that this may not be the ideal way to learn a language).

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

I understand that this question may have little constructive value, but I am just a novice in R, so please bear with me!

推荐答案

这应该是Python代码到R的翻译

This should be a translation of the Python code to 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:

There is a faster way to directly calculate M in R with this code:

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中的二阶微分方程?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！