About Newton Raphson

[Audiority]

Hi guys,

I'm using a very simple NR implementation to find the root of a nonlinear function. Something like

Code: Select all

double A = x * x * x; //x is the input
double dA = 3 * x * x;

double delta = A / dA;
output = output - delta;

This block of code iterates for a predefined number of times, but if the input is too high the NR diverges. I read about the damped NR, but I don't know if it can help.

Any advice?

Cheers,
Luca

[stratum]

I couldn't find anything wrong, what was your exact code?

Code: Select all

#include "math.h"
#include "iostream"

int main()
{
	double x = 1E+10;
	double delta = 0;
	double A = 0.0;
	double dA = 0.0;
	do
	{
		std::cerr << x << std::endl;
		A = x * x * x; //x is the input
		dA = 3 * x * x;

		delta = A / dA;
		x = x - delta;
	} while (fabs(delta) > 1E-5);
	std::cerr << "FINAL ROOT:" << x << std::endl;
	std::cerr << "FINAL RESULT:"<< A << std::endl;
	return 0;
}

[Audiority]

In my case I have A which is a sinh operation, and dA with a cosh. The problem is when the input is over 10, the results of A and dA are very big and the algorithm diverges. I got a better result by using a single NR iteration as my initial estimation, but it still diverge for high gain signals.

[stratum]

sinh(1E+100) returned not-a-number when I've tried.
http://www.efunda.com/math/hyperbolic/d ... ?name=sinh

edit: The one that was not-a-number is delta, not sinh.

Code: Select all

#include "math.h"
#include "assert.h"
#include "iostream"

int main()
{
	double x = 1E+100;
	double delta = 0;
	double A = 0.0;
	double dA = 0.0;
	do
	{
		std::cerr << x << std::endl;
		A = sinh(x); //x is the input
		assert(!isnan(A));

		dA = cosh(x);
		assert(!isnan(dA));

		delta = A / dA;
		assert(!isnan(delta));

		x = x - delta;
	} while (fabs(delta) > 1E-5);
	std::cerr << "FINAL ROOT:" << x << std::endl;
	std::cerr << "FINAL RESULT:" << A << std::endl;
	return 0;
}

[mystran]

This block of code iterates for a predefined number of times, but if the input is too high the NR diverges. I read about the damped NR, but I don't know if it can help.

See proof of quadratic convergence for NR.

Note that your function x^3 in OP has root of multiplicity 3 at zero (so first two derivatives are zero at the root), so you won't be getting quadratic convergence anyway (ie. there's not much point using NR for this particular example), but really the reason it diverges (also with sinh) is probably that with "input too high" the initial guess is not "sufficiently close" (which in this case is a technical term).

[stratum]

If you keep values sane and add a few more assertions to be sure, it converges, maybe by just luck.

Code: Select all

#include "math.h"
#include "assert.h"
#include "iostream"

int main()
{
	double x = 700;
	double delta = 0;
	double A = 0.0;
	double dA = 0.0;
	do
	{
		std::cerr << x << std::endl;
		A = sinh(x); //x is the input
		assert(!isnan(A));

		dA = cosh(x);
		assert(!isnan(dA));
		assert(dA != 0.0);

		delta = A / dA;
		assert(!isnan(delta));

		x = x - delta;
	} while (fabs(delta) > 1E-5);
	std::cerr << "FINAL ROOT:" << x << std::endl;
	std::cerr << "FINAL RESULT:" << A << std::endl;
	return 0;
}

[Audiority]

replacing my for loop with a while loop using an epsilon helps a bit, but in the end the delta becomes NaN because both A and dA are -Inf.

[stratum]

It always converges for any initial value of x between -710 and +710, and considering the shape of sinh, that looks normal to me.

[mystran]

May I ask a stupid question: if you want to solve for x=sinh(y) then why use NR at all?

The analytic solution is y=asinh(x)=log(x+sqrt(x*x+1)) which is cheaper to calculate than a single iteration of Newton.

Or am I completely missing something here?

[stratum]

May I ask a stupid question: if you want to solve for x=sinh(y) then why use NR at all?

Just learning probably. I have started experimenting with the code Audiority has posted for that reason.

[Audiority]

Basically, I'm trying to avoid NaN and Inf while working with high gain signals on diode based circuits. I'm just wondering what is the best approach if a NR solver diverges. On Numerical Recipes it's suggested to use a combination of NR and bisection, but I don't have a predefined range that is required by the book.

[stratum]

Basically, I'm trying to avoid NaN and Inf while working with high gain signals on diode based circuits. I'm just wondering what is the best approach if a NR solver diverges. On Numerical Recipes it's suggested to use a combination of NR and bisection, but I don't have a predefined range that is required by the book.

Solve for the other variable http://www.simulanalog.org/statevariable.pdf

edit: Actually seems like this wasn't your problem. The virtual diode does not have to behave like exactly like a real one though, so the fuction is not necessarily sinh(x). One will never see 700 volts across a pair of real silicon diodes either.

[mystran]

Basically, I'm trying to avoid NaN and Inf while working with high gain signals on diode based circuits. I'm just wondering what is the best approach if a NR solver diverges. On Numerical Recipes it's suggested to use a combination of NR and bisection, but I don't have a predefined range that is required by the book.

The diode current (of the ideal diodes) blows up once the voltage over the diodes gets too high. Real-life diodes would burn in the equivalent situation. In practical circuits, there is always some series resistance of one type or another than eventually becomes the limiting factor and this series resistance is really the only reason you want to use NR to solve these things (because once you add it in, the thing is no longer invertible).

NR can still throw the iterations off the "safe range" though, in which case you might want to use the inverse (of the diode or similar component) to get a more sensible estimate for the next iteration. For this to work though, you need estimates for both sides of the equation though (which is normally not a huge deal in circuit sim, as you can just use the previous Newton linearisation to solve back the current that was implied by the non-sense voltage).

I honestly have no idea if this works for a lumped sinh() component as well as it does for single diodes (even if you put two in parallel to get a bipolar clamp). What I do know is that a pair of transistors (even with just the forward diodes for each) tends to converge more reliably than anything I've been able to do with tanh() directly. YMMV.

[Audiority]

Thanks guys. Let's take a simple diode clipping circuit with a resistor and two diodes in parallel. It goes fine up to 30dB, but once over it blows up. Having gains over 30dB is not so unusual when you want to replicate a clipping stage that features an op-amp. That's clearly something wrong I'm doing, but to figure it out maybe I have to get back to spice first or make some real measurements.

[stratum]

It's not very likely that you'll find an opamp that operates at a voltage that can blow up the solver. No opamp circuit can output 700volts. If you see that voltage at the output then something else is wrong with the model.