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