使用 R 求解 ODE答案

【问题标题】：Using R to solve an ODE使用 R 求解 ODE
【发布时间】：2018-09-23 03:20:35
【问题描述】：

我正在尝试使用 R 解决 ODE，因为我无法访问 matlab

方程是

dh/dt = 0.1*v/(pi*(2*10*h-h^2))

v^2 = (-0.1*v/(pi*(2*10*h-h^2))^2) + 2*9.81*h))

vdpol <- function (h, v, t) ( 
      list(c (
        -0.1*v/(pi*(2*10*h-h^2)),
          (v^2 = (-0.1*v/(pi*(2*10*h-h^2))^2) + 2*9.81*h))

              ))
library(deSolve)
yini <- (c(h = 20, v=0))
nonstiff <- ode(y = yini, func = vdpol,
                times= seq(0, 30, by = 0.01),
                parms = 1)

弹出的问题是：

func() (4) 返回的导数个数必须等于初始条件向量的长度 (2)

我不是为什么当我只输出两个时它建议输入了 4 个导数

【问题讨论】：

请不要交叉发布到 rhelp 和 StackOverflow。

标签： r ode

【解决方案1】：

IIUC，ode 需要具有给定格式的函数。参考here

因此，如果您更改 vdpol() 函数以满足预期格式，它应该会运行。一般格式为func(t, state, parameters)，其中state 控制变量，parameters 控制其他参数。

vdpol <- function (t, state, parameters) ( 
    with(as.list(c(state)), {
    return(list(c (-0.1*v/(pi*(2*10*h-h^2)),
                  (-0.1*v/(pi*(2*10*h-h^2))^2) + 2*9.81*h)))
    })
)

state = c(h = 10, v = 0)
times= seq(0, 30, by = 0.01)
out <- ode(y = state, times = times, func = vdpol, parms = c())
plot(out)

【讨论】：