Nonlinear Equations and Least-Squares Fitting¶

Solve systems of nonlinear equations with np_fsolve and fit overdetermined systems with np_lsq_nonlinear. Both wrap Fortran MINPACK.

Two calling conventions:

  • Expression mode — write equations symbolically; the Jacobian is computed via diff() and compiled automatically
  • Function mode — pass a callback f(x) taking a 1D ndarray; the Jacobian is approximated by finite differences
In [10]:
load("numerics")$
load("numerics-optimize")$
load("ax-plots")$

1. Scalar Root Finding¶

Find $\sqrt{2}$ by solving $x^2 - 2 = 0$.

In [11]:
[x_opt, rnorm, info] : np_fsolve([x^2 - 2], [x], [1.0]);
print("x =", np_ref(x_opt, 0), " residual =", rnorm, " info =", info)$

[ndarray(\#\<ndarray\ 1\ DOUBLE\-FLOAT\ \(1\)\>),8.881784197001252e-16,1]
x = 1.4142135623730947  residual = 8.881784197001252e-16  info = 1

2. System of Equations¶

Find the intersection of the unit circle $x^2 + y^2 = 1$ and the line $y = x$ in the first quadrant.

In [12]:
[xy_opt, rnorm2, info2] : np_fsolve(
  [x^2 + y^2 - 1, y - x], [x, y], [0.5, 0.5]);

print("Intersection:", np_to_list(xy_opt))$
print("Residual:", rnorm2)$

[ndarray(\#\<ndarray\ 2\ DOUBLE\-FLOAT\ \(2\)\>),4.440892098500626e-16,1]
Intersection: [0.7071067811865474,0.7071067811865474]
Residual: 4.440892098500626e-16
In [13]:
/* Visualize */
ax_draw2d(
  color="blue", line_width=2, name="x^2 + y^2 = 1",
  parametric(cos(s), sin(s), s, 0, float(2*%pi)),
  color="orange", line_width=2, name="y = x",
  explicit(x, x, -1.5, 1.5),
  color="red", marker_size=8, name="Solution",
  points([[np_ref(xy_opt, 0), np_ref(xy_opt, 1)]]),
  title="Circle-Line Intersection",
  xlabel="x", ylabel="y", grid=true, showlegend=true
)$
No description has been provided for this image

3. Chemical Equilibrium¶

A reaction $A \rightleftharpoons 2B$ has equilibrium constant $K = [B]^2 / [A]$. Starting from $[A]_0 = 1$, $[B]_0 = 0$, find the equilibrium concentrations (conservation: $[A] + [B]/2 = 1$).

In [14]:
K_eq : 0.5$
[conc, rn, inf_eq] : np_fsolve(
  [cb^2/ca - K_eq, ca + cb/2 - 1],
  [ca, cb], [0.8, 0.4])$

print("[A] =", np_ref(conc, 0))$
print("[B] =", np_ref(conc, 1))$
print("Check K:", np_ref(conc, 1)^2 / np_ref(conc, 0))$

[A] = 0.7034648345913732
[B] = 0.5930703308172536
Check K: 0.5

4. Nonlinear Least Squares: Curve Fitting¶

Fit $y = a \cdot e^{b x}$ to noisy data using np_lsq_nonlinear. We write one residual expression per data point.

In [15]:
/* Generate noisy exponential data: y = 2*exp(0.3*x) + noise */
n_pts : 20$
x_data : np_linspace(0.0, 3.0, n_pts)$
y_true : np_scale(2.0, np_exp(np_scale(0.3, x_data)))$
y_data : np_add(y_true, np_scale(0.15, np_randn([n_pts])))$

/* Build residual list: a*exp(b*x_i) - y_i for each point */
residuals : makelist(
  a_fit * exp(b_fit * np_ref(x_data, i)) - np_ref(y_data, i),
  i, 0, n_pts - 1)$

[ab_opt, res_norm, lsq_info] : np_lsq_nonlinear(
  residuals, [a_fit, b_fit], [1.0, 0.1])$

a_found : np_ref(ab_opt, 0)$
b_found : np_ref(ab_opt, 1)$
print("Fitted: a =", a_found, " b =", b_found)$
print("True:   a = 2, b = 0.3")$

Fitted: a = 1.9740112229040165  b = 0.30184405403101644
True:   a = 2, b = 0.3
In [16]:
ax_draw2d(
  color="red", marker_size=5, name="Noisy data",
  points(np_to_list(x_data), np_to_list(y_data)),
  color="blue", line_width=2, name="Fitted curve",
  explicit(a_found * exp(b_found * x), x, 0, 3),
  color="gray", dash="dash", name="True curve",
  explicit(2 * exp(0.3 * x), x, 0, 3),
  title="Nonlinear Least Squares Fit",
  xlabel="x", ylabel="y", grid=true, showlegend=true
)$
No description has been provided for this image

5. Function Mode¶

When the equations are defined procedurally or involve ndarray operations, use function mode. The callback takes a 1D ndarray and returns a Maxima list.

Expression mode compiles equations and their symbolic Jacobian into native Lisp closures — the Maxima evaluator is not in the solver loop. Function mode calls through the evaluator on each iteration and uses finite-difference Jacobian approximation.

In [17]:
/* np_fsolve function mode: same circle-line problem */
circle_line(w) := [np_ref(w, 0)^2 + np_ref(w, 1)^2 - 1,
                    np_ref(w, 1) - np_ref(w, 0)]$

[xy2, rn2, inf2] : np_fsolve(circle_line, [0.5, 0.5]);
print("Function mode solution:", np_to_list(xy2))$

[ndarray(\#\<ndarray\ 13\ DOUBLE\-FLOAT\ \(2\)\>),4.440892098500626e-16,1]
Function mode solution: [0.7071067811865474,0.7071067811865474]
In [18]:
/* np_lsq_nonlinear function mode */
exp_residuals(p) := block([a_v, b_v],
  a_v : np_ref(p, 0),
  b_v : np_ref(p, 1),
  makelist(
    a_v * exp(b_v * np_ref(x_data, i)) - np_ref(y_data, i),
    i, 0, n_pts - 1))$

[ab2, rn_fn, inf_fn] : np_lsq_nonlinear(exp_residuals, [1.0, 0.1]);
print("Function mode fit:", np_to_list(ab2))$

[ndarray(\#\<ndarray\ 15\ DOUBLE\-FLOAT\ \(2\)\>),0.5070176956707259,1]
Function mode fit: [1.974011223063653,0.30184405401506437]

Summary¶

Function Purpose Expression mode Function mode
np_fsolve Square systems ($n$ equations, $n$ unknowns) Symbolic Jacobian (HYBRJ1) Finite differences (HYBRD1)
np_lsq_nonlinear Overdetermined least squares ($m \ge n$) Symbolic Jacobian (LMDER1) Finite differences (LMDIF1)

Expression mode compiles both the system and its Jacobian into native Lisp closures — the solver loop does not involve the Maxima evaluator. Function mode calls through the evaluator on each iteration.