deSolve error for 0 values

Question

I tried to simulate a pendulum in R, using the package "deSolve" to solve the differential equations. The pendulum moves in two dimensions and includes the most important forces plus coriolisforce and shifting wind from the side. This is the script:

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)

So far it works, as I want it to. BUT if I set the starting values like this:

## 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

I get only NaN-values.

Is this a programming-issue, or did I do something wrong in the math/physiks part?

Answer 1

Try plugging in the new initial conditions to your derivative function like this:

y_strich(0, y_start, parms)

You'll find that you get this:

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

If you dig a little deeper using debug you'll find that the x component of e is zero, so the x component of r is zero. Now, consider this line:

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

You're dividing by zero! Unless you're Chuck Norris, that's not possible and lsoda understandably gets upset.