如何使用 deSolve 在 R 中数值求解 Lane-Emden 方程答案

【问题标题】：How to to solve Lane-Emden equation numerically in R using deSolve如何使用 deSolve 在 R 中数值求解 Lane-Emden 方程
【发布时间】：2019-09-07 16:09:55
【问题描述】：

我想使用 R 数值求解一个称为 Lane-Emden 方程的特定微分方程。 https://en.wikipedia.org/wiki/Lane%E2%80%93Emden_equation#Numerical_solutions

我很困惑如何在 deSolve 中为 ODE 求解器制定问题。我有一些代码，它给出了一个我不理解且不知道如何修复的错误。

我已尝试改编 K Soetart 等人的《Solving Differential Equations in R》一书中的示例。

library(deSolve)
x <- seq(0,20,0.1)
n <- 1
y_ini <- c(phi=0, theta=1)

derivs <- function(x,y,parms){
  with(as.list(y),{
  dtheta <- -phi/(xi^2)
  dphi <- theta^n * xi^2
  list(c(dtheta,dphi))})
}

out <- ode(y=y_ini,times <- x,func = derivs,parms=NULL)

eval(substitute(expr),data,enclos= parent.frame()) 中的错误：找不到对象 'xi' 调用：... Func2 -> func -> with -> with.default -> eval -> eval 执行停止

【问题讨论】：

这个错误来自于你的函数涉及一个对象xi，但是你没有在任何地方定义它。
解决方案就像在derivs <- function(xi,y,parms){ 中添加i 一样简单。
另请参阅stackoverflow.com/q/42404333/3088138，了解在 python 中的成功集成以及使用的解决方法。
"这个错误来自于你的函数涉及一个对象 xi，但是你没有在任何地方定义它。"是的，它确实。这就是为什么我认为我没有正确地为 deSolve 提出问题。看看维基百科中的问题定义，你能告诉我如何重新表述它吗？
"解决方案可以像在 derivs 中添加 i 一样简单 ...withCallingHandlers -> withVisible -> eval -> eval -> ode -> -> ldoda

标签： r ode

【解决方案1】：

执行 cmets 中提到的所有更改：

在derivs的参数列表中使用xi，而不是x，
在状态向量和导数中同步theta和phi的顺序
通过将第一次更改为某个小的正值来避免除以零

给出代码

library(deSolve)
x <- seq(0.0001,20,0.1)
n <- 1
y_ini <- c(theta=1, phi=0 )

derivs <- function(xi,y,parms){
  with(as.list(y),{
  dtheta <- -phi/(xi^2)
  dphi <- theta^n * xi^2
  list(c(dtheta,dphi))})
}

out <- ode(y=y_ini,times <- x,func = derivs,parms=NULL)
print(out)
theta <- out[,"theta"]
plot(times, theta)

然后产生情节

【讨论】：