Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I've only had few previous experiences on Mathematica before. So basically an entry-level question that harasses a rookie like me.

I made the following attempt using bisection method to calculate the real root for the function $f(x)=-x^3+2x^2-2$ with basic while loop and if:

a = -5; b = 5; m = (a + b)/2; While [ Abs[f[m]] > 0.000001, If[f[m]*f[a] < 0, b = m, a = m]; m = (a + b)/2]

The answer turns out to be $-(14080895/16777216)$ that is correct.

Then, I wanted to create a animation that shows the change of position of a dot along the function curve for each iteration of the bisection method. But I tried several times with different codes that all fails to do so. One confusion is I am not sure how to set up a stepsize here which involves not a constant number but "While" and "If".

Would you please help me animate the iteration with while and if?

share|improve this question
    
Two suggestions to make your bisection routine more robust: 1.) compute m as m = a + (b - a)/2, and 2.) change the sign change criterion to Sign[f[m]] != Sign[f[a]]. –  J. M. May 28 '12 at 0:56
add comment

3 Answers

First, what follows is an improvement of your bisection routine, which avoids needless overflow and repeated function evaluations. As with Sjoerd, we use the Sow[]/Reap[] pair for collecting intermediate results:

f[x_] := -2 + 2 x^2 - x^3;
acc = 10.^-6;
a = -5; b = 5;
u = f[a]; v = f[b]; m = $Failed;
pts = Reap[While[NonNegative[u] =!= NonNegative[v],
         m = a + (b - a)/2; h = f[m]; Sow[{m, h}];
         If[NonNegative[u] =!= NonNegative[h],
            (* then *) b = m; v = h,
            (* else *) a = m; u = h];
         If[Abs[h] <= acc, Break[]]; 
         ]][[2, 1]];

I used OP's absolute termination criterion in the code above, but a relative criterion like Abs[b - a] <= (Abs[a] + Abs[b])*tol is usually more preferable in practice.

Having gathered up all required intermediate results in pts, we can elect to either show progress all at once:

Plot[f[x], {x, -5/2, 0}, 
 Epilog -> {Arrow /@ Partition[pts, 2, 1], {AbsolutePointSize[7], Red,
     Point[First[pts]]}, {AbsolutePointSize[7], Blue, 
    Point[Last[pts]]}}, Frame -> True]

bisection plot

or show it gradually:

ListAnimate[
 MapIndexed[
  Plot[f[x], {x, -5/2, 0}, 
    Epilog -> {Arrow[#1], {AbsolutePointSize[7], Red, 
       Point[First[pts]]}, {AbsolutePointSize[7], Blue, 
       Point[Last[pts]]}}, Frame -> True, 
    PlotLabel -> StringForm["Iteration `1`", First[#2]]] &, 
  Partition[pts, 2, 1]]]

bisection animation

In both cases, we use Arrow[] objects to indicate the direction of progress, a red point to denote the midpoint of the initial interval, and a blue point to denote the final result.

share|improve this answer
add comment

Adding a Sow to your loop let's you store the intermediate values of m. You can grab them later with Reap and animate them with Animate:

a = -5; b = 5; m = (a + b)/2; 
steps = 
 Reap[While[Abs[f[m]] > 0.000001, 
    Sow[m] If[f[m]*f[a] < 0, b = m, a = m]; m = (a + b)/2]][[2, 1]];

Animate[
 Plot[f[x], {x, -2.5, -0.5}, 
  Epilog -> {Red, PointSize[0.03], 
    Point@{steps[[i]], f@steps[[i]]}}], {i, 1, Length[steps], 1}
 ]

enter image description here

share|improve this answer
add comment

Perhaps a more Mathematicaish way of doing more or less the same is with NestWhile A possible code could be

NestWhile[
 With[{m = Mean[#]},
   If[
    Negative[f[First@#] f[m]], {First@#, m}, {m, Last@#}
    ]] &, {-5, 5},
 Abs@f@Mean@# > 0.000001 &
 ]

If you have your code written with NestWhile, what you want is as simple as changing NestWhile to NestWhileList.

steps = NestWhileList[
   With[{m = Mean[#]},
     If[
      Negative[f[First@#] f[m]], {First@#, m}, {m, Last@#}
      ]] &, {-5, 5},
   Abs@f@Mean@# > 0.000001 &
   ];

ListLinePlot[Mean /@ steps, PlotRange -> Full]

Mathematica graphics

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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