Lorenz-63

import numpy as np

class L63():
    def __init__(self):
        self.p = 10.0 # prandtl number
        self.r = 32.0 # rayleigh number
        self.b = 8.0/3.0 # aspect ratio
        self.dt = 0.01

    def __call__(self,w):
        x, y, z = w 
        dw = np.zeros_like(w)
        dw[0] = -self.p * x + self.p * y
        dw[1] = -x * z      + self.r * x - y
        dw[2] = x * y       - self.b * z
        return w + dw * self.dt
    
    def tlm(self,wb,wt):
        xb, yb, zb = wb
        xt, yt, zt = wt 
        dw = np.zeros_like(wt)
        dw[0] = -self.p * xt     + self.p * yt
        dw[1] = (self.r - zb)*xt - yt          - xb*zt
        dw[2] = yb*xt            + xb*yt       - self.b * zt
        return wt + dw * self.dt
    
    def adj(self,wb,wt):
        xb, yb, zb = wb
        xt, yt, zt = wt 
        dw = np.zeros_like(wt)
        dw[0] = -self.p * xt + (self.r - zb)*yt + yb*zt
        dw[1] =  self.p * xt - yt               + xb*zt
        dw[2] =              - xb*yt            - self.b * zt
        return wt + dw * self.dt

import matplotlib.pyplot as plt 

model = L63()
nstep = 500

w1 = np.zeros((nstep+1,3))
w2 = np.zeros((nstep+1,3))
w1[0,:] = 1.0,3.0,5.0
w2[0,:] = 1.1,3.3,5.5
for i in range(nstep):
    w1[i+1,] = model(w1[i])
    w2[i+1,] = model(w2[i])

fig = plt.figure(figsize=[8,8])
ax = fig.add_subplot(projection="3d")
ax.plot(*w1.transpose())
ax.plot(*w2.transpose())
plt.show()

TLM and ADJ check

model = L63()
# TLM check
w0 = np.random.randn(3)
dw0 = np.random.randn(3)*0.1
alp = 1.0e-5
wp = model(w0+alp*dw0)
wb = model(w0)
dw = model.tlm(w0,dw0)
ratio = np.sqrt(np.dot((wp-wb),(wp-wb)))/alp/np.sqrt(np.dot(dw,dw))
print(f"|M(x+a*dx)-M(x)|/|a*TLM*dx|={ratio}")
# ADJ check
MtMw = model.adj(w0,dw)
print(f"(Mx)^T(Mx)-x^TM^TMx={np.dot(dw,dw)-np.dot(dw0,MtMw)}")

|M(x+a*dx)-M(x)|/|a*TLM*dx|=1.000000001770137
(Mx)^T(Mx)-x^TM^TMx=-3.469446951953614e-18

Singular vector

## base field
wb = np.array([1.0,3.0,5.0])
nstep = 100
for i in range(nstep):
    wb = model(wb)

Lanczos method

# initialize
alpha=1.0e-5
dw = np.random.randn(3)
dw = dw * alpha / np.sqrt(np.dot(dw,dw))
niter=0
while(True):
    dwp = dw.copy()
    # forward integration
    dw = model.tlm(wb,dw)
    # backward integration
    dw = model.adj(wb,dw)
    # rescaling
    dw = dw * alpha / np.sqrt(np.dot(dw,dw))
    # convergence evaluation
    diff = np.sqrt(np.dot((dw-dwp),(dw-dwp)))
    niter+=1
    print(f"iter:{niter}, diff={diff:.3e}")
    if diff/alpha<1.0e-6: break
print(dw/alpha)

iter:1, diff=1.249e-06
iter:2, diff=1.128e-06
iter:3, diff=1.068e-06
iter:4, diff=1.024e-06
iter:5, diff=9.678e-07
iter:6, diff=8.925e-07
iter:7, diff=8.016e-07
iter:8, diff=7.032e-07
iter:9, diff=6.054e-07
iter:10, diff=5.137e-07
iter:11, diff=4.314e-07
iter:12, diff=3.596e-07
iter:13, diff=2.981e-07
iter:14, diff=2.463e-07
iter:15, diff=2.030e-07
iter:16, diff=1.670e-07
iter:17, diff=1.372e-07
iter:18, diff=1.127e-07
iter:19, diff=9.247e-08
iter:20, diff=7.587e-08
iter:21, diff=6.223e-08
iter:22, diff=5.103e-08
iter:23, diff=4.185e-08
iter:24, diff=3.431e-08
iter:25, diff=2.813e-08
iter:26, diff=2.307e-08
iter:27, diff=1.891e-08
iter:28, diff=1.551e-08
iter:29, diff=1.271e-08
iter:30, diff=1.042e-08
iter:31, diff=8.544e-09
iter:32, diff=7.005e-09
iter:33, diff=5.743e-09
iter:34, diff=4.708e-09
iter:35, diff=3.860e-09
iter:36, diff=3.165e-09
iter:37, diff=2.594e-09
iter:38, diff=2.127e-09
iter:39, diff=1.744e-09
iter:40, diff=1.430e-09
iter:41, diff=1.172e-09
iter:42, diff=9.609e-10
iter:43, diff=7.878e-10
iter:44, diff=6.458e-10
iter:45, diff=5.295e-10
iter:46, diff=4.341e-10
iter:47, diff=3.559e-10
iter:48, diff=2.918e-10
iter:49, diff=2.392e-10
iter:50, diff=1.961e-10
iter:51, diff=1.608e-10
iter:52, diff=1.318e-10
iter:53, diff=1.081e-10
iter:54, diff=8.859e-11
iter:55, diff=7.263e-11
iter:56, diff=5.954e-11
iter:57, diff=4.881e-11
iter:58, diff=4.002e-11
iter:59, diff=3.281e-11
iter:60, diff=2.690e-11
iter:61, diff=2.205e-11
iter:62, diff=1.808e-11
iter:63, diff=1.482e-11
iter:64, diff=1.215e-11
iter:65, diff=9.962e-12
[ 0.59253215  0.76624248 -0.24855202]