Brauche Hilfe bei der Lösung eines zweiter Ordnung nicht-lineare ODE in python

Ich weiß wirklich nicht, wo beginnen mit diesem problem, da hatte ich noch nicht viel Erfahrung mit diesem, aber es ist notwendig, dies zu lösen, Teil des Projektes mit Hilfe eines Computers.

Ich habe einen 2nd order ODE die:

m = 1220

k = 35600

g = 17.5

a = 450000

und b ist zwischen 1000 und 10000 mit Schritten von 500.

x(0)= 0 

x'(0)= 5


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

Brauche ich das kleinste b derart, dass die Lösung ist nie positiv. Ich weiß, was das Diagramm Aussehen soll, aber ich weiß einfach nicht, wie zu verwenden odeint um eine Lösung der Differentialgleichung.
Dies ist der code den ich bisher:

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()

Die nicht wirklich produzieren viel von allem.

Irgendwelche Gedanken? Dank

InformationsquelleAutor Max Rackoff | 2013-11-04
  1. 10

    Zum lösen von second-order ODE mit scipy.integrate.odeint Sie schreiben sollte, es als ein system von first-order ODEs:

    Definiere ich z = [x', x], dann z' = [x'', x'], und das ist Ihr system! Natürlich, stecken Sie in Ihren realen Beziehungen:

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

    wird:

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

    z[1]' = z[0]

    Oder, nennen wir es doch einfach 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]'
                   ))

    Nun können Sie füttern die odeint als solche:

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

    (Die _ ist ein leerer Platzhalter für die x' variable, die wir gemacht)

    Minimieren b unterliegen der Einschränkung, dass das maximum von x ist immer negativ, Sie können scipy.optimize.minimize. Ich werde es implementieren, indem Sie wirklich die Maximierung der maximal x vorbehaltlich der Einschränkung, dass es weiterhin negativ, denn ich kann nicht daran denken, wie zur Minimierung eines Parameters, ohne in der Lage zu invertieren die Funktion.

    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)
    • Cool, wusste gar nicht, dass odeint unterstützt broadcasting der Argumente.
    • Eigentlich ohne Befestigung der return-format d(z) es nicht, ich soll zurück kommen und es beheben, aber habe es vergessen.
    • was bedeutet der "_, x" - Ding? vor allem die _-Teil
    • ist nur ein Variablenname ist, ist es äquivalent zu y, x = odeint(...) oder sogar z = odeint(...), da odeint gibt ein array mit der ersten dimension von 2 kann ich "entpacken" es. Ich entschied mich für die Verwendung _ weil ich habe keine Verwendung für die y Teil der Reihe, denn es stellt x', die ich nicht brauchen. Ich könnte alternativ dazu getan haben: x = odeint(...)[1,:]
    • so ist es im Grunde eine Puppe, die Sie nicht verwenden?
    • genau 🙂 Es ist nur eine Konvention, die nicht offizielle. " nicht benutzte Benennung von Variablen in python, Unterstrich in Python, _als Variablennamen in Python.

    InformationsquelleAutor askewchan
  2. 6

    Ich glaube nicht, können Sie Ihr problem lösen wie gesagt: Ihre ersten Bedingungen, mit x = 0 und x' > 0 implizieren, dass die Lösung positiv sein wird für einige Werte, die sehr nahe an den Ausgangspunkt. Es gibt also keine b für das die Lösung ist nie positiv...

    Abgesehen davon zu lösen zweiter Ordnung differentialgleichungen, müssen Sie zuerst, um zu umschreiben, dass es als ein system von zwei Erster Ordnung differentialgleichungen. Definition y = x' wir können schreiben Sie Ihre einzige Gleichung:

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

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

    Damit Ihre Funktion sollte in etwa so Aussehen:

    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])

    Und Sie können sich lösen und Plotten von Gleichungen zu tun:

    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()

    Brauche Hilfe bei der Lösung eines zweiter Ordnung nicht-lineare ODE in python

    • Ah, ich vergaß, dass die 5 sollte negativ sein. Auch ich verstehe nicht ganz, wie nennen x' = y oder y' = -b/my - k/mx - a/m*x**3 - g. Würde das nicht einfach werden, das beginnt mit strings?
    • Mit dem prime ' ist nur Schreibweise zu erklären, die Gleichung nicht, setzte Sie in Ihrem code. Der eigentliche code ist gegeben durch fun oder (d in meiner Antwort), wo statt der x' = y und y' = blah schreiben Sie fun([x, y]) = [y, blah]
    InformationsquelleAutor Jaime
Schreibe einen Kommentar