带jit码的lyapunov指数

Question

我正在使用JIT代码来计算洛伦兹振子的李雅普诺夫指数

下面是documentation后面的简单脚本：

import numpy as np
import pylab as plt
from jitcode import jitcode_lyap, y


p, r, b = 10.0, 28.0, 8.0/3.0

x0 = np.array([10.0, 10.0, 5.0])
f = [
    p * y(1) - p * y(0),
    -y(0) * y(2) + r * y(0) - y(1),
    y(0) * y(1) - b * y(2)
]

n = len(f)
times = range(10, 1000, 10)
nstep = len(times)
lyaps = np.zeros((nstep, n))

ODE = jitcode_lyap(f, n_lyap=n)
ODE.set_integrator("vode")
ODE.set_initial_value(x0, 0.0)

for time in times:
    print(ODE.integrate(time)[1])

我得到了以下错误

Generating, compiling, and loading C code.
generated C code for f
generated symbolic Jacobian
generated C code for Jacobian
/usr/local/lib/python3.6/dist-packages/scipy/integrate/_ode.py:1009: UserWarning: vode: Excess work done on this call. (Perhaps wrong MF.)
  self.messages.get(istate, unexpected_istate_msg)))
Traceback (most recent call last):
  File "main.py", line 29, in <module>
    print(ODE.integrate(time)[1])
  File "/home/abolfazl/.local/lib/python3.6/site-packages/jitcode/_jitcode.py", line 755, in integrate
    super(jitcode_lyap, self).integrate(*args, **kwargs)
  File "/home/abolfazl/.local/lib/python3.6/site-packages/jitcode/_jitcode.py", line 656, in integrate
    return self.integrator.integrate(*args,**kwargs)
  File "/home/abolfazl/.local/lib/python3.6/site-packages/jitcode/integrator_tools.py", line 131, in integrate
    raise UnsuccessfulIntegration
jitcode.integrator_tools.UnsuccessfulIntegration

我认为这个方程不是问题，因为我已经用jitcode解决了这个问题，但是jitcode_lyap无法解决它。 谢谢你的指导

Answer 1

我从开发商那里得到了答案：

您的问题是集成步骤太长。由于切线向量仅在每个积分步骤后重新归一化，这会导致数值溢出，因为洛伦兹系统的时间尺度比文档中用作示例的时间尺度小得多

使用times = range(1, 100, 1)可以解决这个问题