So I would like to write a code which iteratively finds the double root (a root which is also a maximum or a minumum) of functions. I would like to apply this to more complex functions Mathematica struggle with handling, but I thought I would start with a simple one:

f[x_] = (-1)*x*(x - 4)^2*(x + 6)

Now this function obviously has a double root in x = 4. I wanted to write a loop which would test different values of x until it finds a specific x for which f[x]==0 and f'[x]==0 (the conditions of a double root). I wrote the following:

solution = 0.1
While[f[solution] < 0, solution += 0.1; 
 If[f[solution] == 0 && f'[solution]==0, Return[solution]]

As I see it, this loop should start testing x = 0.1, if f[0.1] < 0 (which I have made sure it is in the whole interval 0.1<x<4) it should increase the tested value with 0.1 and try again until f[x] no longer is negative. Also, when we reach the point where the conditions for the double root are met, the loop should return solution. However, it doesn't seem to work at all. I have let the loop run for half an hour without any result. Could somebody please tell me what's wrong? The only problem I can think about is if for some specific solution we have f[solution] < 0, f[solution+0.1] > 0 meaning that the condition f[solution] == 0 will never be met. This could of course be the case with more complex function, but with the current f[x] this shouldn't be a problem as 4 = 0.1 + 39*0.1.

EDIT: I do know that there are several commands which can find the root to this function. However, all I would like to know is simply why this loop does not work.

  • $\begingroup$ Numerical artifacts; change 0.1 to 1/10. $\endgroup$
    – corey979
    Sep 8 '16 at 10:19
  • 1
    $\begingroup$ It's impossible to answer without a complete example. Please edit the question and post one. Refer to sscce.org for guidance. That said, such a method will never work for finding roots. Say you want to solve x^2 - Sqrt[2] == 0 and you simply test all x values in the interval $[0,2]$ with a certain step size dx. Do you expect to ever find a value for which the equality truly holds? $\endgroup$
    – Szabolcs
    Sep 8 '16 at 10:23
  • $\begingroup$ @szabolcs Firstly, the code included is really the complete code I have written. So It is quite unclear what you are asking for with "complete example". Secondly, as you can see in the question formulation I am aware of that the step size is a problem when applied to other functions, but I believe it can be approximately fixed by changing f[x] == 0 to f[x] >= 0 and something similar with f'[x]. $\endgroup$
    – Jhonny
    Sep 8 '16 at 11:04
  • 1
    $\begingroup$ FindRoots works on any functions. If this is really the complete code, it is good to note that Return is not appropriate here. It should be used in a function and it's almost never needed in Mathematica. This code is not part of a function. $\endgroup$
    – Szabolcs
    Sep 8 '16 at 11:50
  • 1
    $\begingroup$ ... the same way as $1/3$ is not exactly representable in decimal. How many $3$ do you need after the decimal point in $0.33333\dots$? Adding up 0.1 many times will never give an exact 4 because 0.1 is not exactly 1/10. That's the explanation of why this loop doesn't stop. The fundamental issue however is that such an algorithm is not capable of finding roots. The failure you see here is just one manifestation of that. $\endgroup$
    – Szabolcs
    Sep 8 '16 at 12:43

I modified slightly your code to illustrate what its output is (my codes in this post are poor, but they are meant to be a fast and easy demonstration of what does MMA in fact do; they are not meant to be recomendable codes):

solution = 0.1; val = {};
While[solution < 4., solution += 0.1;
 If[f[solution] == 0 && f'[solution] == 0, solution, 
  AppendTo[val, {solution, f[solution], f'[solution]}]]]

When you display what happens when solution=4.:


{4., -1.26218*10^-28, -1.42109*10^-13}

you'll see that your conditions ==0 are never met. Like in my comment, when you change 0.1 into 1/10 then the code works because it finds an exact match:

solution = 1/10;
While[f[solution] < 0, solution += 1/10;


You'd need to change approximate zeros into exact zeros using Chop to have them matched with 0:

solution = 0.1;
 If[Chop[f[solution]] == 0 && Chop[f'[solution]] == 0, 


