multiple unknowns in newton raphson loop

[Marvinh]

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);

[mystran]

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.

[Marvinh]

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.

[Urs]

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.

[Marvinh]

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.

[mystran]

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.

[Marvinh]

I had to chat gpt this one looking through andys code as well but getting closer to figuring this out hopefully do a wavefolder next .

Code: Select all

import math
import numpy as np
import matplotlib.pyplot as plt



def tick(inp):
    sr = 44100.0
    gr1 = 1.0/(2200.0);
    gr2 = 1.0/(6800.0);
    m = 1/2;
    minv = 1/m;
    gc1 = 0.47e-6*sr*minv;
    gc2 = 0.01e-6*sr*minv;
    is1 = is2 = 1e-15;
    vt1 = vt2 = 26e-3;
    ic1eq = 0;
    ic2eq = 0;
    arr = []
    for i in range(0,len(inp)):
        v1 = inp[i]
        v2 = 0.0
        v3 = 0.0
        tolerance = 1e-8
        max_iter = 100
        for j in range(0,max_iter):
            vd1 = (0 - v3);
            ed1 = math.exp(vd1/vt1);
            id1 = is1*ed1 - is1;
            gd1 = is1*ed1/vt1;
            id1eq = id1 - gd1*vd1;
            vd2 = (v3 - 0);
            ed2 = math.exp(vd2/vt2);
            id2 = is2*ed2 - is2;
            gd2 = is2*ed2/vt2;
            id2eq = id2 - gd2*vd2;
        # Calculate functions f1 and f2
            f1 = v2 - (gr1 * v1 + gc1 * v3 - ic1eq) / (gr1 + gc1)
            f2 = v3 - (gc1 * v2 + ic1eq + ic2eq + id1eq - id2eq) / (gc1 + gc2 + gr2 + gd1 + gd2)
        
            # Calculate Jacobian matrix J
            J11 = 1
            J12 = -gc1 / (gr1 + gc1)
            J21 = -gc1 / (gc1 + gc2 + gr2 + gd1 + gd2)
            J22 = 1
        
            # Calculate determinant of J
            det_J = J11 * J22 - J12 * J21
        
            if abs(det_J) < 1e-12:
                print("Jacobian is singular, stopping iteration.")
                break
        
            # Inverse of Jacobian matrix
            J_inv = [[J22 / det_J, -J12 / det_J],
                     [-J21 / det_J, J11 / det_J]]
        
            # Newton-Raphson update: [v2_new, v3_new] = [v2, v3] - J_inv * [f1, f2]
            v2_new = v2 - (J_inv[0][0] * f1 + J_inv[0][1] * f2)
            v3_new = v3 - (J_inv[1][0] * f1 + J_inv[1][1] * f2)
        
            # Check for convergence
            if abs(v2_new - v2) < tolerance and abs(v3_new - v3) < tolerance:
                v2, v3 = v2_new, v3_new
                #print(f"Converged after {i+1} iterations.")
                break
        
            # Update for next iteration
            v2, v3 = v2_new, v3_new
        ic1eq += minv*(gc1*(v3 - v2) - ic1eq);
        ic2eq += minv*(gc2*(v3 - 0) - ic2eq);
        arr.append(v3)
    return arr

# Define the parameters
frequency = 100  # frequency of the sine wave in Hz (cycles per second)
sampling_rate = 44100  # samples per second (sampling rate)
duration = 1  # duration of the signal in seconds
period = 1 / frequency
# Generate time array
t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)

# Generate sawtooth wave
wave = 2 * (t % period) / period - 1
print(len(wave))

plt.plot(wave[0:1000])
plt.plot(tick(wave)[0:1000])
plt.show()

Screenshot 2024-11-09 at 12.49.28 PM.png

[Marvinh]

I’m not sure that solving it this way was even necessary as I could have rearranged for just v1 in v2 and v1 into v3 ill post that here as well

[Marvinh]

this is just solving v3 by itself rather than both v2 and v3 at the same time i dont think it was necessary for multiple but might useful

Code: Select all

import math
import numpy as np
import matplotlib.pyplot as plt



