Potrzebujesz pomocy w rozwiązaniu problemu nieliniowego ODE drugiego rzędu w pytonie

Question

Nie wiem od czego zacząć, ponieważ nie miałem zbyt dużego doświadczenia w tym zakresie, ale konieczne jest rozwiązanie tej części projektu za pomocą komputer.

mam ODE 2-go rzędu, która ma:

m = 1220 

k = 35600 

g = 17.5 

a = 450000 

i b wynosi między 1000 a 10000, w odstępach co 500

x(0)= 0 

x'(0)= 5 


m*x''(t) + b*x'(t) + k*x(t)+a*(x(t))^3 = -m*g 

trzeba znaleźć najmniejszą B tak, że rozwiązanie nigdy nie jest pozytywne. Wiem, jak powinien wyglądać wykres, ale po prostu nie wiem, jak używać odinta, aby uzyskać rozwiązanie równania różniczkowego. Jest to kod mam tak daleko:

from numpy import *  
from matplotlib.pylab import *  
from scipy.integrate import odeint 

m = 1220.0 
k = 35600.0 
g = 17.5 
a = 450000.0 
x0= [0.0,5.0] 

b = 1000 

tmax = 10 
dt = 0.01 

def fun(x, t): 
    return (b*x[1]-k*x[0]-a*(x[0]**3)-m*g)*(1.0/m) 
t_rk = arange(0,tmax,dt) 
sol = odeint(fun, x0, t_rk) 
plot(t_rk,sol) 
show() 

Który tak naprawdę nie produkują dużo niczego.

Jakieś myśli? Dzięki

Answer 1

Nie wydaje mi się, że możesz rozwiązać swój problem zgodnie z opisem: twoje warunki początkowe z x = 0 i x' > 0 oznaczają, że rozwiązanie będzie dodatnie dla niektórych wartości bardzo zbliżonych do punktu początkowego. Tak więc nie ma żadnego b, dla którego rozwiązanie nigdy nie jest dodatnie ...

Pozostawiając to na bok, aby rozwiązać równanie różniczkowe drugiego rzędu, trzeba najpierw przepisać je jako układ dwóch równań różniczkowych pierwszego rzędu. Definiowanie y = x' możemy przepisać jednego równania jako:

x' = y 
y' = -b/m*y - k/m*x - a/m*x**3 - g 

x[0] = 0, y[0] = 5 

Więc funkcja powinna wyglądać mniej więcej tak:

def fun(z, t, m, k, g, a, b): 
    x, y = z 
    return np.array([y, -(b*y + (k + a*x*x)*x)/m - g]) 

I można rozwiązać i wykreślić swoje równania robi:

m, k, g, a = 1220, 35600, 17.5, 450000 
tmax, dt = 10, 0.01 
t = np.linspace(0, tmax, num=np.round(tmax/dt)+1) 
for b in xrange(1000, 10500, 500): 
    print 'Solving for b = {}'.format(b) 
    sol = odeint(fun, [0, 5], t, args=(m, k, g, a, b))[..., 0] 
    plt.plot(t, sol, label='b = {}'.format(b)) 
plt.legend() 

Answer 2

Aby rozwiązać problem ODE drugiej ręki za pomocą scipy.integrate.odeint, powinieneś napisać to jako system ODE pierwszej kolejności:

Zdefiniuję z = [x', x], następnie z' = [x'', x'], a to Twój system!Oczywiście, trzeba podłączyć rzeczywistych stosunków:

x'' = -(b*x'(t) + k*x(t) + a*(x(t))^3 + m*g)/m 

staje:

z[0]' = -1/m * (b*z[0] + k*z[1] + a*z[1]**3 + m*g)
z[1]' = z[0]

Albo po prostu nazwać d(z):

def d(z, t): 
    return np.array((
        -1/m * (b*z[0] + k*z[1] + a*z[1]**3 + m*g), # this is z[0]' 
        z[0]           # this is z[1]' 
        )) 

Teraz można karmić go do odeint jako takie:

_, x = odeint(d, x0, t).T 

(_ jest pusty symbol zastępczy x' zmiennej zrobiliśmy)

w celu zminimalizowania b przy ograniczeniu że maksimum x jest zawsze ujemne, możesz użyć scipy.optimize.minimize. Zaimplementuję go, maksymalizując maksymalnie wartość x, z zastrzeżeniem, że pozostaje ujemny, ponieważ nie mogę wymyślić jak zminimalizować parametr bez możliwości odwrócenia tej funkcji.

from scipy.optimize import minimize 
from scipy.integrate import odeint 

m = 1220 
k = 35600 
g = 17.5 
a = 450000 
z0 = np.array([-.5, 0]) 

def d(z, t, m, k, g, a, b): 
    return np.array([-1/m * (b*z[0] + k*z[1] + a*z[1]**3 + m*g), z[0]]) 

def func(b, z0, *args): 
    _, x = odeint(d, z0, t, args=args+(b,)).T 
    return -x.max() # minimize negative max 

cons = [{'type': 'ineq', 'fun': lambda b: b - 1000, 'jac': lambda b: 1}, # b > 1000 
     {'type': 'ineq', 'fun': lambda b: 10000 - b, 'jac': lambda b: -1}, # b < 10000 
     {'type': 'ineq', 'fun': lambda b: func(b, z0, m, k, g, a)}] # func(b) > 0 means x < 0 

b0 = 10000 
b_min = minimize(func, b0, args=(z0, m, k, g, a), constraints=cons)