# Newton Method of finding roots

I am aware that this is a common question asked, but I want to know some things, which I couldn't find in different Threads about this topic. I am new to mathematica, so the main problem is not the logic behind the solution, but rather the syntax.

So, I wrote the following code, for finding the roots of an arbitrary function, in this case: $$x^2(x-2)$$. I did the following code:

f[x_]=x^2*(x- 2)
p[0]=1.5
Do[p[n+1]=N[p[n]-f[p[n]]/f'[p[n]]],{n,0,30}]
print["n", "p[n]"]
TableForm[Table[{n,p[n]},{n,0,30}]]


I have the following aims and questions:

1. With this method, there is no way to find all the roots at the same time,right?

2. How do I write a function, that I can call it anytime, for any function with real roots. I know how to define a function, but I don't know how to put the arguments, when that argument is a function, plus I don't know how to retrieve the result from it. All I have done is the following:

NewtonRoots[f_,a_,n_]:=

• Why not just use Newton[f_, x0_, n_] := NestList[# - f[#]/f'[#] &, x0, n] then you can do things like: Newton[x |-> x^2 (x - 2), 1.5, 10]. Much simpler. If f doesn't have a simple derivative, you can also look at ResourceFunction["NDerivative"] Commented Jan 4, 2023 at 21:29
• Also explore these functions over at the repository.
– Syed
Commented Jan 4, 2023 at 21:44
• @flinty But this code, if I can say so, and the one Nasser provided, both work well with standard functions, like the one I am trying to solve. The thing is, I thought of functions like x^2 +1 without root, and understandably the end result makes no sense. So I am trying to also implement a condition, that in case the fct has no roots, it lets the user know Commented Jan 4, 2023 at 21:49
• FYI, x^2+1 has complex roots and Newton's method can give the correct answer if you start with an initial guess that is complex. Try Solve[x^2+1==0] and compare it to Newton's method with an initial guess of 1.0+I. Or compare Newton[x|->Sin[x]-2,1.0-I,10] to ArcSin[2.0]. Also "in case the fct has no roots" is undecidable, in general, meaning there is no test. en.wikipedia.org/wiki/Undecidable_problem Commented Jan 5, 2023 at 5:42

Better use := for function definition. Also you need to first take the derivative of f w.r.t. x and after that, substitute the x->p[n] in the result.

ClearAll[f, x, p];
f[x_] := x^2*(x - 2)
p[0] = 1.5
data = Table[p[n + 1] = p[n] - f[p[n]]/f'[x] /. x -> p[n]; {n, p[n]}, {n, 0, 10}]
TableForm[data]


FindRoot[f[x] == 0, {x, 1.5}]


Update

To make it a function do

newtonRoots[f_, x_Symbol, x0_?NumericQ,(nIterations_Integer)?Positive] := Module[{p, data, n},
p[0] = x0;
Table[p[n + 1] = N[p[n] - f[p[n]]/f'[x] /. x -> p[n]]; {n, p[n]}, {n, 0,nIterations}]
];
f[x_] := x^2*(x - 2);
data = newtonRoots[f, x, 1.5, 10];
TableForm[data]


Here is example from help

f[x_] := Sin[x] + Exp[x]
data = newtonRoots[f, x, 0, 10];
TableForm[data]


• But how can i call the function for another function? I am trying something like this: NewtonRoots[f_,a_,n_]:=[Module{}...] where f_ is the argument for a function, a_ is the initial root value and n_ nr. of itterations. Commented Jan 4, 2023 at 20:50
• @imbAF ok, will add that. Commented Jan 4, 2023 at 20:55
• Ok thank you. I don't know the syntax of equalizing the variable that "holds" the function with that of a function within the module. Commented Jan 4, 2023 at 20:56
• Nasser, thank you for your important help. One more thing, when i try to simply find the root of sin(x) and I start with Pi/2 as initial value, after the first interation it goes directly to zero. Does it make sense Commented Jan 4, 2023 at 21:07
• Yes, these ? things in the function definition just adds extra checks on the arguments to make sure they are the correct type. This way if you call the function with input that makes no sense, the function will not be called. It is good to add them, ,your code will be more robust. Commented Jan 4, 2023 at 21:19
Clear[f, x]

f[x_] := x^2*(x - 2)


For multiple roots Map your algorithm onto the starting values.

Row[
Insert[
Grid[
p[0] = #;
Table[
p[n + 1] = p[n] - f[p[n]]/f'[p[n]];
{n, p[n]},
{n, 0, 20}],
Alignment -> ".",
Frame -> All] & /@
{1.5, 1.1}, (* starting values *)
Spacer[10], 2]]


Note that with the first starting value, the algorithm converges more rapidly and there are a lot of unnecessary iterations. It is more efficient to use FixedPointList with the maximum number of iterations limited.

Clear[fixedPt]

fixedPt[p0_, max_ : 50] :=
FixedPointList[# - f[#]/f'[#] &, p0, max]

Row[
Insert[
Grid[
Transpose[
{Range[0, Length[#] - 1], #}] &[fixedPt[#, 20]],
Alignment -> ".",
Frame -> All] & /@
{1.5, 1.1}, (* starting values *)
Spacer[10], 2]]


• But what if the function has no roots, i.e x^2+1. Normally you'd give as info the fact that the function has no roots or that it doesn't converge. Could you use a hidden loop, that iterates again let\s say 1000 times and it generates a single value. Obviously if the function has a root, this value will be make the function have a zero value, but if the function has no roots, even after 1000 iterations, I believe that the end result will be a meaningless value, for which the function won't be close to zero. Would this work, when we consider functions that have no root ? Commented Jan 4, 2023 at 21:46
• To use the Newton method you need a starting value. A good way to get a good starting value is to plot the function. A plot will tell you whether there is a root before you start. Without a plot, if there is no root the sequence for p[n] won't show progress towards convergence. If progress towards convergence is indicated but you have not reached a fixed point, increase the allowed number of iterations. Commented Jan 4, 2023 at 23:08

I understand and respect the fact that you are trying to learn the syntax. But there is more to a programming language than the syntax. Programming languages come with certain biases about how to express ideas. If you really want to learn Mathematica, you should expose yourself to those biases. In Mathematica, do-loops (along with other looping and imperative constructs) are possible, but they aren't considered idiomatic by the Mathematica community at large.

Mathematica has three (at least) "better" ways to approach this problem.

1. Repeatedly apply a function while some condition is satisfied. This concept is called NestWhile. It looks like this [comments are bracketed with (* *)]:

NewtonRoot[initialValue_, function_, tolerance_, iterationMax_] :=
NestWhile[
# - (function[#]/function'[#]) &,(* the core idea in Newton's method *)
initialValue,(* the initial value :) *)
(Abs[f[#]] >= tolerance) &,(* checking whether our value has gotten close to zero *)
1,(* our test function only needs one value at a time *)
iterationMax](* our "safety net" in case the algorithm doesn't terminate normally *)


Test it out:

f[x_] = x^2*(x - 2);
NewtonRoot[1.5, f, .01, 10]


2.00007

If you want to keep track of all intermediate values, then use NestWhileList instead:

NewtonRoot[initialValue_, function_, tolerance_, iterationMax_] :=
NestWhileList[
# - (function[#]/function'[#]) &,
initialValue,
(Abs[f[#]] >= tolerance) &,
1,
iterationMax]

NewtonRoot[1.5, f, .01, 10]


{1.5, 3., 2.4, 2.1, 2.0087, 2.00007}

2. Find a fixed point of a function. This is called FixedPoint. It looks like this:

NewtonRoot[initialValue_, function_, tolerance_, iterationMax_] :=
FixedPoint[# - (function[#]/function'[#]) &, initialValue, iterationMax]


or

NewtonRoot[initialValue_, function_, tolerance_, iterationMax_] :=
FixedPointList[# - (function[#]/function'[#]) &, initialValue, iterationMax]

NewtonRoot[1.5, f, .01, 10] (*using the List version*)


{1.5, 3., 2.4, 2.1, 2.0087, 2.00007, 2., 2., 2.}

3. There is built in functionality for finding roots, FindRoot.

(* probably shouldn't call this NewtonRoot *)
NewtonRoot[initialValue_, function_, tolerance_, iterationMax_] :=
FindRoot[f[\[FormalX]], {\[FormalX], initialValue}, MaxIterations -> iterationMax]

NewtonRoot[1.5, f, .01, 10]


{\[FormalX] -> 2.}

I'm sure there are a bunch of things in these solutions that you don't understand, like the weird syntax with & and # and so forth. I also jumped straight to defining functions with := rather than just demonstrate "raw" solutions. But searching the documentation and digesting these approaches will teach you more about Mathematica than mimicking do loops.

• lericr you are correct, there are plenty of things I don't understand. But what i do is I copy paste a symbol and try to find examples in which this is used. I believe that is a good approach. You believe so? Commented Jan 4, 2023 at 22:41
• Search the documentation for whatever symbol you're interested in. Commented Jan 4, 2023 at 22:45
• Also understand that these weird symbols are syntactic sugar. You can always use the fully expanded forms, and that often makes it easier to understand what's happening. Commented Jan 4, 2023 at 22:46