Lösen nichtlinearer singulärer ODE mit SciPy-Odeint / ODEPACK

8

Ich möchte die Lane-Emden-Isothermengleichung lösen [PDF, Gl. 15.2.9]

d2ψdξ2+2ξdψdξ=eψ

mit den Anfangsbedingungen

ψ(ξ=0)=0dψdξ|ξ=0=0

mit SciPy,odeint() aber wie zu sehen ist, ist die Gleichung am Ursprung singulär. In der Dokumentation wird angegeben, dass ODEPACK verwendet wird.

Ich kenne bereits die Potenzreihen der Lösung in einer Nachbarschaft von ξ=0 ( ref ):

ψ(ξ)ξ26ξ4120+ξ61890

Ich habe versucht Einstellung tcritzu np.array([0.0]), aber nicht Arbeit: Ich habe eine Warnung über ungültige Werte erhalten und dann meine Lösung ist alles NaN. Sollte ich vielleicht ab 0,01 integrieren? Oder gibt es eine andere Lösung?

Astrojuanlu
quelle
t0>0
1
ξt
ξ

Antworten:

7

Gut, diese Antwort ist ein Schuss in die Dunkelheit, aber hier geht es weiter.

Transformieren Sie zunächst die ODE zweiter Ordnung in ein System aus zwei ODEs. Lassen

φ1=ψ,φ2=ψ˙,

ξ

Dann die implizite ODE zweiter Ordnung

ψ¨(ξ)+2ξ1ψ˙(ξ)=eψ(ξ)ψ(0)=0ψ˙(0)=0

kann als explizite ODE erster Ordnung ausgedrückt werden

φ˙1(ξ)=φ2(ξ)φ˙2(ξ)=2ξ1φ2(ξ)+eφ1(ξ)φ1(0)=0φ2(0)=0.

ξ=0ξ0

Erstens wissen wir das

limξ0φ2(ξ)=0,

Da wir davon ausgegangen sind, dass eine Lösung existiert, ist differenzierbar, was bedeutet, dass sie kontinuierlich sein muss. Die Grenze einer stetigen Funktion an einem Punkt ist ihr Wert an diesem Punkt, und wir kennen den Wert von da es sich um eine Anfangsbedingung handelt.φ2φ2(0)

Das wissen wir auch

limξ0eφ1(ξ)=1

aus ähnlichen Gründen; Wir haben angenommen, dass differenzierbar ist, also stetig, und da es sich um eine Anfangsbedingung handelt.φ1φ1(0)=0

Schließlich,

limξ02φ2(ξ)ξ=limξ02φ˙2(ξ),

unter Verwendung der L'Hôpital-Regel auf der unbestimmten Form .0/0

Um weiter fortzufahren, müssen wir eine andere Annahme treffen: ist stetig bei . Dann folgt darausφ˙2ξ=0

limξ02φ˙2(ξ)=2φ2˙(0).

Wenn wir die ODE erster Ordnung erneut betrachten und die rechte Seite bei auswerten , können wir sehen, dass wir:ξ=0

φ˙1(0)=0φ˙2(0)=2φ˙2(0)+1,

woraus folgt, dass .φ˙2(0)=1/3

Mit dieser Analyse können Sie eine ifAnweisung einfügen, die diese Werte der Funktion auf der rechten Seite bei zurückgibt , wodurch Sie die Singularität überwinden sollten. Diese Analyse erfordert jedoch einige Annahmen über die Kontinuität, die möglicherweise zutreffen oder nicht. Nehmen Sie die resultierende Lösung mit einem Salzkorn.ξ=0

Geoff Oxberry
quelle
Geoff, vielleicht du und ? Ich fügte meiner Antwort hinzu, dass ich die Potenzreihen der Lösung am Ursprung bereits kenne (vielleicht hätte ich das vorher tun sollen). Wie auch immer, Sie haben gezeigt, dass , was der Fall ist. Ich werde die Aussage Sache versuchen . φ1=ψφ2=ψ˙ψ¨(0)=1/3if
Astrojuanlu
@ Juanlu001: Guter Anruf; Ich habe den Fehler korrigiert.
Geoff Oxberry
Du hattest Recht! Es war so einfach wie eine einfache ifKlausel. Vielen Dank!
Astrojuanlu
3

Wenn Sie mehr Optionen für Ihren ODE-Solver wünschen, schauen Sie sich das Assimulo- Paket an, das Bindungen an das CVODE-Paket (und RADAU und einige einfache Integratoren) implementiert.

GertVdE
quelle
Sehr wertvolles Paket, wusste es nicht! Vielen Dank.
Astrojuanlu
@GertVdE: Wäre einer dieser Integratoren in der Lage, mit der Singularität umzugehen, ohne eine ifAnweisung zu verwenden, um die richtige Grenze der rechten Seite als auszufüllen ? ξ0
Geoff Oxberry
1
@GeoffOxberry: Ich glaube, der IDA-Löser der Sundials-Suite (den Assimulo für Python umschließt) ermöglicht es, ausgehend von der Vermutung eines Benutzers nach einem konsistenten Anfangswert zu suchen. Dies würde es Juanlu001 ermöglichen, von seiner Serienerweiterung als erste Vermutung auszugehen und IDA nach der richtigen (numerisch, dh) IV auflösen zu lassen.
GertVdE