deSolve爲0值的錯誤

我試圖在R中模擬擺錘，使用軟件包「deSolve」來求解微分方程。擺錘在兩個方向上移動，包括最重要的力量和科里奧利力量，並從側面移動風。這是腳本：deSolve爲0值的錯誤

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值。

這是一個編程問題，還是我在數學/生理科中做錯了什麼？

來源

2017-10-04 S. Gey

嘗試將在新的初始條件的衍生功能是這樣的：

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

你被零除！除非你是Chuck Norris，否則這是不可能的，lsoda會變得不安。

來源

2017-10-06 18:50:22 Lyngbakr

deSolve爲0值的錯誤

回答

相關問題