Nonlinear ZDF SVF part 2
Previously we derived a nonlinear zero delay feedback filter where the nonlinear tanh() elements were placed at the inputs of the integrators. This is not ideal, since it will lead to instabilities at high input amplitudes even with implicit integration.
Another way to model nonlinearities present in active analog integrators is to place the tanh() clipping at the “output” of the integrator, where it will limit the system to the $[1, 1]$ interval.
Nonlinear “derivative”
Expressing tanh()limited integration as a differential equation is not really possible with traditional methods since it breaks the fundamental theorem of calculus and the definition of the derivative as a linear difference quotient.
However, it is easy to express as an equation of discrete integration where for example
taken as the linear slope approximation (Euler integration)
becomes
which is more than enough for the needs of dealing with approximations of differential equations as difference equations.
Nonlinear SVF difference equation
Let’s now use this nonlinear discrete integrator as a basis for building the discretetime version of the statevariable filter.
Taking the statevariable filter differential equation
and applying the trapezoidal rule to discretize it as a difference equation yields the following system.
The latter difference equation can be then nonlinearized as the following system.
If we want to use this as a practical timestepping integration scheme, we need a way to solve for $bp_{t+1}$ to calculate the filter state from. It is possible to solve it from this system numerically with NewtonRaphson iteration as an implicit integration method.
By substitution $a=\frac{1}{2}h$ and gathering terms the equation for $bp_{t+1}$ becomes
Further
and substituting $C_{t} = bp_{t}  alp_{t}a\sqrt{2}bp_t+au_{t}+au_{t+1}$ as a constant
Using tanh summation identity
and gathering terms to LHS
Further substitutions $D_t=lp_{t}+\frac{1}{2}hbp_{t}$, $E_t=tanh(C_{t})$ and $x=bp_{t+1}$ reduce the equation to
This still looks more complicated than it has to be so substituting the functions
and notating it as a function of x
NewtonRaphson needs a derivative of this function
and of the substituted functions
The iteration for solving $x=bp_{t+1}$ can now be carried over as
from where the full filter state can be integrated as
The rough pseudocode for achieving the full implicit timestepping is as follows.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 

Conclusions
This method is implemented in my statevariable filter module for VCV Rack together with other methods from this blog. The performance is excelent with a constant Q factor across the frequency spectrum and practically unconditional stability, although the SVF filter structure together with nonlinear tanh() integration is stable even with the most basic explicit integration methods.
The mathematical derivation can be carried over to other nonlinear filter structures where stability is not so great (such as the Moog ladder filter) with the same basic steps so hopefully this serves the reader mainly as a good example and a learning stepping stone for further experimentation.
Happy hacks.