Die Mathe-Redaktion - 13.11.2019 16:56 - Registrieren/Login
Auswahl
ListenpunktHome
ListenpunktAktuell und Interessant ai
ListenpunktArtikelübersicht/-suche
ListenpunktAlle Links / Mathe-Links
ListenpunktFach- & Sachbücher
ListenpunktMitglieder / Karte / Top 15
ListenpunktRegistrieren/Login
ListenpunktArbeitsgruppen
Listenpunkt? im neuen Schwätz
ListenpunktWerde Mathe-Millionär!
ListenpunktFormeleditor fedgeo
Schwarzes Brett
Aktion im Forum
Suche
Stichwortsuche in Artikeln und Links von Matheplanet
Suchen im Forum
Suchtipps

Bücher
Englische Bücher
Software
Suchbegriffe:
Mathematik bei amazon
Naturwissenschaft & Technik
In Partnerschaft mit Amazon.de
Kontakt
Mail an Matroid
[Keine Übungsaufgaben!]
Impressum

Bitte beachten Sie unsere Nutzungsbedingungen, die Distanzierung, unsere Datenschutzerklärung und
die Forumregeln.

Sie können Mitglied werden. Mitglieder können den Matheplanet-Newsletter bestellen, der etwa alle 2 Monate erscheint.

Der Newsletter Okt. 2017

Für Mitglieder
Mathematisch für Anfänger
Wer ist Online
Aktuell sind 1002 Gäste und 18 Mitglieder online.

Sie können Mitglied werden:
Klick hier.

Über Matheplanet
 
Zum letzten Themenfilter: Themenfilter:
Matroids Matheplanet Forum Index
Moderiert von matroid
Mathematik » Numerik & Optimierung » Python: symplektisches Euler-Verfahren
Druckversion
Druckversion
Autor
Universität/Hochschule J Python: symplektisches Euler-Verfahren
lissy1234567
Aktiv Letzter Besuch: in der letzten Woche
Dabei seit: 01.09.2017
Mitteilungen: 420
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Themenstart: 2019-10-19 13:04


Hallo,

ich programmiere gerade das symplektische Euler-Verfahren am Beispiel des Pendels und habe eine Frage zur Ausgabe.
Hier zunächst  mein Code:
python
import math
import matplotlib.pyplot as plt  
from matplotlib import animation      
 
def approx(t, q0):
    """ Evaluates the analytical approximation. """
    return q0*np.cos(t)
 
def energy (q,p):
    return((p**2)/2 -np.cos(q)+1)
 
def symplektisch(p0,q0,T,n):
    dt=T/float(n)
    t=0
    p,q = p0,q0
    energie = energy(q0,p0)
    motion=[[]for i in range(4)]
    motion[0].append(p)
    motion[1].append(q)
    motion[2].append(t)
    motion[3].append(energie)
    while t<= T:
        p = p + (-np.sin(q))*dt
        q = q + p*dt
        t = t+dt
        energie = energy(q,p)
        motion[0].append(p)
        motion[1].append(q)
        motion[2].append(t)
        motion[3].append(energie)
    return motion
 
p0 = 0
q0 = np.pi/12
T = 40
n = 10000
 
bew = symplektisch(p0,q0,T,n)
 
def kinet(p):
    return(np.array(p)**2/2)
def pot(q):
    return(1-np.cos(q))
kin1 = kinet(bew[0])
pot1 = pot(bew[1])

Bereits geplottet habe ich die analytische Lösung, das Phasenraum-Diagramm, die Gesamtenergie sowie potentielle und kinetische Energie.

Nun würde ich gerne den Energiefehler plotten, weiß aber nicht, wie ich das anstellen soll...
Weiß jemand, wie ich das machen kann?



  Profil  Quote  Link auf diesen Beitrag Link
lissy1234567
Aktiv Letzter Besuch: in der letzten Woche
Dabei seit: 01.09.2017
Mitteilungen: 420
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.1, vom Themenstarter, eingetragen 2019-10-19 13:52


Kann mir jemand helfen?



  Profil  Quote  Link auf diesen Beitrag Link
rlk
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 16.03.2007
Mitteilungen: 10562
Aus: Wien
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.2, eingetragen 2019-10-20 21:24


Hallo Lissy,
bei dem hier angenommenen reibungsfreien Pendel bleibt die Gesamtenergie erhalten, Du kannst daher bew[4] mit energy(q0, p0) vergleichen. Die Plots würden wohl einige hier interessieren.

Gestern schrieb Ken Thompson auf
TUHS darüber, dass die Energie der Raumschiffe in seinem Spiel
Space Travel nicht konstant blieb. Eine leider nicht überlieferte Formel von Richard Hamming behob diesen Fehler.

Servus,
Roland




  Profil  Quote  Link auf diesen Beitrag Link
lissy1234567
Aktiv Letzter Besuch: in der letzten Woche
Dabei seit: 01.09.2017
Mitteilungen: 420
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.3, vom Themenstarter, eingetragen 2019-10-21 10:44


Hallo Roland,

danke für deine Hilfe!
Inzwischen habe ich den Code ein bisschen geändert. Könnte man den Energiefehler auch folgendermaßen darstellen/berechnen?
python
import math
import matplotlib.pyplot as plt  
from matplotlib import animation      
 
# analytische Approximation
def approx(t, q0):
    return q0*np.cos(t)
 
# Symplektisches Euler-Verfahren anwenden
def symplektisch(p0,q0,T,n):
    # Werte setzen
    dt=T/float(n)
    t=0
    p,q = p0,q0
    energie = energy(q0,p0)
 
    # Daten in motion sammeln
    motion=[[]for i in range(3)]
    motion[0].append(p)
    motion[1].append(q)
    motion[2].append(t)
 
    # eigentliches Verfahren
    while t<= T:
        p = p + (-np.sin(q))*dt
        q = q + p*dt
        t = t+dt
 
        motion[0].append(p)
        motion[1].append(q)
        motion[2].append(t)
 
    return motion
 
# Werte wählen
p0 = 0
q0 = np.pi/12
T = 40
n = 10000
 
# Anwendung des Verfahrens und analytische Approximation
bew = symplektisch(p0,q0,T,n)
 
# Umdefinieren
pe, qu, te = bew
analy = approx(te,q0)
 
plt.plot(te, analy, 'black', label = 'Analytische Lösung' )
plt.plot(te, qu, 'green', label = 'Lösung mit Verfahren')
plt.xlabel('Zeit t')
plt.ylabel('Winkel q')
plt.legend(loc = 2)
plt.title('Analytische vs. numerische Lösung')
 
fig, ax = plt.subplots(1,2, figsize=(10,5))
plt.subplots_adjust( right = 1.5, wspace =0.5)
 
ax[0].plot(qu, pe, 'purple')
ax[0].set_xlabel('Winkel q')
ax[0].set_ylabel('Winkelgeschwindigkeit p')
 
ax[1].plot(te, qu, 'green')
ax[1].set_xlabel('Zeit t')
ax[1].set_ylabel('Winkel q')
 
plt.show()
#plt.legend()
plt.show()
 
# kinetische und potentielle Energie definieren (H = K + P)
def kinet(p):
    return(np.array(p)**2/2)
def pot(q):
    return(1-np.cos(q))
 
# Anwendung auf numerische Lösung
kin1 = kinet(pe)
pot1 = pot(qu)
 
# Plot
plt.plot(te, kin1, label = 'Kinetische Energie')
plt.plot(te, pot1, label = 'Potentielle Energie')
plt.plot(te, kin1+pot1, label = 'Gesamtenergie')
plt.title('Gesamtenergie bei $q_0=%.0f^\circ$'%(q02*180/np.pi))
plt.legend()
plt.show()
 
# Anwendung auf analytische Lösung
kin2 = kinet(analy)
pot2 = pot(analy)
 
# Fehlerberechnung
error = ((kin2+pot2)-(kin1+pot1))*100
 
# Plot
plt.plot(te, error, label = 'Fehler')
plt.title('Energiefehler in Prozent')
plt.legend()
plt.show()

Das heißt, zusammengefasst, indem man weiß, dass die Gesamtenergie aus kinetischer und potentieller Energie besteht und jeweils die zwei Formeln auf die numerische Lösung und auf die analytische Lösung anwendet?
Könnte da bitte jemand mal drüberschauen? Ich bin in python oder allgemein im Programmieren nicht wirklich bewandert.

Danke schonmal :))



  Profil  Quote  Link auf diesen Beitrag Link
rlk
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 16.03.2007
Mitteilungen: 10562
Aus: Wien
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.4, eingetragen 2019-10-22 00:08


Hallo Lissy,
ich hatte gemeint, dass Du die Energie der numerischen Lösung mit der Energie im Startpunkt vergleichst.

Die analytische Lösung ist nur eine Näherung, daher geht beim Vergleich der Energie mit der numerischen Lösung auch der Approximationsfehler ein.

Mit error berechnest Du die hundertfache Differenz der Energien, sie hat die Einheit der Energie, und nicht Prozent, wie Du im Titel des Plots (den ich gerne sehen würde) schreibst.

Ist das Dein vollständiges Programm? Ich frage mich nämlich, wo np definiert wird.

Servus,
Roland



  Profil  Quote  Link auf diesen Beitrag Link
lissy1234567
Aktiv Letzter Besuch: in der letzten Woche
Dabei seit: 01.09.2017
Mitteilungen: 420
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.5, vom Themenstarter, eingetragen 2019-10-22 12:01


Vielen Dank für den Hinweis mit den Prozent, das wusste ich nicht :)
Das np soll für numpy stehen. Das Programm würde also ohn 'import numpy as np' nicht funktionieren, was ich aber nicht gemerkt hatte, da ich davor schon andere Jupyter Notebooks offen hatte und das Programm auch so funktionierte.

Hier die Plots, den du gerne sehen würdest :)



Sieht das gut aus?

Dann noch kurz zu deinem Vorschlag, die Start-Energie mit der numerischen zu vergleichen. Wie kann ich das denn plotten? Denn die Startenergie ist ja nur ein Wert, die numerische Lösungen sind ja aber mehrere...ich habs nicht hinbekommen, daher hoffe ich auf deine Hilfe :)



  Profil  Quote  Link auf diesen Beitrag Link
rlk
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 16.03.2007
Mitteilungen: 10562
Aus: Wien
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.6, eingetragen 2019-10-23 00:06


Hallo Lissy,
dass np für numpy steht, hatte ich vermutet. Bitte versuche, die hier veröffentlichten Programme zu testen, um Lesern, die mit den Programmen arbeiten wollen, die Fehlersuche zu ersparen.

Der erste Plot sieht gut aus, wenn Du genau hinsiehst, erkennst Du, dass die Gesamtenergie leichte Schwankungen aufweist, was nicht passieren sollte. Bist Du sicher, dass Dein Integrationsverfahren symplektisch ist?

Zum zweiten Plot schreibe ich später noch etwas.

Servus,
Roland



  Profil  Quote  Link auf diesen Beitrag Link
lissy1234567
Aktiv Letzter Besuch: in der letzten Woche
Dabei seit: 01.09.2017
Mitteilungen: 420
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.7, vom Themenstarter, eingetragen 2019-10-23 11:39


Ich werde in Zukunft drauf achten, hab es diesmal wirklich nicht bemerkt.

Ja, ich sehe die Schwankungen auch...im Buch 3x9 Numerik von Sören Bartels habe ich eine Abbildung gefunden, die allerdings auch solche minimale Schwankungen zeigt und darüber eine mathematische Begründung, auf Seite 238 und 239 (weiß nicht, ob ich Screenshots davon hier reinstellen darf). Beim impliziten Mittelpunktverfahren bspw. ist die Gerade dabei gerade und exakt. Mein Verfahren ist das symplektische Euler-Verfahren, das ich verwendet habe und wie der Name schon sagt ist das symplektisch :D

Danke für deine gute Hilfe :)



  Profil  Quote  Link auf diesen Beitrag Link
