Mike Williams wrote:
> 
>...
> Even if someone tells me what p(x) is, it's still going to be horrible
> to solve for x. I'd be inclined to use some sort of numerical method to
> obtain an approximation.

I too thought there was a need to solve this
numerically, so I started working on some
macros that'll do the work.

See code below. The macros are not finished.
(E.g.: No detection of non-convergence yet.)

I have a suspicion that the RegulaFalsi 
method is the best candidate for the job.


Tor Olav


// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 =
// Copyright 2001 by Tor Olav Kristensen
// mailto:tor### [at] hotmailcom
// http://www.crosswinds.net/~tok
// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 =

#version 3.5;

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 =
// The macros

// Fn: The function to find a root for.
// Xa: x-value that are believed to be close to the root.
// Epsilon: Value that are "close enough" to zero for your use.
#macro NewtonRaphson(Fn, Xa, Epsilon)

  #local DerFn = function(x, H) { (Fn(x + H) - Fn(x - H))/2/H } 
  #local IterFn = function(x, H) { x - Fn(x)/DerFn(x, H) }
  #local hh = Xa/100;
  #local x0 = Xa;
  #local f0 = Fn(x0);
  #while (abs(f0) > Epsilon)
    #local x0 = IterFn(x0, hh);
    #local f0 = Fn(x0);
    #debug ":"
  #end // while
  
  (x0)
  
#end // macro NewtonRaphson


// Fn: The function to find a root for.
// Xa: x-value that are believed to be on one side of the root.
// Xb: x-value that are believed to be on the other side of the root.
// Epsilon: Value that are "close enough" to zero for your use.
#macro RegulaFalsi(Fn, Xa, Xb, Epsilon)

  #local x1 = Xa;
  #local f1 = Fn(x1);
  #local x2 = Xb;
  #local f2 = Fn(x2);
  #local x3 = (x1*f2 - x2*f1)/(f2 - f1);
  #local f3 = Fn(x3);
  #if (f1*f2 > 0)
    #error "Macro RegulaFalsi: Endpoints does not bracket the root"
  #else
    #while (abs(f3) > Epsilon)
      #if (f1*f3 < 0)
        #local x2 = x3;
        #local f2 = f3;
      #else
        #if (f2*f3 < 0)
          #local x1 = x3;
          #local f1 = f3;
        #else
          #error "Macro RegulaFalsi: Maybe Epsilion is negative ?"
        #end // if
      #end // if
      #local x3 = (x1*f2 - x2*f1)/(f2 - f1);
      #local f3 = Fn(x3);
      #debug ":"
    #end // while
  #end // if
  
  (x3)

#end // macro RegulaFalsi


// Fn: The function to find a root for.
// Xa: x-value that are believed to be on one side of the root.
// Xb: x-value that are believed to be on the other side of the root.
// Epsilon: Value that are "close enough" to zero for your use.
#macro Bisection(Fn, Xa, Xb, Epsilon)

  #local x1 = Xa;
  #local f1 = Fn(x1);
  #local x2 = Xb;
  #local f2 = Fn(x2);
  #local x3 = (x1 + x2)/2;
  #local f3 = Fn(x3);
  #if (f1*f2 > 0)
    #error "Macro Bisection: Endpoints does not bracket the root"
  #else
    #while (abs(f3) > Epsilon)
      #if (f1*f3 < 0)
        #local x2 = x3;
        #local f2 = f3;
      #else
        #if (f2*f3 < 0)
          #local x1 = x3;
          #local f1 = f3;
        #else
          #error "Macro Bisection: Maybe Epsilion is negative ?"
        #end // if
      #end // if
      #local x3 = (x1 + x2)/2;
      #local f3 = Fn(x3);
      #debug ":"
    #end // while
  #end // if
  
  (x3)

#end // macro Bisection
        
// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 =
// Sample usage

#declare Fn = function(x) { 3*x^3 - 5*x^2 + x - 1 }

#declare CR = 0.000001; // Convergence criterion

#debug "\n"
#debug "\n"
#declare Root = NewtonRaphson(function(x) { Fn(x) }, -2, CR);
#debug str(Root, 0, -1)
#debug "\n"
#declare Root = RegulaFalsi(function(x) { Fn(x) }, -2, 3, CR);
#debug str(Root, 0, -1)
#debug "\n"
#declare Root = Bisection(function(x) { Fn(x) }, -2, 3, CR);
#debug str(Root, 0, -1)
#debug "\n"
#debug "\n"

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 =

Post a reply to this message

>
> I have a suspicion that the RegulaFalsi
> method is the best candidate for the job.
>

