multiple unknowns in newton raphson loop

DSP, Plugin and Host development discussion.
RELATED
PRODUCTS

Post

How do you converge newton raphson when there is more than one value to solve for? specifically in this example. what part of the equation do I focus on? jacobian matrix necessary?

From andy simpers adc talk.

Code: Select all

tick(input):
v1 = input;
v2 = 0;
v3 = 0;
nr loop:
  vd1 = (0 - v3);
  ed1 = exp(vd1/vt1);
  id1 = is1*ed1 - is1;
  gd1 = is1*ed1/vt1;
  id1eq = id1 - gd1*vd1;
  vd2 = (v3 - 0);
  ed2 = exp(vd2/vt2);
  id2 = is2*ed2 - is2;
  gd2 = is2*ed2/vt2;
  id2eq = id2 - gd2*vd2;
  v2 = solve for v2 using equations; //able to get this
  v3 = solve for v3 using equations; // able to get this
ic1eq += minv*(gc1*(v3 - v2) - ic1eq);
ic2eq += minv*(gc2*(v3 - 0) - ic2eq);

Post

With ordinary scalar-valued NR you basically fit a tangent-line at the current estimated operating point, then use that to solve for a new operating point and fit a new tangent-line, until it converges.

In the case of multiple non-linearities, we can do the same by treating the whole thing as a matrix equation and fitting a tangent (hyper-)plane during each iteration. In practice, when the non-linearities are separate scalar non-linearities, each of them can be separately fit as a separate tangent-line along the relevant dimension, by the magic of linear algebra that's the same as fitting a (hyper-)plane, but it's the partial derivatives (that determine the slopes of the lines/planes) that is the Jacobian.

Here is the method written as typical in math
https://calculus.subwiki.org/wiki/Newto ... r_variable

Yet, observe that even though the formula is written in terms of inverse of the Jacobian, in practice you do not need to actually invert the matrix, simply solving the usual Ax=b matrix-vector equation is sufficient. It is also possible to partition the system matrix into linear and non-linear parts, such that you only need so solve Ax=b for the non-linear part each iteration (and then solve the rest just once), but I would recommend working out the details of the simple version (solve the whole thing directly) first, because understanding how that works is sort of essential for understanding how to make it more efficient (also by not explicitly storing matrixes, but rather just figuring out what arithmetics is necessary).

Unfortunately the path to enlightenment here is to struggle through the derivation from the circuit equations to the resulting code, because the code itself doesn't make any sense, until you understand exactly how it is derived and the only way to really understand it properly is to do the derivation yourself.

Post

Thank you i will try to see what I can do. Looks like I have to solve for everything at the same time not just v2 v3.

Post

Apparently I will never finish my blog post on the topic, but I did once post the code for 4-way vectorial NR:

viewtopic.php?t=459955

Note that this one does not invert the Matrix. Instead of dividing by it, it's moved as product "as is" to the left side of the equation. Then, the equation system is solved using linear algebra.

There is plenty of space for optimisation, one can make the code fly.

Post

Yes forgot to mention that I am looking at it now. helpful resource. the soft clipper above sounds great for a rough triangle to sin just want to hear what it will sound like filtered properly.

Post

Urs wrote: Thu Sep 19, 2024 7:28 am Note that this one does not invert the Matrix. Instead of dividing by it, it's moved as product "as is" to the left side of the equation. Then, the equation system is solved using linear algebra.
It is convenient to mathematically notate a multiplication by inverse, but it's almost always better (faster, more numerically stable) to simply solve the equation directly instead.

So whenever you read multiplication by inverse, unless you're precomputing for hardware accelerated matrix multiply (eg. inverse transforms for a GPU), you can basically read it as "solve this equation" and even if you're solving it multiple times, it's probably better to compute Cholesky/LU/etc instead of precomputing an inverse.

Return to “DSP and Plugin Development”