# SciPy odeint giving different results with matrix multiplication

I've asked this at stackoverflow but maybe this community will have a better idea of the answer.

I'm currently trying to develop a function that performs matrix multiplication while expanding a differential equation with odeint in Python and am seeing strange results.

I have the below matrix of values and function that takes specific values of that matrix:

from scipy.integrate import odeint
x0_train = [2,0]
dt = 0.01
t = np.arange(0, 1000, dt)
matrix_a = np.array([-0.09999975, 1.999999, -1.999999, -0.09999974])
# Function to run odeint with
def f(x, t, a):
return [
a * x + a * x,
a * x - a * x
]
odeint(f, x0_train, t, args=(matrix_a,))

>>> array([[ 2.        ,  0.        ],
[ 1.99760115, -0.03999731],
[ 1.99440529, -0.07997867],
...,
[ 1.69090227,  1.15608741],
[ 1.71199436,  1.12319701],
[ 1.73240339,  1.08985846]])


This seems right, but when I create my own function to perform multiplication/regression, I see the results at the bottom of the array are completely different. I have two sparse arrays that provide the same conditions as matrix_a but with zeros around them.

from sklearn.preprocessing import PolynomialFeatures
new_matrix_a = array([[ 0.        , -0.09999975,  1.999999  ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ],
[ 0.        , -1.999999  , -0.09999974,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
0.        ]])
# New function
def f_new(x, t, parameters):
polynomials = PolynomialFeatures(degree=5)
x = np.array(x).reshape(-1,2)
#x0_train_array_reshape = x0_train_array.reshape(1,2)
polynomial_transform = polynomials.fit(x)
polynomial_features = polynomial_transform.fit_transform(x).T
x_ode = np.matmul(parameters,polynomial_features)
y_ode = np.matmul(parameters,polynomial_features)
return np.concatenate((x_ode, y_ode), axis=None).tolist()

odeint(f_new, x0_train, t, args=(new_matrix_a,))

>>> array([[ 2.00000000e+00,  0.00000000e+00],
[ 1.99760142e+00, -3.99573216e-02],
[ 1.99440742e+00, -7.98188169e-02],
...,
[-3.50784051e-21, -9.99729456e-22],
[-3.50782881e-21, -9.99726119e-22],
[-3.50781711e-21, -9.99722781e-22]])


As you can see, I'm getting completely different values at the end of the array. I've been running through my code and can't seem to find a reason why they would be different. Does anybody have a clear reason why or if I'm doing something wrong with my f_new? Ideally, I'd like to develop a function that can take any values in that matrix_a, which is why I'm trying to create this new function.

• Not an answer, but try just sampling your two functions with random inputs. While f and f_new seem to match perfectly for any value of x0_train, if you change x0_train to anything other than 0, the second output differs significantly. It looks like the difference is -0.2*x0_train, though I haven't dug into your second function enough to see why that would happen. Nov 19, 2021 at 16:48
• @Tyberius strange, I've written a new function for f that performs multiplication rather than how it appears above, so now it's more similar to f_new, and I get similar results to f_new.
– AW27
Nov 19, 2021 at 16:52

I focused on the second function initially, assuming the first was correct, but I realize the issue is actually in the first function.

Your matrix multiplication is incorrect in the first function. It should look like this:

def f(x, t, a):
return [
a * x + a * x,
a * x + a * x
]


You were subtracting a * x rather than adding it. This was exactly the -0.2*x0_train I pointed out in the comments; since $$a_3\approx -0.1$$, and you subtracted it rather than adding it, that would result in the second component of the output being off by 2*a*x at each step. Making this change, both produce the same result when passed into odeint.

• Thanks, I was wondering if that might be an issue. Do you know why the values are incredibly small as the ode increases?
– AW27
Nov 19, 2021 at 17:05
• @AW27 That should be easy to solve. Your system is $\dot x=Ax$ with $A=\pmatrix{-0.1&2\\-2&-0.1}$. This has eigenvalues $-0.1\pm 2i$, as the matrix has the structure of a complex number representation. Thus at the end you get values of size $e^{-0.1·1000}=e^{-100}=3.72·10^{-44}$. Apparently you get interference from floating point errors or the odeint-internal step size control so that the observed values do not sink that low. Nov 19, 2021 at 17:19
• @AW27 Lutz has the right idea here. Since you are testing a fairly simple case here, always see if you can solve an analytical version of the problem to use as reference. In general, even if you can't get an analytical solution, you should try to get a sense of roughly what the solution should be (or at least some properties of it) before calculating it numerically. Nov 19, 2021 at 17:26