erreur deSolve pour 0 valeurs

Question

J'ai essayé de simuler un pendule dans R, en utilisant le paquet "deSolve" pour résoudre les équations différentielles. Le pendule se déplace en deux dimensions et comprend les forces les plus importantes, plus coriolisforce et le vent qui se déplace du côté. Voici le 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) 

Jusqu'à présent, cela fonctionne, comme je veux. MAIS si je règle les valeurs de départ comme ceci:

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

Je reçois seulement des valeurs NaN.

Est-ce un problème de programmation, ou ai-je fait quelque chose de mal dans la partie maths/kinésithérapeutes?

Answer 1

Essayez de brancher dans les nouvelles conditions initiales à votre fonction dérivée comme ceci:

y_strich(0, y_start, parms) 

Vous trouverez que vous obtenez ceci:

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

Si vous creusez un peu plus profond à l'aide que vous debug trouveront que le composant x de e est zéro, donc le composant x de r est zéro. Maintenant, considérez cette ligne:

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

Vous divisez par zéro! Sauf si vous êtes Chuck Norris, ce n'est pas possible et lsoda se fâche naturellement.