archimedes.implicit¶

archimedes.implicit(func: Callable, static_argnames: str | Sequence[str] | None = None, solver: str = 'newton', name: str | None = None, **options) → FunctionCache¶

Construct an explicit function from an implicit relation.

Given a relation F(x, p) = 0 that implicitly defines x as a function of p, this function creates a solver that computes x given p and an initial guess. This effectively transforms the equation F(x, p) = 0 into an explicit function x = G(x0, p) that returns the solution x for parameters p starting from initial guess x0.

Parameters:

func (callable) – The implicit function, with signature func(x, *args). Must return a residual with the same shape and dtype as the input x.
static_argnames (tuple of str, optional) – Names of arguments that should be treated as static (non-symbolic) parameters. Static arguments are not differentiated through and the solver will be recompiled when their values change.
solver (str, default="newton") –
The root-finding method to use. Options are:
- "newton" : Newton’s method (default), best for general problems
- "fast_newton" : Simple Newton iterations with no line search
- "kinsol" : KINSOL solver from SUNDIALS, robust for large systems
name (str, optional) – Name for the returned function. If None, a name will be generated based on the input function name.
tol (float, optional) – Absolute for convergence. If None, the default tolerance for the chosen method will be used.
**options (dict, optional) –
Common additional options specific to the chosen method:

For "newton" and "fast_newton":
- max_iter : int, maximum iterations (default: 100)
For "kinsol":
- max_iter : int, maximum iterations
- strategy : str, globalization strategy ("none", "linesearch", "picard", "fp")
See the CasADi documentation for more details on the available options for each method.

Returns:

solver – A function that, when called with signature solver(x0, *args), returns the solution x to F(x, *args) = 0 starting from the initial guess x0. This function can be evaluated both numerically and symbolically.

Return type:

FunctionCache

Notes

When to use this function:

When you have an equation F(x, p) = 0 that you need to solve repeatedly for different values of parameters p
For implementing equations of motion for constrained mechanical systems
For implicit numerical schemes in simulation
For embedding root-finding operations within larger computational graphs

The solver automatically computes the Jacobian of the residual function using automatic differentiation.

For simple one-off root-finding problems where parameters don’t change, consider using root() instead which provides a simpler interface.

The function generates a callable that behaves as a pure function, making it suitable for embedding in larger computational graphs or differentiating through.

Examples

>>> import numpy as np
>>> import archimedes as arc
>>>
>>> # Example 1: Simple nonlinear equation x^2 = p
>>> def f(x, p):
...     return x**2 - p
>>>
>>> # Create a solver for x given p
>>> sqrt = arc.implicit(f)
>>>
>>> # Solve for sqrt(2) starting from initial guess 1.0
>>> x = sqrt(1.0, 2.0)
>>> print(f"Solution: {x:.10f}")  # Should be close to 1.4142135624
Solution: 1.4142135624
>>>
>>> # Example 2: Implicit equation with vector input and parameters
>>> def nonlinear_system(x, a, b):
...     # System: a*x[0]^2 + b*x[1] = 1, x[0] + x[1]^2 = 4
...     return np.array([
...         a * x[0]**2 + b * x[1] - 1,
...         x[0] + x[1]**2 - 4
...     ], like=x)
>>>
>>> # Create a solver with static parameter 'a'
>>> solver = arc.implicit(nonlinear_system, static_argnames=('a',))
>>>
>>> # Solve the system for a=2, b=3
>>> x0 = np.array([0.0, 0.0])  # Initial guess
>>> solution = solver(x0, 2.0, 3.0)
>>> print(solution)  # Should converge to a valid solution
>>>
>>> # Example 3: Using a different solver and options
>>> def kepler(x, e):
...     # Kepler's equation: x - e*sin(x) = M where M is fixed at 0.5
...     return x - e * np.sin(x) - 0.5
>>>
>>> # Create solver with more iterations and higher accuracy
>>> kepler_solver = arc.implicit(
...     kepler,
...     solver="newton",
...     max_iter=50,
...     tol=1e-12
... )
>>>
>>> # Solve for eccentric anomaly with eccentricity e=0.8
>>> anomaly = kepler_solver(0.5, 0.8)
>>> print(f"Eccentric anomaly: {anomaly:.10f}")