Python Compute Jacobian numerically

Question

Hi I have a class with the purpose to solve system of Differential equations, this class contains a class named Rhs (right hand side) that contains all the feature of the differential problem: init_time , final_time , Initial value , function(that is a numpy.array containing lambda function)

when I have a system of differential equation the functions are managed in this way :

eq1 = lambda t,u : a*(u[0]-u[0]*u[1]);
eq2 = lambda t,u : -c*(u[1]-u[0]*u[1]);

func1 = np.array([eq1,eq2])

y0      = np.array([2.,1.])
system1 = rhs.Rhs(func1, 0,10,y0,500 )

using the parameter above the class Rhs is completed ! it have a method to return the value of each function, and a method for compute the derivate of each of function:

def f(self,ti,ui):
       return  np.array([function(ti,ui) for function in self.func])     


def Df(self,ti,ui):
      eps = 10e-12
      return (self.f(ti,ui+eps) - self.f(ti,ui-eps) )/(2*eps)

Now I have a big problem, in order to create a class to solve an implicit method I have to compute the Jacobian of the function ! but I have no idea how to do this !

EDIT no I need to define the jacobian matrix yes I wrote the method for derivative but I really have not idea how to define J[i][j] @saulspatz may you help me pls ? I don't need to compute it's determinant just to write the terms $J_{i,j}$

I wrote this but doesn't works :

   def J(self,ti,ui):   
     self.t = ti 
     self.u = ui
     Jac = np.zeros((len(ui),len(ui)))
     n = len(ui)
     eps = 1e-12
     for i in range(n):
         for j in range(n):
             Jac[i,j] = self.Df(ti,ui) 

     return Jac

EDIT thank you very much @caverac how can I generalize for any system ? is that straightforward ?

Answer 1

If you have a problem of the form

\begin{eqnarray} \frac{{\rm d}x_1}{{\rm d}t} &=& f_1(x_1, x_2, \cdots) \\ \frac{{\rm d}x_2}{{\rm d}t} &=& f_2(x_1, x_2, \cdots) \\ &\vdots& \end{eqnarray}

then the Jacobian you are looking for is

$$ J = \pmatrix{ \partial f_1/\partial x_1 & \partial f_1/\partial x_2 & \cdots \\ \partial f_2/\partial x_1 & \partial f_2/\partial x_2 & \cdots \\ \vdots & \vdots & \ddots } $$

or

$$ J_{ij} = \frac{\partial f_i}{\partial x_j} $$

in your case

$$ f_1(x_1, x_2) = a(x_1 - x_1 x_2) ~~~\mbox{and}~~ f_2(x_1, x_2) = -c (x_2 - x_1 x_2) $$

So the Jacobian takes the form

$$ J = \pmatrix{ a(1 - x_2) & -ax_1 \\ cx_2 & -c(1 - x_1) } $$

def f(t, x, **params):

    a = params['a']
    c = params['c']

    f1 = a * (x[0] - x[0] * x[1])
    f2 = -c * (x[1] - x[0] * x[1])   

    return np.array([f1, f2], dtype = np.float)

def df(t, x, **params):

    eps = 1e-10
    J = np.zeros([len(x), len(x)], dtype = np.float)

    for i in range(len(x)):
        x1 = x.copy()
        x2 = x.copy()

        x1[i] += eps
        x2[i] -= eps

        f1 = f(t, x1, **params)
        f2 = f(t, x2, **params)

        J[ : , i] = (f1 - f2) / (2 * eps)

    return J

t = 0
x = np.array([1, 2], dtype = np.float)
print df(t, x, a = 1, c = 1)

EDIT Wrapped methods into a class

class ODEModel(object):

    def __init__(self, eps = 1e-10):

        self.eps = eps
        self.farray = []

    def add_function(self, f):

        self.farray.append(f)

    def f(self, t, x):

        return np.array([fi(t, x) for fi in self.farray], dtype = np.float)

    def df(self, t, x):

        J = np.zeros([len(x), len(x)], dtype = np.float)

        for i in range(len(x)):
            x1 = x.copy()
            x2 = x.copy()

            x1[i] += self.eps
            x2[i] -= self.eps

            f1 = self.f(t, x1)
            f2 = self.f(t, x2)

            J[ : , i] = (f1 - f2) / (2 * self.eps)

        return J


F = ODEModel(eps = 1e-12)

eq1 = lambda t,u : u[1]
eq2 = lambda t,u : u[2]
eq3 = lambda t,u : u[3]
eq4 = lambda t,u : -8 * u[0] + np.sin(t) * u[1] - 3 * u[3] + t**2

F.add_function(eq1)
F.add_function(eq2)
F.add_function(eq3)
F.add_function(eq4)


t = 0
x = np.array([1, 2, 3, 4], dtype = np.float)
print F.df(t, x)