Ошибка решения для 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-значения.

Это проблема программирования, или я сделал что-то не так в части математики / физики?

Answer 1

Попробуйте подключить новые начальные условия к производной функции следующим образом:

y_strich(0, y_start, parms)

Вы обнаружите, что вы получаете это:

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

Если вы копаете немного глубже, используя debug вы обнаружите, что x составляющая e ноль, поэтому x составляющая r это ноль. Теперь рассмотрим эту строку:

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

Вы делите на ноль! Если вы не Чак Норрис, это невозможно и lsoda понятно расстраивается.