Re: How to formulate this while- loop in mathematica ?

*To*: mathgroup at christensen.cybernetics.net*Subject*: [mg1516] Re: How to formulate this while- loop in mathematica ?*From*: villegas (Robert Villegas)*Date*: Wed, 21 Jun 1995 21:04:08 -0400*Organization*: Wolfram Research, Inc.

> i would like to implement a newton-like approximation in mathematica, > > which involves the following program fragment in C-like notation: > > while( (condition1(a,b) > eps) && condition2(a,b) > eps ) { > > /* some numerical computations */ > > a = ..... ; > b = ..... ; > } > > eps gives the tolerance when the loop can terminate. > > My problem is to how to formulate this algorithmic skeleton in mathematica ? ================ Method 1: While ================ You could map this directly to Mathematica using While: In[1]:= ? While While[test, body] evaluates test, then body, repetitively, until test first fails to give True. Something like {a, b} = {a0, b0}; While[condition1[a, b] > eps && condition2[a, b] > eps, (* some numerical computations *) a = aNew; b = bNew; ] {a, b} ===================== Method 2: FixedPoint ===================== A different implementation can be gotten from FixedPoint. First, formulate the body of the 'while' loop as a function that takes an ordered pair {a, b}, does the computations, and returns an ordered pair containing the new values of a and b: successor[{a_, b_}] := {aNew[...], bNew[...]} Do not assign new values to a and b, just formulate the new pair in terms of the previous pair {a, b} and let it be the return value of 'successor'. Really all you're doing is extracting the body of the 'while' loop into an ordered pair formula. Then formulate the conditions as functions that take an ordered pair {a, b} and return True or False: condition1[{a_, b_}] := <some testing> condition2[{a_, b_}] := <more testing> With the functions defined, the procedure is a one-liner: FixedPoint[successor, {a0, b0}, SameTest->(condition1[#2] <= eps || condition2[#2] <= eps &) ] That's all! As for how to use SameTest, and why I chose the function I did, read on if you are interested in details. The spirit of FixedPoint is to iterate a function until the result no longer changes. More precisely, what it does is iterate until two consecutive terms satisfy some condition, which is usually a function to decide if the two are close enough to be considered equivalent. go until test[ x(n-1), x(n) ] returns True SameTest specifies this comparison function (if you don't want to use the default one, which is essentially SameQ). The comparison function will be given two arguments: the current term and the previous term, or x(n-1) and x(n) as I wrote them above. Very often the most convenient way to express the comparison is as a pure function, using #1 and #2 for the variables (see p. 207 of the Mathematica book). For instance, FixedPoint[f, x, SameTest -> (Abs[#1 - #2] < 10^-4 &) ] would go until two consecutive terms were less than 10^-4 distant. But the function doesn't have to perform a comparison of the two arguments; i.e. it isn't required to use both arguments that are given to it. It could simply test the second argument (= most recent term) for some property. That is what I did using just #2 above: #2 will represent {a, b} at the current step, and will be run through both conditions. FixedPoint used this way becomes an Until loop, in effect. =============================== Newton's method with FixedPoint =============================== Here's a quick and dirty implementation of a minimum-finder using FixedPoint (this doesn't use any special SameTest) (* This has a local minimum at x = Pi, y = E. *) In[1]:= g[x_, y_] =(x - Pi)^2 + (y - E)^2 + 1 2 2 Out[1]= 1 + (-Pi + x) + (-E + y) (* We want to know where the gradient vanishes: *) In[2]:= grad[x_, y_] = {D[g[x, y], x], D[g[x, y], y]} Out[2]= {2 (-Pi + x), 2 (-E + y)} In[3]:= jacob[x_, y_] = Together @ Outer[D, grad[x, y], {x, y}] Out[3]= {{2, 0}, {0, 2}} (* This finds the next approximation of the root of the gradient: *) In[4]:= next[{x_, y_}] := {x, y} + LinearSolve[jacob[x, y], -grad[x, y]] (* Iterate this approximation until it stops moving: *) In[5]:= FixedPoint[next, {1., 1.3}] Out[5]= {3.14159, 2.71828} (* Compare to the result of FindMinimum: *) In[6]:= FindMinimum[g[x, y], {x, 1.}, {y, 1.3}] Out[6]= {1., {x -> 3.14159, y -> 2.71828}} Robby Villegas