SciPy "odeint"-Befehl, wie operiert Python dabei

Wenn du dir nicht sicher bist, in welchem der anderen Foren du die Frage stellen sollst, dann bist du hier im Forum für allgemeine Fragen sicher richtig.
Antworten
alex2007
User
Beiträge: 40
Registriert: Montag 14. April 2014, 10:08

Ich soll eine Aufgabe erledigen, bei der eine DGL 2. Ortnung gelöst werden soll. Diese überführe ich in ein System aus zwei DGLn 1. Ordnung und möchte diese mit Hilfe der 'odeint'-Operation lösen. Soweit so gut.

Nun habe ich aber folgendes Problem, dass die Lösung mir teilweise willkürlich erscheint. Betrachtet man folgende fast identische Codes:

Code: Alles auswählen

from numpy import *
from matplotlib import pyplot as plt
from scipy.integrate import odeint # Integrationsroutine importieren

def abl(y, t, A, B, omega):
    x_dot=y[1]                                     #x_dot=p
    p_dot=(y[0]**3)*4 -2*y[0]+A+B+sin(omega * t)   #p_dot=4x^3-2x+A+B+sin(wt)
    return array([x_dot, p_dot])
    
B = 0.1 
A = 0.1
omega = 1
y0 = array([0.0, 0.0]) # Anfangsbedingung: (x0, p0)
zeiten = linspace(0.0, 400*pi, 20) # Zeiten

y_t = odeint(abl, y0, zeiten, args=(A, B, omega)) # Integration der DGL
x_t = y_t[:, 0] # Auslesen der Spalten mit x
p_t = y_t[:, 1] # und p

print x_t
print p_t
plt.figure(1) 
plt.subplot(111, xscale="log", yscale="log", autoscale_on=False)
plt.axis([1e-300, 1e300, 1e-300, 2e300])
plt.plot(x_t, p_t)
plt.show()

Code: Alles auswählen

from numpy import *
from matplotlib import pyplot as plt
from scipy.integrate import odeint # Integrationsroutine importieren

def abl(y, t, A, B, omega):
    x_dot=y[1]                                     #x_dot=p
    p_dot=(y[0]**3)*4 -2*y[0]+A+B+sin(omega * t)   #p_dot=4x^3-2x+A+B+sin(wt)
    return array([x_dot, p_dot])
    
B = 0.1 
A = 0.1
omega = 1
y0 = array([0.0, 0.0]) # Anfangsbedingung: (x0, p0)
zeiten = linspace(0.0, 400*pi, 20) # Zeiten

y_t = odeint(abl, y0, zeiten, args=(A, B, omega)) # Integration der DGL
x_t = y_t[:, 0] # Auslesen der Spalten mit x
p_t = y_t[:, 1] # und p

plt.figure(1) 
plt.subplot(111, xscale="log", yscale="log", autoscale_on=False)
plt.axis([1e-300, 1e300, 1e-300, 2e300])
plt.plot(x_t, p_t)
plt.show()
Diese unterscheiden sich lediglich darin, dass einmal Werte berechnet werden, welche dann geplottet werden und einmal die Werte berechnet werden und dann erst ausgegeben werden, bevor sie geplottet werden. Theoretisch dürfte das am erhaltenden Ergebnis nichts ändern. Allerdings ergeben sich zwei unterschiedliche Plots. Warum?

PS: Ich weiß Sternchenimporte sind nicht so gut und die Kommentierung diente auch nur dem schnellen zweck und nicht der Form, die ich sonst hege. Darüber bitte ich hinwegzusehen :)

Danke für eure Hilfe!
BlackJack

@alex2007: Starte doch mal ein und das selber Programm mehrfach hintereinander. Dann wirst Du feststellen das selbst das nicht immer den selben Plot ergibt. Da wird also irgendwo etwas durch Zufall ausgewählt und kann dadurch zu unterschiedlichen Ergebnissen führen.
alex2007
User
Beiträge: 40
Registriert: Montag 14. April 2014, 10:08

