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值。
這是一個編程問題,還是我在數學/生理科中做錯了什麼?