Search code examples
pythonnumpyscipylinear-algebradifferential-equations

Ensure trivial solution is found to matrix equation


I'm trying to solve a matrix equation $Ax = b$ with numpy (actually I'm using scipy.sparse but I believe the question remains the same). In my setup, $b$ is the derivative of some function.

In my system, it's possible that the equation may sometimes be $Ax = 0$, in which case I want to get the zero vector back. When I try to do this with a vector of floats, the result I get back is a vector of values ~e-22. This causes havoc at future timesteps and causes weird numerical artefacts when it should be a steady state. When I manually make the right hand side a vector of integer 0s, I get a vector of exactly zeros back.

How can I do this without having to manually check if the vector is all zeros, and convert to integers if so?

Many thanks


Solution

  • I would use np.all to check if every value in the array is below a certain tolerance value, and if it is, convert it to an actual zero-value vector:

    import numpy as np
    
    a = np.array([1e-22, 2e-21])
    print(a)  # [1.e-22, 2.e-21]
    
    tolerance = 1e-20
    if np.all(np.absolute(a) < tolerance):
        a = np.zeros(a.shape)
    
    print(a)  # [0., 0.]