Oh yeah, that's what I was thinking :)

Seriously Tor, thanks for the efforts, but can you help a guy who was
educated in Design understand what RegulaFalsi is?  What does it mean to
solve something "numerically"?

Thanks,
Arie

Post a reply to this message

"Arie L. Stavchansky" wrote:

> Seriously Tor, thanks for the efforts, but can you help a guy
> who was educated in Design understand what RegulaFalsi is? 
> What does it mean to solve something "numerically"?

To abandon the search for a symbolic expression like "pi-ln(3)" and
settle for "well, it's between 2.0429 and 2.0430"; to try a rough
solution, see how far off it is, adjust proportionately and try again.

There are several versions of `adjust proportionately', applicable to
different kinds of problems.  In this case I'd use the secant method
(which is similar to regula falsi). 

Start with
	x0 = the lowest reasonable value of x (zero?)
	x1 = the highest reasonable value of x (r1?)
Do this repeatedly:
	q0 = x0 - (r1 + (R - x0)/(R * r1)) * sin(atan(x0/p(x0))
	q1 = x1 - (r1 + (R - x1)/(R * r1)) * sin(atan(x1/p(x1))
	x2 = x1 - (x1-x0)*q1/(q1-q0)	//the fun bit
	x0 = x1; x1 = x2
until abs(q1) is small enough.

The line I've marked as "the fun bit" amounts to imagining that <x0,q0>
and <x1,q1> are plotted on a graph, drawing a line through them and
marking where it crosses the axis.  If q is well-behaved, the
crossing-point (x2) is a better guess than x0 or x1, so throw away the
oldest guess (x0).


-- 
Anton Sherwood

Post a reply to this message

"Arie L. Stavchansky" wrote:
> While muddling around, I have come up with this.  You can start
> see the effect I am after, but this is truly a linear inclination. 
> I am looking to acheive a curved inclination all the while my
> spheres remain tangential.  I just can't get enough of this stuff.

Hooked, eh?

If you want to keep the spheres one size (or roughly so),
you'll have to abandon the equal spacing of the green lines.

	r0 = radius of the spheres of the previous tier
	R0 = radius of the circle containing their centers
	n0 = how many are in that tier

	p(R) = your profile curve
	P(R) = its first derivative

In the next tier, as a first approximation,
	r1 = r0
	R1 = R0 - 2*r0 / sqrt(1+(P(R0))^2)	// i hope; i'm tired!
	n1 = round(n0*R1/R0)

This won't be an exact fit, but there are techniques for improving it
(which, unfortunately, are beyond my intelligence at this hour).


...
What is your desired profile curve, by the way?
For a different look, I might try adapting what others have called a
"Fibonacci spiral" (really the Fermat spiral with golden sector
sampling).


-- 
Anton Sherwood

Post a reply to this message

Arie, I have thought about how to explain
what numerical methods are, but I find it
difficult to give a good explanation.

But here's my try at it anyway:

When one are solving a mathematical problem
analytically, one often manipulate SYMBOLS
that are describing the problem. The mani-
pulations are done according to a set of 
rules that preserves the exact description
of the problem. (And by doing it this way,
any found solutions will be exact too.)

If the manipulations are done in an intelli-
gent way, each step of them will simplify the
problem, or parts of it, until simple
expressions containing the symbols are found.

But if one are solving a problem numerically,
one may start with a symbolic description of
the problem. And, sometimes, one thereafter
rewrites (and maybe simplifies) that
description a bit. Then one substitute
NUMBERS for some or all of the symbols in
the (modified) description and evaluate it.
(These numbers can be results of analysis,
guessing or measurements obtained from "the
real world".)

Often one can use the result of this
evaluation to determine how "far off" those
initial numbers are from a solution (i.e.
estimate the errors).

Some knowledge of similar problems can then
often be used to try to guess or estimate how
to alter those numbers in order to get closer
to a solution. If so, one may repeat the
process all over again and hope that a good-
enough solution appears (i.e. has an
acceptable error.)

If one find a solution to a problem, one
sometimes have enough information
to get closer to other solutions analytically.
If one succeeds doing this, one can then go
back and apply further numerical methods to
seek out other solutions (or more accurate
ones).

It is important to notice that with numerical
methods one does not get exact results if one
are using computers. This is because the
numbers involved are always either truncated
or rounded off. And also because some of the
calculations the computers perform are done
by numerical methods themselves or by
approximations of the mathematical operations.


I'll try to explain more about the mentioned
methods later.


Tor Olav

Btw:
Arie, did you read my reply to your question in
the "Not a sphere" thread (by Zebu 2. Oct.) ?


Here's a web page that you hopefully will
find interesting:

http://www.damtp.cam.ac.uk/user/fdl/people/sd/lectures/index.html

It contains links to some "Numerical Methods"
lecture notes by Stuart Dalziel.

There's a html version here:
http://www.damtp.cam.ac.uk/user/fdl/people/sd/lectures/nummeth98/index.htm

And if you look at this page;
http://www.damtp.cam.ac.uk/user/fdl/people/sd/lectures/nummeth98/roots.htm

- you'll find these sub-chapters(?):

"3.3 Linear interpolation (regula falsi)"
(Read the part "3.2 Bisection" first.)

"3.5 Secant (chord)"
(Read the part "3.4 Newton-Raphson" first.

- where he talks about the methods I and
Anton has mentioned.

It is often quite useful to study the graph-
images supplied within, while you are
figuring out what is _really_ going on.

Post a reply to this message

Hi Tor,

Thank you so much for your help and response.  I feel like this group
provides kinda of like an "underground" educational curriculum.  It's kinda
nice knowing that science and mathematics have something to do with the
visual arts.

I did read your answer to my question on "Not a sphere. . ." and understood
it.  Thank you for that response as well.  I have visited the sites you have
given me and will say that (for me) it will take a while to absorb the
material.  Thank you for answering my questions, it is greatly appreciated.

Arie

in article 3BD8B1DF.82914778@hotmail.com, Tor Olav Kristensen at
tor### [at] hotmailcom wrote on 10/25/01 7.44 PM:

> 
> Arie, I have thought about how to explain
> what numerical methods are, but I find it
> difficult to give a good explanation.
> 
> But here's my try at it anyway:
> 
> When one are solving a mathematical problem
> analytically, one often manipulate SYMBOLS
> that are describing the problem. The mani-
> pulations are done according to a set of
> rules that preserves the exact description
> of the problem. (And by doing it this way,
> any found solutions will be exact too.)
> 
> If the manipulations are done in an intelli-
> gent way, each step of them will simplify the
> problem, or parts of it, until simple
> expressions containing the symbols are found.
> 
> But if one are solving a problem numerically,
> one may start with a symbolic description of
> the problem. And, sometimes, one thereafter
> rewrites (and maybe simplifies) that
> description a bit. Then one substitute
> NUMBERS for some or all of the symbols in
> the (modified) description and evaluate it.
> (These numbers can be results of analysis,
> guessing or measurements obtained from "the
> real world".)
> 
> Often one can use the result of this
> evaluation to determine how "far off" those
> initial numbers are from a solution (i.e.
> estimate the errors).
> 
> Some knowledge of similar problems can then
> often be used to try to guess or estimate how
> to alter those numbers in order to get closer
> to a solution. If so, one may repeat the
> process all over again and hope that a good-
> enough solution appears (i.e. has an
> acceptable error.)
> 
> If one find a solution to a problem, one
> sometimes have enough information
> to get closer to other solutions analytically.
> If one succeeds doing this, one can then go
> back and apply further numerical methods to
> seek out other solutions (or more accurate
> ones).
> 
> It is important to notice that with numerical
> methods one does not get exact results if one
> are using computers. This is because the
> numbers involved are always either truncated
> or rounded off. And also because some of the
> calculations the computers perform are done
> by numerical methods themselves or by
> approximations of the mathematical operations.
> 
> 
> I'll try to explain more about the mentioned
> methods later.
> 
> 
> Tor Olav
> 
> Btw:
> Arie, did you read my reply to your question in
> the "Not a sphere" thread (by Zebu 2. Oct.) ?
> 
> 
> Here's a web page that you hopefully will
> find interesting:
> 
> http://www.damtp.cam.ac.uk/user/fdl/people/sd/lectures/index.html
> 
> It contains links to some "Numerical Methods"
> lecture notes by Stuart Dalziel.
> 
> There's a html version here:
> http://www.damtp.cam.ac.uk/user/fdl/people/sd/lectures/nummeth98/index.htm
> 
> And if you look at this page;
> http://www.damtp.cam.ac.uk/user/fdl/people/sd/lectures/nummeth98/roots.htm
> 
> - you'll find these sub-chapters(?):
> 
> "3.3 Linear interpolation (regula falsi)"
> (Read the part "3.2 Bisection" first.)
> 
> "3.5 Secant (chord)"
> (Read the part "3.4 Newton-Raphson" first.
> 
> - where he talks about the methods I and
> Anton has mentioned.
> 
> It is often quite useful to study the graph-
> images supplied within, while you are
> figuring out what is _really_ going on.

Post a reply to this message