python--实现Lorenz 63模式

通过python实现Lorenz 63模式(Lorenz,1963)：

其中 σ , ρ 和 β 是参数，分别设置为10,28,8/3。而x,y,z是模式状态变量，在模式中记为矢量 x→ 的三个元素 x1 , x2 , x3 可以使用数值方法求解常微分方程(欧拉法或者Runge-Kutta都可以，这里使用RK45方法)，从初值求得步长 δt 后的x状态：

def RK45(x,func,h):K1=func(x);K2=func(x+h/2*K1);K3=func(x+h/2*K2);K4=func(x+h*K3);x1=x+h/6*(K1+2*K2+2*K3+K4);return x1def L63_rhs(x):# ODE右端项import numpy as npdx=np.ones_like(x);sigma=10.0; rho=28.0;beta=8/3;   # default parametersdx[0]=sigma*(x[1]-x[0])dx[1]=rho*x[0]-x[1]-x[0]*x[2];dx[2]=x[0]*x[1]-beta*x[2]return dxdef L63_adv_1step(x0,delta_t):# 使用RK45求解ODE，从初值x0求得步长delta_t后的x状态x1=RK45(x0,L63_rhs,delta_t)return x1

调用 L63_adv_1step 可以积分模式，测试以x0为初值积分5000步，图像如下，称为洛伦茨吸引子（蝴蝶效应）

import numpy as np
# 模式积分
x0 = np.array([1.508870, -1.531271, 25.46091])
Xtrue = np.zeros([3000,3]);Xtrue[0]=x0
delta_t=0.01
for j in range(1,3000):Xtrue[j] = L63_adv_1step(Xtrue[j-1],delta_t)
# 画图
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
mpl.rcParams['legend.fontsize'] = 10
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(Xtrue[:,0], Xtrue[:,1], Xtrue[:,2],'r', label='Lorenz 63 model')
ax.legend()
plt.xlabel('x');plt.ylabel('y');
ax.set_zlabel('z')

Lorenz63模式具有强非线性，即使初值进行微小的扰动，也能对积分的结果造成巨大影响

# 模式积分
x0p = x0+0.001
Xctl = np.zeros([3000,3]);Xctl[0]=x0p
for j in range(1,3000):Xctl[j] = L63_adv_1step(Xctl[j-1],delta_t)
# 画图部分
fig = plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.plot(Xtrue[range(1000),1], Xtrue[range(1000),2],'r', label='Truth')
plt.plot(Xtrue[0,1], Xtrue[0,2],'bx',ms=10,mew=3)
plt.plot(Xtrue[1000,1], Xtrue[1000,2],'bo',ms=10)
plt.ylim(0,50);plt.title('True',fontsize=15);plt.ylabel('z');plt.xlabel('y')
plt.text(5,25,r'$x_0$',fontsize=14)
plt.grid()
plt.subplot(1,2,2)
plt.plot(Xctl[range(1000),1], Xctl[range(1000),2],'g')
plt.plot(Xctl[0,1], Xctl[0,2],'bx',ms=10,mew=3,label='t0')
plt.plot(Xctl[1000,1], Xctl[1000,2],'bo',ms=10,label='t')
plt.ylim(0,50);plt.title('Control',fontsize=15);plt.ylabel('z');plt.xlabel('y')
plt.grid();plt.legend()
plt.text(5,25,r'$x_0+0.001$',fontsize=14)

"X"代表初始状态，左右两边的初值只相差0.001。模式积分不久，两边的差异就非常大了。正是在于模式具有这种混沌性质，微小的初始偏差能够造成巨大的预报误差，需要使用同化手段利用现实的观测纠正模式，或者提供更精准的初值。

https://dafi.readthedocs.io/en/latest/tutorial_lorenz.html

python--实现Lorenz 63模式相关推荐

USF MSDS501 计算数据科学中文讲义 2.4 Python 中的编程模式
来源:ApacheCN『USF MSDS501 计算数据科学中文讲义』翻译项目原文:Programming Patterns in Python 译者:飞龙协议:CC BY-NC-SA 4.0 现 ...
Python设计模式-装饰器模式
Python设计模式-装饰器模式代码基于3.5.2,代码如下; #coding:utf-8 #装饰器模式class Beverage():name = ""price = 0.0 ...
Python设计模式-中介者模式
Python设计模式-中介者模式代码基于3.5.2,代码如下; #coding:utf-8 #中介者模式class colleague():mediator = Nonedef __init__(s ...
Python设计模式-职责链模式
Python设计模式-职责链模式代码基于3.5.2,代码如下; #coding:utf-8 #职责链模式class Handler():def __init__(self):self.success ...
Python设计模式-享元模式
Python设计模式-享元模式基于Python3.5.2,代码如下 #coding:utf-8class Coffee:name = ""price = 0def __init_ ...
python的编程模式-Python设计模式之状态模式原理与用法详解
本文实例讲述了Python设计模式之状态模式原理与用法.分享给大家供大家参考,具体如下: 状态模式(State Pattern):当一个对象的内在状态改变时允许改变其行为,这个对象看起来像是改变了其类 ...
python策略模式包含角色_详解Python设计模式之策略模式
虽然设计模式与语言无关,但这并不意味着每一个模式都能在每一门语言中使用.<设计模式:可复用面向对象软件的基础>一书中有 23 个模式,其中有 16 个在动态语言中"不见了,或者简 ...
python文件合法模式组合_以下选项中,不是Python文件二进制打开模式的合法组合是...
以下选项中,不是Python文件二进制打开模式的合法组合是答:\"x+\" 建立良好的谈判气氛主要是在( )阶段答:开局 the ruling class had long b ...
[转载] Python中的Phyllotaxis模式| 算法植物学的一个单位
参考链接: Python中的Phyllotaxis模式| 算法植物学的单位简介| 叶底 Phyllotaxis / phyllotaxy是植物茎上叶子的排列,Phyllotactic螺旋形成自然 ...

python--实现Lorenz 63模式

python--实现Lorenz 63模式相关推荐

最新文章

热门文章