def tick(inp):
    sr = 44100.0
    gr1 = 1.0/(2200.0);
    gr2 = 1.0/(6800.0);
    m = 1/2;
    minv = 1/m;
    gc1 = 0.47e-6*sr*minv;
    gc2 = 0.01e-6*sr*minv;
    is1 = is2 = 1e-15;
    vt1 = vt2 = 26e-3;
    ic1eq = 0;
    ic2eq = 0;
    arr = []
    for i in range(0,len(inp)):
        v1 = inp[i]
        v2 = 0.0
        v3 = 0.0
        tolerance = 1e-8
        max_iter = 100
        for j in range(0,max_iter):
            vd1 = (0 - v3);
            ed1 = math.exp(vd1/vt1);
            id1 = is1*ed1 - is1;
            gd1 = is1*ed1/vt1;
            id1eq = id1 - gd1*vd1;
            vd2 = (v3 - 0);
            ed2 = math.exp(vd2/vt2);
            id2 = is2*ed2 - is2;
            gd2 = is2*ed2/vt2;
            id2eq = id2 - gd2*vd2;
        # Calculate functions f1 and f2
            v2=(-((gc2+gd1+gd2+gr2)*(ic1eq-gr1*v1))+gc1*(ic2eq+id1eq-id2eq+gr1*v1))/(gr1*(gc2+gd1+gd2+gr2)+gc1*(gc2+gd1+gd2+gr1+gr2))
            nx=(gc1*(ic2eq+id1eq-id2eq)+gr1*(ic1eq+ic2eq+id1eq-id2eq+gc1*v1))/(gr1*(gc2+gd1+gd2+gr2)+gc1*(gc2+gd1+gd2+gr1+gr2))

            er = abs(nx-v3)
            if er < tolerance:
                v3 = nx
                break
            v3 = nx

        ic1eq += minv*(gc1*(v3 - v2) - ic1eq);
        ic2eq += minv*(gc2*(v3 - 0) - ic2eq);
        arr.append(v3)
    return arr

def tick2(inp):
    sr = 44100.0
    gr1 = 1.0/(2200.0);
    gr2 = 1.0/(6800.0);
    m = 1/2;
    minv = 1/m;
    gc1 = 0.47e-6*sr*minv;
    gc2 = 0.01e-6*sr*minv;
    is1 = is2 = 1e-15;
    vt1 = vt2 = 26e-3;
    ic1eq = 0;
    ic2eq = 0;
    arr = []
    for i in range(0,len(inp)):
        v1 = inp[i]
        v2 = 0.0
        v3 = 0.0
        tolerance = 1e-8
        max_iter = 100
        for j in range(0,max_iter):
            vd1 = (0 - v3);
            ed1 = math.exp(vd1/vt1);
            id1 = is1*ed1 - is1;
            gd1 = is1*ed1/vt1;
            id1eq = id1 - gd1*vd1;
            vd2 = (v3 - 0);
            ed2 = math.exp(vd2/vt2);
            id2 = is2*ed2 - is2;
            gd2 = is2*ed2/vt2;
            id2eq = id2 - gd2*vd2;
        # Calculate functions f1 and f2
            f1 = v2 - (gr1 * v1 + gc1 * v3 - ic1eq) / (gr1 + gc1)
            f2 = v3 - (gc1 * v2 + ic1eq + ic2eq + id1eq - id2eq) / (gc1 + gc2 + gr2 + gd1 + gd2)
        
            # Calculate Jacobian matrix J
            J11 = 1
            J12 = -gc1 / (gr1 + gc1)
            J21 = -gc1 / (gc1 + gc2 + gr2 + gd1 + gd2)
            J22 = 1
        
            # Calculate determinant of J
            det_J = J11 * J22 - J12 * J21
        
            if abs(det_J) < 1e-12:
                print("Jacobian is singular, stopping iteration.")
                break
        
            # Inverse of Jacobian matrix
            J_inv = [[J22 / det_J, -J12 / det_J],
                     [-J21 / det_J, J11 / det_J]]
        
            # Newton-Raphson update: [v2_new, v3_new] = [v2, v3] - J_inv * [f1, f2]
            v2_new = v2 - (J_inv[0][0] * f1 + J_inv[0][1] * f2)
            v3_new = v3 - (J_inv[1][0] * f1 + J_inv[1][1] * f2)
        
            # Check for convergence
            if abs(v2_new - v2) < tolerance and abs(v3_new - v3) < tolerance:
                v2, v3 = v2_new, v3_new
                #print(f"Converged after {i+1} iterations.")
                break
        
            # Update for next iteration
            v2, v3 = v2_new, v3_new
        ic1eq += minv*(gc1*(v3 - v2) - ic1eq);
        ic2eq += minv*(gc2*(v3 - 0) - ic2eq);
        arr.append(v3)
    return arr

