Re: A question on Newton's Method
- From: Jon Harrop <usenet@xxxxxxxxxxxxxx>
- Date: Sun, 03 Apr 2005 16:15:24 +0100
Roman Werpachowski wrote:
> On the Sat, 2 Apr 2005 20:24:58 -0700, James Van Buskirk wrote:
>> I can't believe you're wasting your time debating someone who
>> tinks that a lisp-like language is suitable for a beginner in
>> numerical analysis.
What is Lisp-like about Mathematica and OCaml?
>> How many times have you seen Mathematica
>> just go into the tank and never come back? Not to mention
>> that the syntax is just sooo obscure, for e.g. D[f[x],x] doesn't
>> yield the derivative of function f.
Yes it does:
In[1]:= f[x_] = x^3 - x - 1
Out[1]= -1 - x + x^3
In[2]:= D[f[x], x]
Out[2]= -1 + 3 x^2
It even represents function derivatives in conventional mathematical
notation:
In[3]:= D[g[x],x]
Out[3]= g'[x]
In[4]:= f''[x]
Out[4]= 6 x
>> It takes way more time to
>> figure out what the syntax is to ask Mathematica for the answer
>> than it does to just work out the answer by hand. It's a total
>> waste of time and math departments should just forget about
>> Mathematica and TI-85s that they seem to be so hung up on and
>> go back to teaching math or whatever it was they used to do.
I wouldn't hesitate to recommend Mathematica, both to beginners and to
people wanting to do serious work.
> Yes, but what's the point of discussing Mathematica in the context of
> NUMERICAL analysis? It's not its main strength, is it?
I assume you're thinking of Mathematica as a symbolic maths package. It
certainly does symbolic maths, and is very good at it, but it also does
lots of other useful stuff. It contains a huge repository of information on
maths as well, including numerical analysis.
Mathematica contains native support for a wide variety of useful numerical
types (arbitrary-precision integers, floating-point real and complex,
rationals etc.) as well as a huge number of numerical algorithms (FindRoot,
FindMinimum, NSolve etc.). These make it a great tool for someone learning
numerical analysis.
OCaml isn't as well suited to numerical analysis as Mathematica because it
wasn't designed for this. However, OCaml is a much more appropriate
language for a beginner to learn than Fortran because it lets you
concentrate on the maths and forget about arbitrary restrictions imposed by
the language.
For example, in Mathematica you can find the root of f(x) symbolically:
In[5]:= First[Solve[f[x] == 0, x]]
Out[5]= {x -> (27/2 - (3*Sqrt[69])/2)^(1/3)/3 + ((9 +
Sqrt[69])/2)^(1/3)/3^(2/3)}
but you can also find the root numerically:
In[5]:= FindRoot[f[x]]
Out[5]= {x -> 1.32472}
This is the easiest way to find a root so it is great for researchers. But a
beginners wants to learn how to implement FindRoot themselves. So let's
define our own root finding algorithm:
In[6]:= NewtonRaphson = Function[{f}, Function[{x}, x - f[x]/f'[x]]]
applying it once to an initial estimate of the root of f(x) gets closer to
the actual root:
In[7]:= NewtonRaphson[f][2.]
Out[7]= 1.54545
Applying our root finder until the result stops changing (i.e. to fixed
point) gets the correct position of the root:
In[8]:= FixedPoint[NewtonRaphson[f], 2.]
Out[8]= 1.32472
You can do the same thing in OCaml, it is just a little more work. You'd
want to start by defining your own numerical derivative function:
# let d f x =
let delta = sqrt epsilon_float in
(f (x +. delta) -. f (x -. delta)) /. (2. *. delta);;
val d : (float -> float) -> float -> float = <fun>
Then you'd define your own version of FixedPoint:
# let rec fixed_point f x1 =
let x2 = f x1 in if x2 = x1 then x2 else fixed_point f x2;;
val fixed_point : ('a -> 'a) -> 'a -> 'a = <fun>
Now you can define Newton-Raphson:
# let newton_raphson f x = x -. f x /. d f x;;
val newton_raphson : (float -> float) -> float -> float = <fun>
So that first iteration gives you:
# let f x = x ** 3. -. x -. 1.;;
val f : float -> float = <fun>
# newton_raphson f 2.;;
- : float = 1.54545454545454541
and iterating until fixed point gives you:
# fixed_point (newton_raphson f) 2.;;
- : float = 1.32471795724474606
IMHO, this illustrates the relevant numerical analysis very well (in
Mathematica you can even visualise the algorithm). We've actually used many
features of these modern languages:
1. The value returned by the subexpressions "NewtonRaphson[f]" and
"newton_raphson f" is actually a function.
2. We didn't have to explicitly deallocate anything because everything is
deallocated for us.
3. We've used an interactive environment to get feedback on the definitions
we're making. OCaml even checks the definitions and tells us the types it's
inferred.
4. The "fixed_point" function is polymorphic => it works for any type 'a.
5. The actual output of Mathematica is typeset mathematics, I've just
converted it to ASCII to post it here.
A Fortran implementation of this kind of stuff would be much longer but
wouldn't actually convey any more useful information, indeed it would
distract the student from the point of the exercise. As the difficulty
level increases, the amount of code required to solve a given problem in
old fashioned languages explodes. Hence, languages like OCaml are much
better for more difficult problems as well.
--
Dr Jon D Harrop, Flying Frog Consultancy
http://www.ffconsultancy.com
.
- Follow-Ups:
- Re: A question on Newton's Method
- From: James Van Buskirk
- Re: A question on Newton's Method
- References:
- A question on Newton's Method
- From: David M
- Re: A question on Newton's Method
- From: Jon Harrop
- Re: A question on Newton's Method
- From: beliavsky
- Re: A question on Newton's Method
- From: James Van Buskirk
- Re: A question on Newton's Method
- From: Roman Werpachowski
- A question on Newton's Method
- Prev by Date: Mathematica for numerical analysis (was Re: A question on Newton's Method)
- Next by Date: Re: can somebody verify this C program which calls dsaupd_ ? (longish)
- Previous by thread: Mathematica for numerical analysis (was Re: A question on Newton's Method)
- Next by thread: Re: A question on Newton's Method
- Index(es):
Relevant Pages
|