使用 R 求解 ODE

Question

我正在尝试使用 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)

弹出的问题是：

The number of derivatives returned by func() (4) must equal the length of the initial conditions vector (2)

我不是为什么它提示输入了 4 个导数，而我只输出了两个

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

使用 R 求解 ODE

Using R to solve an ODE

r

ode