getting a weird error while using solve_ivp from scipy.integrate after qr from numpy.linalg

Question

I am writing a code to calculate Lyapunov exponent for Lorenz system using qr decomposition. I am getting this wired error message which i am unable to figure out. Below piece of code run without any error.

import numpy as np
from scipy.integrate import solve_ivp
from numpy.linalg import qr, norm

def Lorenz(t, x):
    xdot = [s*(x[1]-x[0]), x[0]*(r-x[2])-x[1], x[0]*x[1]-b*x[2]]
    return xdot

s=10; r=28; b = 8/3
y0 = np.random.randn(3)
T = 100
sol = solve_ivp(Lorenz, [0, T], y0, dense_output = True, 
                atol = 1e-9, rtol = 1e-6)

I need to perform a qr decomposition of another numpy array before the final solve_ivp step. So i inserted following lines before the final step.

A = np.random.randn(3, 3)
q, r = qr(A)

Above block of code have no direct reference to any variable or function used in the solve_ivp step but still I am getting error message.

runfile('E:/Codes/ChuaDiode/untitled0.py', wdir='E:/Codes/ChuaDiode')
Traceback (most recent call last):

  File "E:\Codes\ChuaDiode\untitled0.py", line 23, in <module>
    sol = solve_ivp(Lorenz, [0, T], y0, dense_output = True,

  File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\integrate\_ivp\ivp.py", line 542, in solve_ivp
    solver = method(fun, t0, y0, tf, vectorized=vectorized, **options)

  File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\integrate\_ivp\rk.py", line 94, in __init__
    self.f = self.fun(self.t, self.y)

  File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\integrate\_ivp\base.py", line 138, in fun
    return self.fun_single(t, y)

  File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\integrate\_ivp\base.py", line 20, in fun_wrapped
    return np.asarray(fun(t, y), dtype=dtype)

  File "C:\ProgramData\Anaconda3\lib\site-packages\numpy\core\_asarray.py", line 102, in asarray
    return array(a, dtype, copy=False, order=order)

ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (3,) + inhomogeneous part.

I have tried the same code for a 2-dimensional Duffing system which runs smoothly.

def Duffing(t, x):
    xdot = [x[1], x[0] - x[0]**3-c*x[1] + F*np.sin(w*t)]
    return xdot

F = 1.3; w = 0.6; c = 0.2

is this a problem with my python installation. Kindly help. Thanks in advance.

Answer 1

When you write:

def Lorenz(t, x):
    xdot = [s*(x[1]-x[0]), x[0]*(r-x[2])-x[1], x[0]*x[1]-b*x[2]]
    A = np.random.randn(3, 3)
    q, r = qr(A)  # Override the r global variable
    return xdot

You override the r global variable the system use to compute its quantities. Just change for another names:

def Lorenz(t, x):
    xdot = [s*(x[1]-x[0]), x[0]*(r-x[2])-x[1], x[0]*x[1]-b*x[2]]
    A = np.random.randn(3, 3)
    Q, R = qr(A)  # Don't conflict anymore
    return xdot

Even better, don't use global variables at all:

def Lorenz(t, x, s, r, b):
    xdot = [s*(x[1]-x[0]), x[0]*(r-x[2])-x[1], x[0]*x[1]-b*x[2]]
    A = np.random.randn(3, 3)
    Q, R = qr(A)
    return xdot

And pass it to the solver:

sol = solve_ivp(
    Lorenz, [0, T], y0, dense_output=True, args=(s, r, b), 
    atol=1e-9, rtol=1e-6
)

That will solve your issue.