0

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.

Szabolcs
  • 234,956
  • 30
  • 623
  • 1,263
Jhonny
  • 213
  • 1
  • 2
  • 7
  • Numerical artifacts; change 0.1 to 1/10. – corey979 Sep 08 '16 at 10:19
  • 1
    It's impossible to answer without a complete example. Please edit the question and post one. Refer to http://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? – Szabolcs Sep 08 '16 at 10:23
  • @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]. – Jhonny Sep 08 '16 at 11:04
  • There are built-in commands well suited for this task. Try, e.g., reading in the documentation about Roots and related topics. – corey979 Sep 08 '16 at 11:15
  • @corey979 The problem is that the Roots command doesn't work on the more complicated function I would like to apply a similar loop on later. Right now I would just like to find the fault in the loop described above. – Jhonny Sep 08 '16 at 11:20
  • 1
    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. – Szabolcs Sep 08 '16 at 11:50
  • @Szabolcs And how do I make sure it is a double root? Should I solve f[x] == 0 and f'[x]==0 simultaneously? – Jhonny Sep 08 '16 at 11:56
  • You have two equations but you are solving for a single variable. Usually there will be no solution. If it's a special case where it's know to be a double root in advance, does it matter which one you're solving for, f or f'? Normally you'd have an extra parameter which you are also solving for (makes for a second variable). You might want to review the involved mathematics as well as some basic numerical methods for equation solving before starting this project. In the least take the time to understand the bisection method, Newton's method, secant method. – Szabolcs Sep 08 '16 at 12:23
  • @Szabolcs In case there are more roots then one, i think it matters if I solve for f or f'. If I for instance was only solving for f in the example abow, I could get three roots: -6, 0 and 4. Even if I solve for f', there is more than one point where that condition is met (which is very clear if you plot the graph). All this said I could review the math, but this discussion is going off-topic because all I asked for was an explanation why my loop doesn't do what I want it to do, and not alternative methods. – Jhonny Sep 08 '16 at 12:39
  • Solve just for f[x]==0 to find a set of solutions, then compute f'[x] on these solutions and check which derivatives are equal to zero. – corey979 Sep 08 '16 at 12:41
  • @Alettix It doesn't work for the simple reason that your function is strictly less than zero for any $x>0$ except for $x==4$. And your algorithm never tests f[x] for x==4. Thus the test in the While never fails. You are incrementing in steps of 0.1, you you may have expected that eventually solution will be equal to 4. But 0.1 is an inexact number (floating point). Computers represent floating point numbers in binary. 0.1 is not exactly representable in binary ... – Szabolcs Sep 08 '16 at 12:42
  • 1
    ... 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. – Szabolcs Sep 08 '16 at 12:43
  • Read also this – corey979 Sep 08 '16 at 12:51

1 Answers1

0

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.:

Last@val

{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;
 If[f[solution] == 0 && f'[solution] == 0, Return[solution]]]

Return[4]

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

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

Return[4.]

corey979
  • 23,947
  • 7
  • 58
  • 101