|
![](/i/fill.gif) |
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] hotmail com
// 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/fill.gif) |