# Define the parameters
frequency = 1000  # frequency of the sine wave in Hz (cycles per second)
sampling_rate = 44100  # samples per second (sampling rate)
duration = 1  # duration of the signal in seconds
period = 1 / frequency
# Generate time array
t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)

# Generate sawtooth wave
wave = 2 * (t % period) / period - 1
print(len(wave))

plt.plot(wave[0:100])
plt.plot(tick(wave)[0:100])
plt.plot(tick2(wave)[0:100])
plt.show()

[mystran]

I would highly recommend always adding a bit of parallel conductance (what Spice calls "gMin" on the order of 1e-12) and some series resistance (a few ohms typically) to every pn-junction. This kinda makes the model more realistic too, because real devices have some reverse leakage and series resistance... but that's (kinda) not the main reason.

Rather the main reason is that the parallel conductance essentially forces a minimum derivative and the series resistance caps the maximum. Together they make it a lot less likely that Newton gets stuck or runs into singularity problems. The cost of doing this is basically nothing and even in circuits where it isn't strictly necessary it can still speed up the iteration a bit.

This is one of those little things where making things "just a tiny bit more realistic" happens to also make things numerically more robust, by avoiding some of the numerically problematic special cases that don't really happen in the real-world anyway. Similar things can be said about capacitors and lumped vs. discrete LTPs and .. what not... but Newton in particular doesn't really like things to be too ideal.

[Marvinh]

Thanks for the tips but basically means just add some offset to the resistors even though they maybe the same values or ohms ?

The vt1 vt2 here is the diode thresholds I believe why this is around -.55 floating point to .5 floating point

I change the vt1 vt2 value to get closer to the -1 to 1 threshold I wanted
The capacitors gc1 gc2 are time constants for like a cutoff I assuming according to the sample rate still working through it

Edit : vt1,vt2 doesnt do much as changing the resistor values and i guess you can apply post gain to get it into -1 1 range anyway

[Marvinh]

I took andys diode clipper and a single stage from ken stone for this wave folder. it needed a capacitor but i think it okay for now will probably alias like crazy

edit: code to fix some values like vt (spice default diode is 35.5e-3) and gr2 calculation.

Code: Select all

import math
import numpy as np
import matplotlib.pyplot as plt


def tick(inp):
    arr = []
    for i in range(0, len(inp)):
        is1 = is2 = 1e-15;
        vt1 = vt2 = 35.5e-3;
        gr1 = 1.0/33.2e3;
        gr2 = 1.0/100.0e3;
        v1 = inp[i];
        v2 = 0;
        v3 = 0
        for j in range(0,50):
            vd1 = (0 - v2);
            ed1 = np.exp(vd1/vt1);
            id1 = is1*ed1 - is1;
            gd1 = is1*ed1/vt1;
            id1eq = id1 - gd1*vd1;
            vd2 = (v2 - 0);
            ed2 = np.exp(vd2/vt2);
            id2 = is2*ed2 - is2;
            gd2 = is2*ed2/vt2;
            id2eq = id2 - gd2*vd2
            nx = (id1eq - id2eq + gr1*v1)/(gd1+gd2+gr1)
            er = abs(nx-v2)
            if er < 1.0e-6:
                v2 = nx
                break
            v2 = nx
        arr.append((v2+((v2-v1)*gr1)/(gr2)))
    return arr

# Define the parameters
frequency = 100  # frequency of the sine wave in Hz (cycles per second)
sr = 44100  # samples per second (sampling rate)
duration = 1  # duration of the signal in seconds
period = 1 / frequency
# Generate time array
t = np.linspace(0, duration, int(sr * duration), endpoint=False)

# Generate sawtooth wave
wave = ((np.sin(2 * np.pi * frequency * t)))
print(len(wave))
plt.plot(wave[0:1000])
plt.plot(tick(wave)[0:1000])
plt.show()
        
        

[mystran]

