0值的deSolve错误

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

我只得到南的价值。

这是程序问题，还是我在数学/物理部分做错了什么？

Answer 1

尝试将新的初始条件插入到派生函数中，如下所示：

y_strich(0, y_start, parms)

你会发现你得到了这个：

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

如果使用debug深入挖掘，您会发现e的x组件为零，因此r的x组件为零。现在，考虑一下这一行：

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

你被零除了！除非你是查克·诺里斯，否则这是不可能的，lsoda也会心烦意乱。