code: def d_operator(self, P):
          r
          C, theta = self.C, self.theta
          I = np.identity(self.j)
          S1 = dot(P, C)
          S2 = dot(C.T, S1)
  
          dP = P + dot(S1, solve(theta * I - S2, S1.T))
  
          return dP
docs: r"""
          The D operator, mapping P into
  
          .. math::
  
              D(P) := P + PC(\theta I - C'PC)^{-1} C'P.
  
          Parameters
          ----------
          P : array_like(float, ndim=2)
              A matrix that should be n x n
  
          Returns
          -------
          dP : array_like(float, ndim=2)
              The matrix P after applying the D operator