Secant method implementation problems

[shoulderman]

Hello!

I have been trying to create an implementation of the secant method in C/C++, but I seem to be coming across a few problems. I don't know if my code is erroneous or I am misunderstanding the algorithm itself.

This is the code:

        const double EPSILON = 0.000001;
	
	double f(double x) { return log(x); }
	
	double FindSecant(double a, double b)
	{
		double fa = f(a), fb = f(b);
		double c = ((b * fa) - (a * fb)) / (fa - fb);
		
		double dx, x1, x2;
		if (fabs(c - a) < fabs(c - b)) { dx = fabs(c - a); x1 = a; x2 = c; }
		else { dx = fabs(c - b); x1 = b; x2 = c; }
		
		if (dx > EPSILON) return FindSecant(x1,  x2);
		return c;
	}

I would greatly appreciate any clarification on this matter.

[E0p9LyTq]

I don't know if this would help, but I found a "Secant Method" in C:

http://www.dailyfreecode.com/code/secant-method-2388.aspx

[ne555]

> I seem to be coming across a few problems.
¿could you be a little more descriptive? ¿what problems?


> I don't know if my code is erroneous or I am misunderstanding the algorithm itself
It would help if you comment your code with your intentions.

[shoulderman]

I am attempting to find the root of ln(x) with a = 0.1 and b = 50. However, c goes through the following values:

c = 18.588534485316014
c = 8.247438135248577
c = 4.351608345118318
c = -4.60894735209671

ln(x) is undefined for x <= 0. It is never even supposed to reach this point.

[Duthomhas]

You are doing something weird with a, b and c on lines 10..12.

You should be recursing on FindSecant( b, c );.

There are a couple things you should be looking out for as well (but you are not):
- Division by zero!
- The method might not converge.

You can check for non-convergence with a simple counter -- too many iterations indicate that you are not converging.

Hope this helps.

[ne555]

> ln(x) is undefined for x <= 0. It is never even supposed to reach this point.
you are throwing a line by two points of the curve.
If you keep in mind that the derivative of your function is `1/x' you may realise that the derivative decreases rapidly (close to 0) and if your two test points are too big you'll intersect the x axis in the negative values.

Your initial approximation should be "close enough" to the root (by instance, use binary search), in your case you don't converge because you fail outside of the domain of the function (yet another check)

[shoulderman]

@Duoas
I was told previously by someone else that I should be recursing with c as well as a or b (whichever is closest to c).
What do you mean by not converging? This has just four iterations until the error appears.

@ne555
I can't keep in mind the derivative since the method isn't supposed to need it. What do you mean by "use binary search"? How would that help?

[TheIdeasMan]

Hi,

Just a couple of things:

on line 7 you make use of the comma operator so I wonder if fa = f(a) won't be evaluated, only the variable fa might be created. I am not sure about this though, but you may be able to humour me for a minute, and put those 2 statements on separate lines. It is good practice to declare and initialise to something 1 variable per line.

Also, won't you need some logic to cover if your function f() returns negative, or at least cover if the sign changes?

A minor thing - you are using fabs(), c++ has std::abs() (in <cmath> )which is overloaded to handle various types.

[Duthomhas]

@shoulderman
You were told wrong. Don't waste time or energy playing with a or b.
Remember, a=x_n-2 and b=x_n-1.
The next iterations's (a,b) should be x_n-1, x_n... which are (b,c).

The secant method has limitations.
https://mat.iitm.ac.in/home/sryedida/public_html/caimna/transcendental/iteration%20methods/secant/secant.html
(Pay attention to the convergence criteria listed as C1 and C2.)

Also, a prerequisite is that f be differentiable, or continuous. The secant method does not require you to know f′, but it must exist.

And x₀ and x₁ need to be well-chosen -- somewhere that f is more-or-less linear and crosses the x-axis.

@TheIdeasMan
Line 7 is fine, though I agree that, stylistically speaking, it is better to put each declaration on its own line. [edit] It's not the comma operator in this case. [/edit]

It doesn't matter if f() returns negative. (Most functions with roots do return negative.)

[ne555]

> I can't keep in mind the derivative since the method isn't supposed to need it.
I was trying to explain how you may fall outside the domain


> What do you mean by "use binary search"? How would that help?
You need a good initial approximation, x0 and x1 should be near the root
https://en.wikipedia.org/wiki/Bisection_method may help you refine them

[shoulderman]

@Duoas
C1 is how many iterations are allowed and C2 is epsilon. For my implementation, only C2 is applicable.
I don't really understand what the actual problem is.

This is some code I found on another website implementing the secant method:

double secant(double a, double b)
{
	do
	{
		double fa = f.f(a), fb = f.f(b);
		if (fa == fb) exit(-1);
		double c = b - ((fb * (b - a)) / (fb - fa));
		a = b;
		b = c;
	} while (fabs(b - a) > EPSILON);

	return b;
}

This does not work for ln(x) between [0.1, 50]. Why is this? What do I do when c becomes negative?

I apologize for my lack of understanding, but I've never really worked with this method before.

[Duthomhas]

Well, no, it won't. Try an (a,b) of (1,2).

(Look at a graph of ln(x) and see what happens when you choose b=50.)

[shoulderman]

@Duoas
It's basically a straight line between ~15 and 50, but there must be some way to make the method work for long intervals, right?

[Duthomhas]

That's beside the point.

Remember, the function does not constrain itself to working only inside the interval [a,b]. That's just the starting point.

If you use 50 as a starting point, then the line crosses the x-axis somewhere to the left of 0, which is outside your function's domain.

As ne555 indicated, the typical way of implementing a root finder with this method is to start with the bisection method to get a pretty good choice for (a, b), and then use the secant method to quickly refine it.

Hope this helps.

[shoulderman]

Thanks for all the help, I think I've found a decent solution now by combining the secant and bisection methods.