Thanks for the tips but basically means just add some offset to the resistors even though they maybe the same values or ohms ?

What I'm trying to say is that Newton sometimes doesn't like multiple non-linearities next to each other without any resistance between them and when this is the case, adding a little bit of additional series resistance helps. If there is already some anyway, then adding more is not really necessary.

[Music Engineer]

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

Oh - that's great! Thanks a lot!

[Marvinh]

Updated with code for: Virtual Analog Models of the Lockhart and Serge Wavefolders

And then another version i figured out with the paper avoiding LambertW

Code: Select all

struct OddSerge {
    
    
        const double Vt = 25.864e-3;
        const double Is = 2.52e-9;
        const double n = 1.752;
        const double R1 = 33e3;
        const double R2 = 100e3/2;
        const double RL = R2;
        double vp[7];
        double vi[7];
    
        int stage = 6;
        
        OddSerge() {
            for(int i =0; i < 7; i++) {
                vp[i] = 0.0;
                vi[i] = 0.0;
            }
        }

        double process(const double vin) {
            vi[0] = vin*24;
            for(int i = 0; i < stage; i++)
            {
                const double vin = vi[i];
                const double si = boost::math::sign(vi[i]);
                const double v1 = vin - si*n*Vt*boost::math::lambert_w0((RL*Is)/(n*Vt)*exp((si*vin)/(n*Vt)));
                const double v2 = (1+((1/R2)/(1/R1)))*v1 - vin*((1/R2)/(1/R1));
                vp[i] = vi[i];
                vi[i+1] = v2;
            }
            return vi[stage];
        }
   

};

Edit: AADA 1st order version

Code: Select all

struct OddSerge {
        
        double vp[7];
        double vi[7];
        double vpp[7];
    
        const double n = 1.752;
        const double t = 25.864e-3;
        const double s = 2.52e-9;
        const double R = 100e3 / 2.0;
        const double n_t = n * t;
        const double n_t_squared = n_t * n_t;
    
        int stage = 6;
        
        OddSerge() {
            for(int i =0; i < 7; i++) {
                vp[i] = 0.0;
                vi[i] = 0.0;
                vpp[i] = 0.0;
            }
        }
    
    

        double process(const double vin) {
            vi[0] = vin*6;
            for(int i = 0; i < stage; i++)
            {
                if (std::abs(vi[i] - vp[i]) < 1e-7) {
                    double x = vi[i];
                    double sgn = (x > 0) - (x < 0);
                    double lambert = boost::math::lambert_w0((R * s) / n_t * std::exp(sgn * x / n_t));
                    vi[i+1] = 2 * (x - sgn * n_t * lambert) - vi[i];
                    vp[i] = vi[i];
                    
                    
                    // Exponential terms
                    double exp_neg_x = std::exp(-x / n_t);
                    double exp_pos_x = std::exp(x / n_t);
                    
                    // Lambert W terms
                    double W_neg = boost::math::lambert_w0(exp_neg_x * R * s / n_t);
                    double W_pos = boost::math::lambert_w0(exp_pos_x * R * s / n_t);
                    
                    
                    double result1 =  n_t_squared * (
                                                     (sgn - 1) * (0.25 * W_neg * W_neg + 0.5 * W_neg) -
                                                     (sgn + 1) * (0.25 * W_pos * W_pos + 0.5 * W_pos)
                                                     ) + (x * x) / 2;
                    vpp[i] = result1;
                }else{
                    
                    double x = vi[i];
                    double sgn = (x > 0) - (x < 0);  // Equivalent to sign function
                    
                    // Exponential terms
                    double exp_neg_x = std::exp(-x / n_t);
                    double exp_pos_x = std::exp(x / n_t);
                    
                    // Lambert W terms
                    double W_neg = boost::math::lambert_w0(exp_neg_x * R * s / n_t);
                    double W_pos = boost::math::lambert_w0(exp_pos_x * R * s / n_t);
                    
                    
                    double result1 =  n_t_squared * (
                                                     (sgn - 1) * (0.25 * W_neg * W_neg + 0.5 * W_neg) -
                                                     (sgn + 1) * (0.25 * W_pos * W_pos + 0.5 * W_pos)
                                                     ) + (x * x) / 2;
                    
                    vi[i+1] =  2 * (result1 - vpp[i]) / (vi[i] - vp[i]) - vi[i];
                    vp[i] = vi[i];
                    vpp[i] = result1;
                }
            }
            return vi[stage];
        }
};

