2
$\begingroup$

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_]:=

$\endgroup$
4
  • 1
    $\begingroup$ 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"] $\endgroup$
    – flinty
    Jan 4, 2023 at 21:29
  • 2
    $\begingroup$ Also explore these functions over at the repository. $\endgroup$
    – Syed
    Jan 4, 2023 at 21:44
  • $\begingroup$ @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 $\endgroup$
    – imbAF
    Jan 4, 2023 at 21:49
  • $\begingroup$ 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 $\endgroup$
    – user6552
    Jan 5, 2023 at 5:42

3 Answers 3

7
$\begingroup$

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]

Mathematica graphics

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

Mathematica graphics

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

enter image description here

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

Mathematica graphics

$\endgroup$
9
  • $\begingroup$ 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. $\endgroup$
    – imbAF
    Jan 4, 2023 at 20:50
  • $\begingroup$ @imbAF ok, will add that. $\endgroup$
    – Nasser
    Jan 4, 2023 at 20:55
  • $\begingroup$ 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. $\endgroup$
    – imbAF
    Jan 4, 2023 at 20:56
  • $\begingroup$ 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 $\endgroup$
    – imbAF
    Jan 4, 2023 at 21:07
  • 1
    $\begingroup$ 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. $\endgroup$
    – Nasser
    Jan 4, 2023 at 21:19
3
$\begingroup$
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]]

enter image description here

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]]

enter image description here

$\endgroup$
2
  • $\begingroup$ 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 ? $\endgroup$
    – imbAF
    Jan 4, 2023 at 21:46
  • $\begingroup$ 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. $\endgroup$
    – Bob Hanlon
    Jan 4, 2023 at 23:08
2
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$
    – imbAF
    Jan 4, 2023 at 22:41
  • $\begingroup$ Search the documentation for whatever symbol you're interested in. $\endgroup$
    – lericr
    Jan 4, 2023 at 22:45
  • $\begingroup$ 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. $\endgroup$
    – lericr
    Jan 4, 2023 at 22:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.