rlk
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 16.03.2007
Mitteilungen: 10562
Aus: Wien
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.8, eingetragen 2019-10-27 21:40


Hallo Lissy,
2019-10-23 11:39 - lissy1234567 in Beitrag No. 7 schreibt:
Ich werde in Zukunft drauf achten, hab es diesmal wirklich nicht bemerkt.
kein Problem, aber wenn Du in Zukunft darauf achtest, hilft das allen Beteiligten. Du bist nicht die Einzige, die solche Schwierigkeiten hat, vielleicht helfen Dir einige der Ideen, die auf reproducibleresearch.net beschrieben werden?

2019-10-23 11:39 - lissy1234567 in Beitrag No. 7 schreibt:
Ja, ich sehe die Schwankungen auch...
diese Schwankungen der Gesamtenergie $W(t)$ sind ja bis auf die vertikale Verschiebung um $W(0)$ der Fehler, den ich in Beitrag No. 4 gemeint hatte. Im Vergleich zu der Schrittweite kommt er mir recht groß vor (mehr dazu weiter unten).
2019-10-23 11:39 - lissy1234567 in Beitrag No. 7 schreibt:
im Buch 3x9 Numerik von Sören Bartels habe ich eine Abbildung gefunden, die allerdings auch solche minimale Schwankungen zeigt und darüber eine mathematische Begründung, auf Seite 238 und 239 (weiß nicht, ob ich Screenshots davon hier reinstellen darf).
Die Begründung würde mich schon interessieren, aber ich bin nicht sicher, ob zwei Seiten noch als fair use gelten. Gibt es auch Formeln, die die Abhängigkeit der Schwankung von der Schrittweite $h=T/N$ angeben? Falls ja, könntest Du diese ja nachprüfen.
2019-10-23 11:39 - lissy1234567 in Beitrag No. 7 schreibt:
Beim impliziten Mittelpunktverfahren bspw. ist die Gerade dabei gerade und exakt. Mein Verfahren ist das symplektische Euler-Verfahren, das ich verwendet habe und wie der Name schon sagt ist das symplektisch biggrin
Das setzt voraus, dass Du das symplektische Euler-Verfahren richtig implementiert hast. Ein erster Vergleich mit der Beschreibung auf Wikipedia, Geometric integrator
$\begin{align*}
p_{k+1} &= p_k - h\sin(q_k) \\
q_{k+1} &= q_k + h p_{k+1} \\
\end{align*}$
sieht gut aus.

Zum zweiten Plot in Beitrag No. 5: die Amplitude des Fehlers, also der Differenz der Gesamtenergien der numerischen Lösung und der analytischen Näherung ist etwa doppelt so groß wie die der kinetischen oder potentiellen Energie, die im ersten Plot zu sehen sind. Das sieht für mich nach einem Vorzeichenfehler aus. Wie sieht der Plot von kinetischer und potentieller Energie der analytischen Näherung aus?

2019-10-23 11:39 - lissy1234567 in Beitrag No. 7 schreibt:
Danke für deine gute Hilfe smile
Es freut mich, dass ich Dir helfen konnte. Ich muss selbst noch viel zu diesem Thema lernen.

Servus,
Roland



  Profil  Quote  Link auf diesen Beitrag Link
lissy1234567
Aktiv Letzter Besuch: in der letzten Woche
Dabei seit: 01.09.2017
Mitteilungen: 420
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.9, vom Themenstarter, eingetragen 2019-10-29 13:40


Hallo Roland,

inzwischen hatte ich schon eine Besprechung mit meinem Seminar-Betreuer, da ich das Thema für einen Vortrag erarbeite. Ein paar Vorzeichen-Fehler, wie du es auch schon gesagt hast, waren drin und die Plots sehen nun doch etwas anders aus, allerdings wurde mir bestätigt dass diese Schwankungen 'normal' sind, einfach wegen der numerischen Inexaktheit. Ob es für h = T/N Formeln gibt weiß ich noch nicht - ich werde mal recherchieren :) Danke für den Tipp!

Hier sind nochmal meine 'neuen' Plots - ich enthalte den Code mal vor, falls es dich/jemand anderes interessiert, kann ich ihn aber gerne reinstellen.

