0 값에 대한 deSolve 오류

"deSolve"패키지를 사용하여 R의 진자를 미분 방정식으로 풀려고했습니다. 진자는 2 차원으로 움직이며, 가장 중요한 힘과 코리올리 스 (coriolisforce) 및 측면으로부터의 바람의 이동을 포함합니다. 이 스크립트입니다 : 내가 원하는대로0 값에 대한 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

만 얻을 수 있습니다.

프로그래밍 문제입니까? math/physiks 부분에서 잘못된 점이 있습니까?

출처

2017-10-04 S. Gey

이처럼 파생 함수에 새로운 초기 조건에 연결해보십시오 :

y_strich(0, y_start, parms)

당신은 당신이 얻을 찾을 수 있습니다 : 당신이 발굴되면 조금 더 깊이 당신을 debug를 사용

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

e의 x 구성 요소는 0이므로 r의 x 구성 요소는 0입니다. 이제 다음 줄을 고려해보십시오.

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

0으로 나누고 있습니다. 니가 척 노리스가 아니라면 그건 불가능하니 lsoda이 당연히 당황해.

출처

2017-10-06 18:50:22 Lyngbakr

0 값에 대한 deSolve 오류

답변

관련 문제