deSolve 0 值的错误

Question

我尝试在 R 中模拟钟摆，使用包 "deSolve" 来求解微分方程。摆在两个维度上移动，包括最重要的力加上科里奥利力和从侧面移动的风。这是脚本：

install.packages("deSolve")
library("deSolve")

#parameters
 parms=c( 
  xs=0.0,  #x-coordinate at rest
  ys=0.0,  #y-coordinate at rest
  kz=0.005,  #backwards-coefficient [N/m]
  m =0.01,  #mass pendulum [kg]
  kr=0.001,   #friction-coefficient [N/(m/s²)]
  wE=7.292115*10^-5, # angular speed earth (source: IERS)
  kw=0.002 # wind-coefficient
  )

  tmax=80  #end time [s]
  delta_t=0.05   #time steps [s]

# Initialisation
  t=seq(0,tmax,by=delta_t)    # time

## variable
 y=cbind(                                
  x=array(0,length(t)),     #x-coordinate         [m]
  y=array(0,length(t)),     #y-coordinate
  vx=array(0,length(t)),    #x-velocity   [m/s]
  vy=array(0,length(t))    #y-velocity
  )

## starting values
  y_start=c( 
  x=0.1,     #x-coordinate
  y=0.2,     #y-coordinate
  vx=0.1,     #x-velocity
  vy=-0.2      #y-velocity
  )

  y[1,]=y_start               #set start parameter



## function
y_strich=function(t, y_i,parms) 
{
  s = y_i[c(1,2)] # position at t
  v = y_i[c(3,4)] # velocity at t

  s_strich = v  # derivation of position

  e  = s - parms[c(1,2)]  # difference of position and rest = radius
  r = e  

  # WIND
  vw = parms["kw"]*(sin(t*0.3)) # windspeed
  Fw = y_i[3] * vw # windforce

  # CORIOLISFORCE
  rw = ((s/(2*pi*r))*360)*(pi/180) # rotation angle
  wg = rw / delta_t # angular velocity [in rad/s]
  Fc = (2*parms["m"]*(parms["wE"]*wg)) # Coriolisforce  

  # FRICTION AND BACKWARDS FORCE
  Fr = -v * parms["kr"] # friction
  Fz = -e * parms["kz"] # backwards force

  # sum of forces and velocity
  Fges = Fr + Fz + Fw + Fc    # sum of forces
  a = Fges / parms["m"] # accelariation


  v_strich = a

  return (list(c(s_strich, v_strich)))
}


# lsoda
y = lsoda(y=y_start, times=t, func=y_strich, parms=parms)

到目前为止，一切如我所愿。但是如果我像这样设置起始值：

## starting values
  y_start=c( 
  x=0.0,     #x-coordinate
  y=0.2,     #y-coordinate
  vx=0.0,     #x-velocity
  vy=-0.2      #y-velocity

我只得到 NaN 值。

这是编程问题，还是我在 math/physiks 部分做错了什么？

Answer 1

尝试像这样将新的初始条件插入到您的导数函数中：

y_strich(0, y_start, parms)

你会发现你得到了这个：

[[1]]
        vx          vy          vx          vy 
0.00000000 -0.20000000         NaN -0.07708315

如果使用 debug 进行更深入的挖掘，您会发现 e 的 x 分量为零，因此 [=17] 的 x 分量=] 为零。现在，考虑这一行：

rw = ((s/(2*pi*r))*360)*(pi/180) # rotation angle

你正在除以零！除非你是 Chuck Norris，否则这是不可能的，lsoda 可以理解地感到不安。

deSolve 0 值的错误

deSolve error for 0 values

physics

r

differential-equations