Compute the Jacobian of a discipline with finite differences

In this example, we will compute the Jacobians of some outputs of a MDODiscipline with respect to some inputs using the finite difference method.

from __future__ import annotations

from gemseo.disciplines.auto_py import AutoPyDiscipline
from numpy import array

First, we create a discipline, e.g. an AutoPyDiscipline:

def f(a=0.0, b=0.0):
    y = a**2 + b
    z = a**3 + b**2
    return y, z

discipline = AutoPyDiscipline(f)

We can execute it with its default input values:

{'a': array([0.]), 'b': array([0.]), 'y': array([0.]), 'z': array([0.])}

or with custom ones:

discipline.execute({"a": array([1.0])})
{'a': array([1.]), 'b': array([0.]), 'y': array([1.]), 'z': array([1.])}

Then, we use the method MDODiscipline.linearize() to compute the derivatives:

jacobian_data = discipline.linearize()

There is no Jacobian data because we need to set the input variables against which to compute the Jacobian of the output ones. For that, we use the method add_differentiated_inputs(). We also need to set these output variables: with the method add_differentiated_outputs(). For instance, we may want to compute the derivative of "z" with respect to "a" only:

jacobian_data = discipline.linearize()
{'z': {'a': array([[1.e-14]])}}

We can have a quick look at the values of these derivatives and verify that they are equal to the analytical results, up to the numerical precision.

By default, GEMSEO uses MDODiscipline.default_inputs as input data for which to compute the Jacobian one. We can change them with input_data:

jacobian_data = discipline.linearize(input_data={"a": array([1.0])})
{'z': {'a': array([[3.0000003]])}}

We can also force the discipline to compute the derivatives of all the outputs with respect to all the inputs:

jacobian_data = discipline.linearize(compute_all_jacobians=True)
{'y': {'a': array([[1.e-07]]), 'b': array([[1.]])}, 'z': {'a': array([[1.e-14]]), 'b': array([[1.e-07]])}}

Lastly, we can change the approximation type to complex step and compare the results:

jacobian_data = discipline.linearize(compute_all_jacobians=True)
/home/docs/checkouts/ ComplexWarning: Casting complex values to real discards the imaginary part
  processed_data[key] = float(val[0])
{'y': {'a': array([[0.]]), 'b': array([[0.]])}, 'z': {'a': array([[0.]]), 'b': array([[0.]])}}

Total running time of the script: ( 0 minutes 0.018 seconds)

Gallery generated by Sphinx-Gallery