Man sieht, dass die Amplitude nun deutlich geringer ist. Grund waren wie gesagt VZ-Fehler sowie einmal das Vertauschen von p und q im Code :D

Hier auch nochmal analytische und numerische Gesamtenerige:


Sieht jetzt erstmal nach großen Schwankungen aus, aber man erkennt ja an der y-Achse, dass es sich nur um minimale Änderugen handelt. Denkst du das ist 'normal'?



  Profil  Quote  Link auf diesen Beitrag Link
rlk
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 16.03.2007
Mitteilungen: 10562
Aus: Wien
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.10, eingetragen 2019-11-04 18:54


Hallo Lissy,
der relative Energiefehler
$$\frac{\Delta W}{W} \approx \frac{2.5\cdot 10^{-4}}{3.4\cdot 10^{-2}} \approx 7.35e-3$$
ist ungefähr fünf mal größer als die relative Schrittweite
$$\frac{h}{T_P}\approx \frac{40}{10000\cdot 3}\approx 1.33e-3$$
die Größenordnung ist also richtig.

Der Energiefehler der analytischen Näherung sollte nur von der Auslenkung, aber nicht von der Schrittweite abhängen.

Servus,
Roland



  Profil  Quote  Link auf diesen Beitrag Link
lissy1234567
Aktiv Letzter Besuch: in der letzten Woche
Dabei seit: 01.09.2017
Mitteilungen: 420
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.11, vom Themenstarter, eingetragen 2019-11-05 12:31


Hallo Roland,

danke für das Vorrrechnen, aber wie kommst du denn auf die Zahlen, also vor allem auf T = 10000*3 (die 3???) und auf die 2.5 im relativen Energiefehler? :o



  Profil  Quote  Link auf diesen Beitrag Link
rlk
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 16.03.2007
Mitteilungen: 10562
Aus: Wien
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.12, eingetragen 2019-11-05 21:51


Hallo Lissy,
das Maximum $\Delta W \approx 2.5\cdot 10^{-4}$ des Energiefehlers habe ich vom ersten Plot in Beitrag  9 abgelesen. Der Wert 3 ist die Periode des Pendels, de ich zuerst $T$ genannt habe. Das war aber ungeschickt, denn in Deinem Programm ist T die Länge des Zeitintervalls, in dem die Lösung berechnet wird. Die Schrittweite ergibt sich aus der Intervalllänge $T$ und der Anzahl $n$ der Punkte:
$$h = \frac{T}{n}$$

Die Rechnung dient nur zur Plausibilitätskontrolle. Hast Du inDeinen Büchern schon Formeln zur Abschätzung des Energiefehlers gefunden?

Servus,
Roland



  Profil  Quote  Link auf diesen Beitrag Link
lissy1234567
Aktiv Letzter Besuch: in der letzten Woche
Dabei seit: 01.09.2017
Mitteilungen: 420
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.13, vom Themenstarter, eingetragen 2019-11-06 11:02


Ah okay, danke - dann ist alles klar :)

ja, ich habe das ein oder andere Theorem zur Abschätzung gefunden, allerdings noch nicht angewandt :D



  Profil  Quote  Link auf diesen Beitrag Link
lissy1234567 hat die Antworten auf ihre/seine Frage gesehen.
lissy1234567 hat selbst das Ok-Häkchen gesetzt.
Neues Thema [Neues Thema]  Druckversion [Druckversion]

 


Wechsel in ein anderes Forum:
 Suchen    
 
All logos and trademarks in this site are property of their respective owner. The comments are property of their posters, all the rest © 2001-2019 by Matroids Matheplanet
This web site was made with PHP-Nuke, a web portal system written in PHP. PHP-Nuke is Free Software released under the GNU/GPL license.
Ich distanziere mich von rechtswidrigen oder anstößigen Inhalten, die sich trotz aufmerksamer Prüfung hinter hier verwendeten Links verbergen mögen.
Lesen Sie die Nutzungsbedingungen, die Distanzierung, die Datenschutzerklärung und das Impressum.
[Seitenanfang]