洛倫茲吸引子是如何畫出的？

01-12

Python

代碼如下

# -*- coding: utf-8 -*- from scipy.integrate import odeint import numpy as np


def lorenz(w, t, p, r, b):

    # 給出位置矢量w，和三個參數p, r, b計算出

    # dx/dt, dy/dt, dz/dt的值

    x, y, z = w

    # 直接與lorenz的計算公式對應

    return np.array([p*(y-x), x*(r-z)-y, x*y-b*z]) 
t = np.arange(0, 30, 0.01) # 創建時間點

# 調用ode對lorenz進行求解, 用兩個不同的初始值

track1 = odeint(lorenz, (0.0, 1.00, 0.0), t, args=(10.0, 28.0, 3.0))

track2 = odeint(lorenz, (0.0, 1.01, 0.0), t, args=(10.0, 28.0, 3.0)) 
# 繪圖

from mpl_toolkits.mplot3d import Axes3D

import matplotlib.pyplot as plt

fig = plt.figure() ax = Axes3D(fig) ax.plot(track1[:,0], track1[:,1], track1[:,2]) ax.plot(track2[:,0], track2[:,1], track2[:,2]) plt.show()

互動式三維視圖

代碼

洛倫茲吸引子軌跡

from scipy.integrate import odeint import numpy as np from mayavi import mlab def lorenz(w, t, a, b, c): # 給出位置矢量w，和三個參數a, b, c計算出 # dx/dt, dy/dt, dz/dt的值 x, y, z = w.tolist() # 直接與lorenz的計算公式對應 return np.array([a * (y - x), x * (b - z) - y, x * y - c * z]) t = np.arange(0, 30, 0.01) # 創建時間點 # 調用ode對lorenz進行求解, 用兩個不同的初始值 track1 = odeint(lorenz, (0.0, 1.00, 0.0), t, args=(10.0, 28.0, 3.0)) track2 = odeint(lorenz, (0.0, 1.01, 0.0), t, args=(10.0, 28.0, 3.0)) #繪製圖形 mlab.plot3d(track1[:, 0], track1[:, 1], track1[:, 2], color=(1, 0, 0), tube_radius=0.1) mlab.plot3d(track2[:, 0], track2[:, 1], track2[:, 2], color=(0, 0, 1), tube_radius=0.1)

我寫過一個軟體可以生成洛倫茲圖像,見:YChaos生成混沌圖像

混沌圖像---洛倫茲的蝴蝶

混沌數學之Lorenz(洛倫茨)吸引子

在matlab里用ode45解出x,y,z；然後直接用plot畫圖，只要係數是標準洛倫茲係數就沒有問題；

function dy=lorenzf(t,y)

%%parameter definition%%

p=10;

r=28;

b=(8/3);

dy=zeros(3,1);

dy(1,:)=p*(-y(1)+y(2));

dy(2,:)=(r-y(3))*y(1)-y(2);

dy(3,:)=y(1)*y(2)-b*y(3);

[t,yy]=ode45("lorenzf",[0:0.01:300],[1;1;1]);

x=yy(:,1);

y=yy(:,2);

z=yy(:,3);

figure(1);

plot3(x,y,z);

xlabel("x(t)")

ylabel("y(t)")

zlabel("z(t)")

title("Lorenz")

洛倫茲模型微分方程如下

這只是其中一個例子，描述二維大氣流動。

這個模型隨參數變化而變化，不同的參數有可能不會形成洛倫茲模型，狀態空間向量x的元素小於3也不會形成洛倫茲模型。

洛倫茲模型是指一類對初始值極敏感的模型，所以天氣預報很難預測長時間的天氣，除非初始值無窮精確，但是這個是做不到的，因為數值計算的精度受內存的限制。

同樣的情況也存在於三體或者多體運動，初始值偏差一點長時間的結果就會不一樣。

用Dymola對這個微分方程進行分析，畫出以x1 x2 x3為坐標的時間軌跡就是洛倫茲吸引子。

大家都在用ode，我就用RK4吧

%%%%%%% MATLAB %%%%%%% %%%%%%% main.m %%%%%%% dt = 0.001; t = 0:dt:100; x = zeros(3,length(t));


p = 10;

r = 28;

b = 8/3;
x(:,1) = [-8; 8; r-1];

for tn = 1:1:length(t)-1

    k1 = Lorenz( x(:,tn) , p, r, b);

    k2 = Lorenz( x(:,tn)+dt*k1/2 , p, r, b);

    k3 = Lorenz( x(:,tn)+dt*k2/2 , p, r, b);

    k4 = Lorenz( x(:,tn)+dt*k3, p, r, b);

    x(:,tn+1) = x(:,tn) + dt/6*(k1+2*k2+2*k3+k4);

end
plot3(x(1,:),x(2,:),x(3,:),"-");

box on;

grid on;

drawnow;

%%%%%%% Lorenz.m %%%%%%% function xd = Lorenz( x, p, r, b ) xd = zeros(3,1); xd(1) = - p * x(1) + p * x(2); xd(2) = - x(1) * x(3) + r * x(1) - x(2); xd(3) = x(1) * x(2) - b * x(3); end

洛倫茲吸引子wiki有介紹，連matlab、mathematica、python代碼都貼出來了。

Lorenz system - Wikipedia

前面已貼出Matlab和Python,我就貼一下Mathematica代碼：

Mathematica代碼

sigma = 3; rho = 26.5; beta = 1; soln = {x[t], y[t], z[t]} /. NDSolve[{x"[t] == sigma (y[t] - x[t]), y"[t] == rho x[t] - x[t]z[t] - y[t], z"[t] == x[t] y[t] - beta z[t], x[0] == 0, y[0] == 1, z[0] == 1}, {x[t], y[t], z[t]}, {t, 0, 100}, MaxSteps -&> 10000][[1]]; ParametricPlot3D[soln, {t, 0, 100}, PlotPoints -&> 10000, Compiled -&> False]

以下視頻是wiki上的，原址Lorenz，只是感覺很酷炫。

Plus

也可以用Pov-Ray進行渲染，這是臨時找到的代碼，原址Lorenz Attractor

//****************************************************************************************** //* Lorenz Attractor, by Paul Nylander, bugman123.com, 7/15/07 //******************************************************************************************


camera{location 34*&<1,-1,1&> look_at &<-2,0,26.5&> up z right x*image_width/image_height sky &<0,0,1&>}

light_source{34*&<1,-1,1&>,1}

#declare rtube=0.25; #declare dt=0.01; #declare sigma=3; #declare rho=26.5; #declare beta=1; #macro dpdt(p) & #end #declare p1=&<0,1,1&>; #declare T=0; #while(T&<100) #declare k1=dt*dpdt(p1); #declare k2=dt*dpdt(p1+k1/2); #declare k3=dt*dpdt(p1+k2/2); #declare k4=dt*dpdt(p1+k3); #declare p2=p1+(k1+2*k2+2*k3+k4)/6; // 4th order Runge-Kutta method union{sphere{p2,rtube} #if(vlength(p2-p1)&>0.001) cylinder{p1,p2,rtube} #end pigment{rgb 1} finish{phong 0.8 phong_size 1.5} } #declare p1=p2; #declare T=T+dt; #end