Lösen von 4D gekoppelten Systems von Eulers Verfahren mit

Question

Ich möchte das System implementieren unten in Code gegeben, aber wenn ich es bis 1500 Wiederholungen erhöhen, dann bekomme ich folgende Fehler:

Warning (from warnings module): 
    File "D:\python test files\sys1.py", line 16 
    dy = c*x- x*z + w 
RuntimeWarning: overflow encountered in double_scalars 

Warning (from warnings module): 
    File "D:\python test files\sys1.py", line 17 
    dz = -b*z + x*y 
RuntimeWarning: overflow encountered in double_scalars 

Warning (from warnings module): 
    File "D:\python test files\sys1.py", line 18 
    du = -h*u - x*z 
RuntimeWarning: overflow encountered in double_scalars 

Warning (from warnings module): 
    File "D:\python test files\sys1.py", line 42 
    zs[i+1] = zs[i] + (dz * t) 
RuntimeWarning: invalid value encountered in double_scalars 

Warning (from warnings module): 
    File "D:\python test files\sys1.py", line 15 
    dx = a*(y-x) + u 
RuntimeWarning: invalid value encountered in double_scalars 

Warning (from warnings module): 
    File "D:\python test files\sys1.py", line 19 
    dw = k1*x - k2*y 
RuntimeWarning: invalid value encountered in double_scalars 

Warning (from warnings module): 
    File "C:\Python27\lib\site-packages\mpl_toolkits\mplot3d\proj3d.py", line 156 
    txs, tys, tzs = vecw[0]/w, vecw[1]/w, vecw[2]/w 
RuntimeWarning: invalid value encountered in divide 

Mein Code:

from __future__ import division 
import numpy as np  
import math  
import random 
import matplotlib.pyplot as plt  
from mpl_toolkits.mplot3d import Axes3D   
# import pdb 
# pdb.set_trace() 

def sys1(x, y, z, u, w , a=10, b=8.0/3.0, c=28, k1=0.4, k2=8, h=-2): 

    dx = a*(y-x) + u 
    dy = c*x- x*z + w  
    dz = -b*z + x*y  
    du = -h*u - x*z  
    dw = k1*x - k2*y  
    return dx, dy, dz, du, dw 

t = 0.01  
itera = 2500 

# Need one more for the initial values 

xs = np.empty((itera+1,))  
ys = np.empty((itera+1,))  
zs = np.empty((itera+1,))  
us = np.empty((itera+1,))  
ws = np.empty((itera+1,)) 

# Setting initial values 

xs[0], ys[0], zs[0], us[0], ws[0] = (0.1, 0.1, 0.1, 0.1, 0.1) 


# Stepping through "time". 

for i in range(itera):  
# Derivatives of the X, Y, Z state 
    dx, dy, dz, du, dw = sys1(xs[i], ys[i], zs[i], us[i], ws[i])  

    xs[i+1] = xs[i] + (dx * t) 
    ys[i+1] = ys[i] + (dy * t) 
    zs[i+1] = zs[i] + (dz * t) 
    us[i+1] = us[i] + (du * t) 
    ws[i+1] = ws[i] + (dw * t) 

fig = plt.figure() 

ax = fig.gca(projection='3d') 
ax.plot(xs, ys, zs) 
ax.set_xlabel("X Axis") 
ax.set_ylabel("Y Axis") 
ax.set_zlabel("Z Axis") 
ax.set_title("Lorenz Attractor") 
plt.show() 

Answer 1

Sie versuchen, ein berühmten sensibles System von nichtlinearen Differentialgleichungen zu simulieren (eigentlich ein aufgerauht Version von a famously sensitive system) mit einem berühmten einfachen numerischen Schema. Ihre Lösung divergiert zu einem bestimmten Zeitschritt, der sich in Ihren Lösungswerten manifestiert, die zuerst inf werden (was Sie nicht bemerken), dann nan (was Sie immer noch nicht bemerken), dann erzeugt die Skalierung in Axes3D.plot eine Division durch Null, während du mit deinen Unendlichkeiten jonglierst.

Hier ist die Ausgabe, wie sie ist:

Beachten Sie die Achsen begrenzen Waage: alle oben 1e100 sind. Das könnte dir gesagt haben, dass du Unendlichkeiten hast, die herumlaufen.

Die gute Nachricht ist, dass das gleiche Programm Sie vernünftige Ausgabe geben können, nur durch die Verringerung der Zeit Schritt, die immer Ihre erste Schätzung mit einer Methode erster Ordnung wie Euler sein sollte (vor allem, dass Sie von einem überzeugt sind MATLAB-Implementierung, dass Ihr Algorithmus korrekt ist).

Beispiel Ausgänge mit t=0.001; itera=25000 (links) und t=0.0001; itera=250000 (rechts):

Erstens beachten Sie, dass die beiden Parzellen ziemlich vernünftig sind, und offensichtlich endlich. Zweitens, beachten Sie, dass die beiden Bahnen, obwohl sie eine allgemein ähnliche Form haben, sehr unterschiedlich sind. Wenn Sie längere Iterationen verwenden würden (mit denen ich insgesamt länger meine, t*itera), würden die Unterschiede allmählich pronouned werden, und schließlich würden sich die zwei Trajektorien vollständig aufteilen.

Sie sollten sicherstellen, dass Sie verstehen, dass Sie eine sehr ungenaue Methode verwenden, um die Trajektorie eines sehr sensiblen (um genau zu sein: chaotischen) Systems zu zeichnen. Selbst mit einer sehr genauen Methode werden Sie möglicherweise einige Fehler sammeln, und Sie werden von der tatsächlichen Lösung zu Ihrem Anfangswertproblem abweichen. Alles, worauf Sie hoffen können, ist die grobe Form des Attraktors, um den herum die Bahnen im Zickzack verlaufen. Aber ich bin mir ziemlich sicher, dass das dein Ziel ist.