BlackJack hat geschrieben:@alex2007: Starte doch mal ein und das selber Programm mehrfach hintereinander. Dann wirst Du feststellen das selbst das nicht immer den selben Plot ergibt. Da wird also irgendwo etwas durch Zufall ausgewählt und kann dadurch zu unterschiedlichen Ergebnissen führen.
Die Frage ist aber, warum?

Ich meine die Reihenfolge der Berechnung ist klar festgelegt. Demzufolge werden die Werte in einem Array abgespeichert, was aus Wertepaaren besteht. Siehst du irgendwo eine Stelle, wo die Problematik herühren könnte?

Ich meine eigentlich ist die Zuordnung klar.
BlackJack

@alex2007: Ich nehme mal an das die Integration nicht deterministisch ist und das da irgendwo Zufallswerte verwendet werden. Da müsste man sich halt durch die Dokumentation zur Implementierung durcharbeiten, was wenn ich das richtig sehe auf eine FORTRAN-Bibliothek `odepack` hinaus läuft.
alex2007
User
Beiträge: 40
Registriert: Montag 14. April 2014, 10:08

@ Blackjack: Könnte das http://mail.scipy.org/pipermail/scipy-u ... 21640.html der entscheidende Hinweis sein?

Es erscheint ja die im Link genannte Fehlermeldung, die auf die Genauigkeit der Integration abzielt.
ändert man die Anzahl der Zeiten in meinem Code von 20 auf bspw. 200, dann gibt es keine Fehlermeldung mehr. Kannst du das bestätigen, dass bei dir mehrere Durchläufe des Programms, dann stabile Lösungen anzeigen?
BlackJack

@alex2007: Das scheint ein stabiles Ergebnis zu geben. Allerdings immer noch mit Fehler/Warnmeldungen:

Code: Alles auswählen

$ python forum9.py 
 lsoda--  at current t (=r1), mxstep (=i1) steps   
       taken on this call before reaching tout     
      in above message,  i1 =       500
      in above message,  r1 =  0.2623602097886D+01
Excess work done on this call (perhaps wrong Dfun type).
Run with full_output = 1 to get quantitative information.
/usr/lib/pymodules/python2.7/matplotlib/path.py:136: RuntimeWarning: invalid value encountered in isfinite
  self.has_nonfinite = not np.isfinite(vertices).all()
/usr/lib/pymodules/python2.7/matplotlib/scale.py:130: RuntimeWarning: invalid value encountered in multiply
  a = self._handle_nonpos(a * 10.0)
Und der Plot sieht für mich auch eher nach moderner Kunst als nach etwas aus was mir was sagen möchte. :-)
alex2007
User
Beiträge: 40
Registriert: Montag 14. April 2014, 10:08

BlackJack hat geschrieben:@alex2007: Das scheint ein stabiles Ergebnis zu geben. Allerdings immer noch mit Fehler/Warnmeldungen:

Code: Alles auswählen

$ python forum9.py 
 lsoda--  at current t (=r1), mxstep (=i1) steps   
       taken on this call before reaching tout     
      in above message,  i1 =       500
      in above message,  r1 =  0.2623602097886D+01
Excess work done on this call (perhaps wrong Dfun type).
Run with full_output = 1 to get quantitative information.
/usr/lib/pymodules/python2.7/matplotlib/path.py:136: RuntimeWarning: invalid value encountered in isfinite
  self.has_nonfinite = not np.isfinite(vertices).all()
/usr/lib/pymodules/python2.7/matplotlib/scale.py:130: RuntimeWarning: invalid value encountered in multiply
  a = self._handle_nonpos(a * 10.0)
Und wenn du die Anzahl der Zeiten noch größer machst? Bei mir kommt auch bei 200 keine Fehlermeldung mehr!
Und der Plot sieht für mich auch eher nach moderner Kunst als nach etwas aus was mir was sagen möchte. :-)
Das ist bei theoretischer Physik nichts Unnormales :wink: (PS: wenn du Punkte statt eine Linie zeichnen lässt ist es nicht ganz so schlimm)
BlackJack

Die Warnungen sind dann weg.
Antworten