This one is blazing fast avoids lambertw
can just stack them in series and it doesnt alias as much as even the first order adda of lambert
CPU usage is reduced in half when letting v1 have the previous value

OddSerge2 saturator;

Code: Select all

/*
a = saturator.process(x)
b = saturator.process(a)
c = saturator.process(b)
d = saturator.process(c)
e = saturator.process(d)
f = saturator.process(e)
return f
*/

struct OddSerge2 {
    
    const double n = 1.752;
    const double t = 25.864e-3;
    const double s = 2.52e-9;
    const double R = 1.0/(100e3/2.0);
    const double n_t = n * t;
    const double n_t_1 = 1.0/n_t;
    const double n_t_squared = n_t * n_t;
    
    const double st1 = 1/6.0;
    const double st2 = 1/120.0;
    const double st3 = 1/5040.0;
    const double ct1 = 0.5;
    const double ct2 = 1.0/24.0;
    const double ct3 = 1.0/720.0;
    
    double taylor_sinh(double x){
        double xx = x*x;
        double x5 = xx*xx*x;
        return x + (xx*x)*st1 + (x5)*st2 + (x5*xx)*st3;
    }

    double taylor_cosh(double x){
        double xx = x*x;
        double x4 = xx*xx;
        return 1 + (xx)*ct1 + (x4)*ct2 + (xx*x4)*ct3;
    }
    
    double x,fk,gk,gkeq,v1_new,v1=0.0;
    
    double process(double vi) {
        for(int i = 0; i < 50; i++) {
            x = v1*n_t_1;
            fk = 2*s*taylor_sinh(x);
            gk = 2*s*taylor_cosh(x)*n_t_1;
            gkeq = fk - gk*v1;
            v1_new = (-gkeq + vi*R)/(gk+R);
            if(abs(v1_new-v1) < 1e-6){
                v1 = v1_new;
                break;
            }
            v1 = v1_new;
        }
        return 2*v1-vi;
    }
};

Code: Select all

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import sawtooth

def taylor_sinh(x):
    """5-term Taylor series approximation for sinh(x)."""
    return x + (x**3)/6 + (x**5)/120 + (x**7)/5040

def taylor_cosh(x):
    """5-term Taylor series approximation for cosh(x)."""
    return 1 + (x**2)/2 + (x**4)/24 + (x**6)/720
    
def foldw(vin):
    n = 1.752 
    t = 25.864e-3
    s = 2.52e-9
    R = 1.0/(100e3/2.0)
    R2 = 1.0/(33e3/2.0)
    sr = 44100.0
    C = 1
    arr = np.array([])
    ic1eq = 0
    avg_iter = 0
    for vi in vin:
        x = i
        v1 = 0.0
        j = 0
        for j in range(0,50):
            x = v1/(t*n)
            fk = 2*s*taylor_sinh(x)
            gk = 2*s*taylor_cosh(x)/(t*n)
            gkeq = fk - gk*v1
            v1_new = (-gkeq + (vi+ic1eq)*R)/(gk+R)
            if abs(v1_new-v1) < 1e-6:
                v1 = v1_new
                break
            else:
                v1 = v1_new
        vout =  2*v1-vi
        avg_iter += j
        arr = np.append(arr,vout)
    print("avg_iter %f",avg_iter/len(vin))
    return arr

# Define the parameters
frequency = 12  # frequency of the sine wave in Hz (cycles per second)
sr = 44100  # samples per second (sampling rate)
duration = 1  # duration of the signal in seconds
period = 1 / frequency
# Generate time array
t = np.linspace(0, duration, int(sr * duration), endpoint=False)
amplitude = 2.0
si = amplitude * np.sin(t*frequency*2*np.pi)[0:8820]
saw = amplitude * sawtooth(2 * np.pi * frequency * t)[0:8820]

pp = foldw(si)
pp = foldw(pp)
pp = foldw(pp)
print(max(pp))
plt.plot(pp,"black")



# bb = foldw(saw+1)
# bb = foldw(bb)
# bb = foldw(bb)
# print(max(bb))
# plt.plot(bb,"green")


plt.show()