オイラー法

Question

を使用して、4D結合されたシステムを解く私はコードで以下に与えられたシステムを実装したいが、私は1500回の繰り返しにそれを増やしたときに、私は次のエラーを取得：

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 

私のコード：

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

あなたが実際に（非線形微分方程式の有名敏感なシステムをシミュレートしようとしている、有名な単純な数値スキームでバージョンのa famously sensitive systemをバフにしてください）。あなたの解は、あなたの解決値が最初にinf（あなたが気づいていない）になり、nan（これはまだ気づいていません）に現れた特定の時間ステップで分かれ、最後にAxes3D.plotの倍率はあなたの無限の周りをジャグリングしながらゼロ。そのまま

はここにあなたの出力です：

お知らせ軸はスケールを制限：すべては1e100を超えています。これは、周りに無限を走っているとあなたに語ったことがあります。

良いニュースは、あなたが持っている非常に同じプログラムは、タイムステップを減らすだけで妥当な出力を与えることができるということです。これは、オイラーのような一次メソッドで常に最初の推測でなければなりません。あなたのアルゴリズムが正しいMATLAB実装）。

まず

二つのプロットはかなり合理的、および露骨有限であることに注意：

例（右）t=0.001; itera=25000を使用して（左）とt=0.0001; itera=250000出力します。第2に、2つの軌跡は、概して同様の形状を有するが、非常に異なることに留意されたい。より長い反復を使用すると（これは、合計がt*iteraであることを意味します）、その差は徐々に大きくなり、最終的には2つの軌道が完全に分割されます。

非常に敏感な（正確にはカオス的な）システムの軌道をプロットするには、非常に不正確な方法を使用していることを理解する必要があります。非常に正確な方法であっても、最終的にはいくつかのエラーが蓄積され、初期値問題の実際の解決策から逸脱します。軌道が必然的にジグザグになるアトラクタの大まかな形をプロットすることだけが期待できます。しかし、私はそれがあなたの目標であることはかなり確信しています。