Разработать и реализовать метод для расчета минимальной скорости ветра W, для которой небольшое возмущение θ (0) = 10–3 имеет коэффициент увеличения 100 или более в пределах 0,5 × 10–3 км / ч.
Здесь много дополнительной информации. У меня все это есть, но я думаю, что мой мозг немного жарен, потому что я не могу найти правильный способ сделать эту простую часть.
Я решил для диапазона скорости ветра, который имеет коэффициент увеличения (ЭДС) 100 (57,4 - 57,5). Я думал о попытке использовать метод деления пополам, но я не совсем уверен, как реализовать его в этом случае.
Моя проблема с методом деления пополам состоит в том, что я вызываю так много функций, что я не уверен, что разделить пополам.
Я думал, смогу ли я каким-то образом пройти через ЭДС как допуск и изменять скорость ветра (W) от 57,4 до тех пор, пока он не достигнет значения tol. Я просто застрял, кодблок вещь?
Даже малейший толчок в правильном направлении будет полезен на этом этапе.
test = 10**-3
wind = [55,56,57,57.2,57.4]
def rk4step0(t,w,h):
#one step of the Runge-Kutta order 4 method
s1 = ydot0(t,w)
s2 = ydot0(t+h/2,w+h*s1/2.)
s3 = ydot0(t+h/2,w+h*s2/2.)
s4 = ydot0(t+h,w+h*s3)
out = w + h*(s1+2*s2+2*s3+s4)/6.
#print(out)
return out
def ydot0(t,y):
leng = 6; a = 0.2; W = wind[0]; omega = 2*pi*38/60.
K=1000
M=2500
Koverm=K/M
#Koverm=omega**2
linear = False #toggle to switch model
if linear:
F1 = y[0]-leng*sin(y[2])
F2 = y[0]+leng*sin(y[2])
else:
F1 = exp(a*(y[0]-leng*sin(y[2])))-1
F1 /= a
F2 = exp(a*(y[0]+leng*sin(y[2])))-1
F2 /= a
ydot = empty(4)
ydot[0] = y[1]
ydot[1] = -0.01*y[1]-Koverm*(F1+F2)+0.2*W*sin(omega*t)
ydot[2] = y[3]
ydot[3] = -0.01*y[3]+3*Koverm*cos(y[2])*(F1-F2)/leng
return ydot
def rk4step1(t,w,h):
#one step of the Runge-Kutta order 4 method
s1 = ydot1(t,w)
s2 = ydot1(t+h/2,w+h*s1/2.)
s3 = ydot1(t+h/2,w+h*s2/2.)
s4 = ydot1(t+h,w+h*s3)
out = w + h*(s1+2*s2+2*s3+s4)/6.
#print(out)
return out
def ydot1(t,y):
leng = 6; a = 0.2; W = wind[1]; omega = 2*pi*38/60.
K=1000
M=2500
Koverm=K/M
#Koverm=omega**2
linear = False #toggle to switch model
if linear:
F1 = y[0]-leng*sin(y[2])
F2 = y[0]+leng*sin(y[2])
else:
F1 = exp(a*(y[0]-leng*sin(y[2])))-1
F1 /= a
F2 = exp(a*(y[0]+leng*sin(y[2])))-1
F2 /= a
ydot = empty(4)
ydot[0] = y[1]
ydot[1] = -0.01*y[1]-Koverm*(F1+F2)+0.2*W*sin(omega*t)
ydot[2] = y[3]
ydot[3] = -0.01*y[3]+3*Koverm*cos(y[2])*(F1-F2)/leng
return ydot
def rk4step2(t,w,h):
#one step of the Runge-Kutta order 4 method
s1 = ydot2(t,w)
s2 = ydot2(t+h/2,w+h*s1/2.)
s3 = ydot2(t+h/2,w+h*s2/2.)
s4 = ydot2(t+h,w+h*s3)
out = w + h*(s1+2*s2+2*s3+s4)/6.
#print(out)
return out
def ydot2(t,y):
leng = 6; a = 0.2; W = wind[2]; omega = 2*pi*38/60.
K=1000
M=2500
Koverm=K/M
#Koverm=omega**2
linear = False #toggle to switch model
if linear:
F1 = y[0]-leng*sin(y[2])
F2 = y[0]+leng*sin(y[2])
else:
F1 = exp(a*(y[0]-leng*sin(y[2])))-1
F1 /= a
F2 = exp(a*(y[0]+leng*sin(y[2])))-1
F2 /= a
ydot = empty(4)
ydot[0] = y[1]
ydot[1] = -0.01*y[1]-Koverm*(F1+F2)+0.2*W*sin(omega*t)
ydot[2] = y[3]
ydot[3] = -0.01*y[3]+3*Koverm*cos(y[2])*(F1-F2)/leng
return ydot
def rk4step3(t,w,h):
#one step of the Runge-Kutta order 4 method
s1 = ydot3(t,w)
s2 = ydot3(t+h/2,w+h*s1/2.)
s3 = ydot3(t+h/2,w+h*s2/2.)
s4 = ydot3(t+h,w+h*s3)
out = w + h*(s1+2*s2+2*s3+s4)/6.
#print(out)
return out
def ydot3(t,y):
leng = 6; a = 0.2; W = wind[3]; omega = 2*pi*38/60.
K=1000
M=2500
Koverm=K/M
#Koverm=omega**2
linear = False #toggle to switch model
if linear:
F1 = y[0]-leng*sin(y[2])
F2 = y[0]+leng*sin(y[2])
else:
F1 = exp(a*(y[0]-leng*sin(y[2])))-1
F1 /= a
F2 = exp(a*(y[0]+leng*sin(y[2])))-1
F2 /= a
ydot = empty(4)
ydot[0] = y[1]
ydot[1] = -0.01*y[1]-Koverm*(F1+F2)+0.2*W*sin(omega*t)
ydot[2] = y[3]
ydot[3] = -0.01*y[3]+3*Koverm*cos(y[2])*(F1-F2)/leng
return ydot
def rk4step4(t,w,h):
#one step of the Runge-Kutta order 4 method
s1 = ydot4(t,w)
s2 = ydot4(t+h/2,w+h*s1/2.)
s3 = ydot4(t+h/2,w+h*s2/2.)
s4 = ydot4(t+h,w+h*s3)
out = w + h*(s1+2*s2+2*s3+s4)/6.
#print(out)
return out
def ydot4(t,y):
leng = 6; a = 0.2; W = wind[4]; omega = 2*pi*38/60.
K=1000
M=2500
Koverm=K/M
#Koverm=omega**2
linear = False #toggle to switch model
if linear:
F1 = y[0]-leng*sin(y[2])
F2 = y[0]+leng*sin(y[2])
else:
F1 = exp(a*(y[0]-leng*sin(y[2])))-1
F1 /= a
F2 = exp(a*(y[0]+leng*sin(y[2])))-1
F2 /= a
ydot = empty(4)
ydot[0] = y[1]
ydot[1] = -0.01*y[1]-Koverm*(F1+F2)+0.2*W*sin(omega*t)
ydot[2] = y[3]
ydot[3] = -0.01*y[3]+3*Koverm*cos(y[2])*(F1-F2)/leng
return ydot
def mag_factor(y):
test = []
#for i in range(4):
one=y[1:None,1]
two=y[1:None,3]
temp = max(two)
# test.append(temp)
return(temp)
a,b = tacoma([0,1000],[0,0,test,0],25000,time_int=rk4step0)
maxB=mag_factor(b)
EMF = maxB/test
print("Wind = 55 EMF",EMF)
print("Initial θ=10^-3 Max θ= ", maxB)
a,b = tacoma([0,1000],[0,0,test,0],25000,time_int=rk4step1)
maxB=mag_factor(b)
EMF = maxB/test
print("Wind = 56 EMF for",EMF)
print("Initial θ=10^-3 Max θ= ", maxB)
a,b = tacoma([0,1000],[0,0,test,0],25000,time_int=rk4step2)
maxB=mag_factor(b)
EMF = maxB/test
print("Wind = 57 EMF",EMF)
print("Initial θ=10^-3 Max θ= ", maxB)
a,b = tacoma([0,1000],[0,0,test,0],25000,time_int=rk4step3)
maxB=mag_factor(b)
EMF = maxB/test
print("Wind = 57.2 EMF",EMF)
print("Initial θ=10^-3 Max θ= ", maxB)
a,b = tacoma([0,1000],[0,0,test,0],25000,time_int=rk4step4)
maxB=mag_factor(b)
EMF = maxB/test
print("Wind = 57.4 EMF",EMF)
print("Initial θ=10^-3 Max θ= ", maxB)
print("Minimum wind is somewhere between [57.4,57.5]")