Skip to content

fix(indLin): restore the inductive-linearization iteration #1185

Description

@mattfidler

method="indLin" currently performs one relinearization per hmax-capped
substep. The fixed-point iteration that makes it inductive linearization was
deleted.

Current behavior

Commit bf85a053e ("add indLin method", 2026-06-09) removed case 3: and
case 4: from switch(doIndLin) in src/expm.cpp -- 59 lines. Recover them
with git show bf85a053e^:src/expm.cpp.

Their structure:

  • y0 holds the interval-start state and is never updated; w is the iterate.
  • Each pass calls meOnly(cSub, w.memptr(), y0.memptr(), ...). meOnly passes
    its yc_ argument to ME as __zzStateVar__ (src/expm.cpp:251), so the
    matrix is built at the latest iterate while propagation always starts from
    y0.
  • Loop bound maxsteps = op->mxstep; convergence when
    fabs(w[j]-wLast[j]) < rtol[j]*fabs(w[j]) + atol[j] for every j in
    op->indLin[].

rtol, atol and maxsteps are still declared at src/expm.cpp:329-331 and
are unused in the live body. The current default: branch is
stop(_("unsupported indLin code: %d"), doIndLin), so codes 3/4 error rather
than degrade.

Change

  • Reinstate case 3:/case 4: nested inside the hmax substep loop at
    src/expm.cpp:364-396: subdivide [tp,tf] by hmax, iterate to convergence
    within each substep. Pass the state pointer to IndF.
  • src/genModelVars.c -- flip fullIndLin to length(wIndLin) > 0.
  • src/rxData.cpp -- no change.
  • Surface non-convergence. The deleted code returned 1 unconditionally after
    exhausting maxsteps.

Correct value after this issue

A Michaelis-Menten model written with an indLin() forcing converges to the
LSODA answer. Linear models and state-free forcings keep codes 1/2 and are
bit-identical to before.

Test

In tests/testthat/test-ind-lin.R:

  • MM via an indLin() forcing matches liblsoda to 1e-5 at the default hmax
  • a linear matExp() model is bit-identical to before the change
  • hmax still refines the answer at iteration convergence (this pins the order
    of convergence, and is what makes the deferred linear-ramp issue measurable)
  • a deliberately non-converging model reports rather than silently returning the
    last iterate

Blocked by #1183 and #1184.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions