Bairstow's method, Durand-Kerner method and other root-finding methods

  • 103 Replies
  • 2047 Views

0 Members and 1 Guest are viewing this topic.

Offline ottomagus

  • *
  • Fractal Phenom
  • ****
  • Posts: 56
  • Otto Magus
    • Otto Magus Digital Art
« Reply #90 on: August 10, 2020, 12:37:29 PM »
I couldn't make much headway with the Bisection method. In the paper linked to by FractalAlex, below, there is something called the Midpoint method. This does, at least, have a kind of similar name :). I've played around with this one before, it gives shapes reminiscent of centipedes. See first image below. The formula is z=z-f/(f1*(z-(f/f1)/2)), where f is the expression of z (z^3-1 here) and f1 the first derivative.

The other images are of a Newton/Midpoint hybrid. Newton iteration alternating with Midpoint. The first one starts with Newton, the second with Midpoint.

Edit: The formula given above had one too many brackets. Now corrected.
« Last Edit: August 12, 2020, 01:04:04 PM by ottomagus »

Offline ottomagus

  • *
  • Fractal Phenom
  • ****
  • Posts: 56
  • Otto Magus
    • Otto Magus Digital Art
« Reply #91 on: August 11, 2020, 09:26:19 AM »
Here's two more, the Nova version of the Midpoint method. Just adding +pixel to the UF code. The first is the unzoomed fractal, the second is rotated and zoomed to show the detail better.

Offline FractalAlex

  • *
  • Fractal Flamingo
  • ****
  • Posts: 344
  • Experienced root-finding method expert
« Reply #92 on: August 11, 2020, 04:46:16 PM »
I love it! Would you mind sharing the parameters? Thanks, ottomagus!
"I am lightning, the rain transformed."
- Raiden, Metal Gear Solid 4: Guns of the Patriots

Offline ottomagus

  • *
  • Fractal Phenom
  • ****
  • Posts: 56
  • Otto Magus
    • Otto Magus Digital Art
« Reply #93 on: August 12, 2020, 01:20:20 AM »
Glad you like it, Alex!

Parameters are here https://onedrive.live.com/?id=10065AB7502DF479%2152982&cid=10065AB7502DF479

You may have to update your public formulas, as one formula was uploaded to the database yesterday...

Offline ottomagus

  • *
  • Fractal Phenom
  • ****
  • Posts: 56
  • Otto Magus
    • Otto Magus Digital Art
« Reply #94 on: August 13, 2020, 09:36:22 PM »
Schröder's method:
z=z-f*f1/(f1^2-f*f2)
where f is an expression of z, f1 is the first derivative and f2 is the second derivative.

This can be formed by applying Newton's method to f/f1.
According to WolframMathWorld "Almost all other root-finding methods can be considered as special cases of Schröder's method".

First image: Schröder's method for f=(z^2+z)^3-1
Second image: Schröder's method for f=z^3+z^2+z+(z^2+z)^3-1

Parameters available via link in my previous post.

Some Schröder hybrids to follow....

Offline FractalAlex

  • *
  • Fractal Flamingo
  • ****
  • Posts: 344
  • Experienced root-finding method expert
« Reply #95 on: August 13, 2020, 10:06:18 PM »
Ah, yes. Schroder's method. When the function is z^3 - 1 = 0, you see structures that are like the standard Newton's method. The formula is similar to Halley's method's formula, but with the 2's in the nominator and denominator of the Newton step. This slows the order of convergence to 2, which gives more interesting patterns and fractals, unlike Halley's method. The reason is that Schroder's method has more extremum points than Halley's method. Also, if you turn it into a Nova fractal, you'll see that the order of convergence is equal to the power of the minibrots. For Newton, Schroder and Steffensen, it is equal to 2. For cubic minibrots, it's Halley, Householder, Parhalley and Laguerre's methods. For quartic minibrots, it's the third order Householder method and King's family of iterative multipoints.

Offline ottomagus

  • *
  • Fractal Phenom
  • ****
  • Posts: 56
  • Otto Magus
    • Otto Magus Digital Art
« Reply #96 on: August 14, 2020, 12:45:11 PM »
I have experimented with all the formulas you mention, Alex, and agree that Schröder gives much more interesting shapes than Halley. I didn't know why, though! Thanks!

Offline ottomagus

  • *
  • Fractal Phenom
  • ****
  • Posts: 56
  • Otto Magus
    • Otto Magus Digital Art
« Reply #97 on: August 14, 2020, 01:02:47 PM »
The first two images combine Schröder and Midpoint.
The next two combine Schröder for f=z^3-1 and Schröder for f=z^3-z^2-z-1. I'm calling this a semi-hybrid :)

Offline ottomagus

  • *
  • Fractal Phenom
  • ****
  • Posts: 56
  • Otto Magus
    • Otto Magus Digital Art
« Reply #98 on: August 15, 2020, 01:04:33 PM »
Jarratt's Method is one I had not tried before. The formula (using standard notation this time) is:

u(z) = f(z)/f'(z)
z = z - u(z)/2 + f(z)/(f'(z)-3*f'(z-2*u(z)/3))

These three images utilize the "relaxation parameter" a, so the formula becomes:

u(z) = f(z)/f'(z)
z = z - a*(u(z)/2 - f(z)/(f'(z)-3*f'(z-2*u(z)/3)))

Offline FractalAlex

  • *
  • Fractal Flamingo
  • ****
  • Posts: 344
  • Experienced root-finding method expert
« Reply #99 on: August 15, 2020, 01:29:06 PM »
Nice! Also, Oscar Veliz will cover this method in one of his future videos. He recently covered videos about minimization, including ternary search, dichotomous search, Fibonacci search and golden-section search methods. After Jarratt's method, he will cover Brent's minimization method, and plans to do a two-parter about the Jenkins-Traub algorithm.

Offline FractalAlex

  • *
  • Fractal Flamingo
  • ****
  • Posts: 344
  • Experienced root-finding method expert
« Reply #100 on: August 15, 2020, 07:41:30 PM »
Does anyone know how to turn the following:
Code: [Select]
-----------------------------------------
-- Example Newton-Bisection Hybrid
-- Root-Finding Method created by
-- Oscar Veliz based off of NewtSafe from
-- https://www.ldeo.columbia.edu/~mspieg/e4300/BlankPDFs/Lecture06_blank.pdf
-- and rtsafe from Numerical Recipies in C p.366-367 and
-- Numerical Methods That Work by Acton p.51-52
-- www.youtube.com/OscarVeliz
-----------------------------------------
with Ada.Text_IO; use  Ada.Text_IO;
with Ada.Integer_Text_IO; use  Ada.Integer_Text_IO;
with Ada.Float_Text_IO; use  Ada.Float_Text_IO;
with Ada.Task_Identification;  use Ada.Task_Identification;
with Ada.Numerics.Generic_Elementary_Functions;

procedure NewtSafe is
---------------------
-- f(x) = x^3 - x^2 - x - 1
---------------------
function f (x: Float) return Float is
begin
return x**3 - x**2 - x - 1.0;
end;
-------------------
-- f'(x) = 3x^2 - 2x - 1
-------------------
function fp (x: Float) return Float is
begin
return 3.0*x*x - 2.0*x - 1.0;
end;
----------
-- Regular Newton's Method (unsafe)
----------
function newt (xinit: Float) return Float is
x: Float := xinit;
eps: Float := 10.0**(-6);
begin
while abs(f(x)) > eps loop
put(x,2,6,0);put_line("");
x := x - f(x) / fp(x);
end loop;
return x;
end;
-- header so forceBis can call hybrid
function hybrid (inita, initb: Float) return Float;
----------
-- Bisection between a and b then restart Hybrid
----------
function forceBis(a, b: Float) return Float is
c: Float := (a + b) / 2.0;
begin
--put_line("forced bisection");
if f(a) * f(c) < 0.0 then -- fa & fc have different signs
return hybrid(a, c);
else
return hybrid(c, b);
end if;
end;
----------
-- Newton-Bisection Hybrid, inita and initb specify interval
----------
function hybrid (inita, initb: Float) return Float is
x, fa, fb, fx, fpx, oldx, lastStep: Float;
a: Float := inita;
b: Float := initb;
eps: Float := 10.0**(-6);
begin
if(a > b) then --ensure a is on left, swap
a := initb;
b := inita;
end if;
if a = b then
return a;
end if;
-- check interval points
fa := f(a);
fb := f(b);
if fa = 0.0 then
return a;
end if;
if fb = 0.0 then
return b;
end if;
if fa * fb > 0.0 then -- invalid interval
put_line("error: interval must contain root.");
Abort_Task (Current_Task);
end if;
lastStep := b - a;
oldx := b; -- intitial oldx is on boundary
x := (a + b) / 2.0; -- start from midpoint
fx := f(x);
while abs(fx) > eps loop
put(x,2,6,0);put_line("");
fpx := fp(x);
if fpx = 0.0 then -- avoid divide by zero
return forceBis(a,b);
end if;
lastStep := abs(x - oldx);
if abs(fx * 2.0) > abs(lastStep * fpx) then -- decreasing slowly
--put_line("too slow");
return forceBis(a,b);
end if;
oldx := x;
x := x - fx / fpx;
if x = oldx then -- converged
return x;
end if;
if x <= a or else x >= b then -- avoid going outside interval
--put_line("outside interval");
return forceBis(a,b);
end if;
--shrink interval
fx := f(x);
if fx * fa < 0.0 then
b:= x;
fb:= fx;
else
a:= x;
fa:= fx;
end if;
end loop;
return x;
end;
begin --main--
put_line("Newton's Method");
put(newt(1.5),2,6,0); put_line("");
put_line("Newt-Safe");
put(hybrid(1.0,2.0),2,6,0); put_line("");
end NewtSafe;
into a viable Ultra Fractal formula? Some help would be very appreciated!

Offline FractalAlex

  • *
  • Fractal Flamingo
  • ****
  • Posts: 344
  • Experienced root-finding method expert
« Reply #101 on: August 15, 2020, 08:56:39 PM »
I managed to successfully make inverse quadratic interpolation work in Ultra Fractal today. The bad news is that there's no fractal structure, only the basins of attraction for the individual roots of the function (here is z^3 - z^2 - z - 1). Fun fact: it has the same order of convergence as Muller's method, being 1.84 (as can be proved by secant method analysis). That's because both methods rely on Lagrange polynomials.
Outside coloring used: continuous potential.
Code: [Select]
IQI3_1 {
init:
  z = pixel
  ooz = p2
  oz = p1
  i = sqrt(-1)
loop:
  f2 = ooz^3 - ooz^2 - ooz - 1
  f1 = oz^3 - oz^2 - oz - 1
  f = z^3 - z^2 - z - 1
  z = ooz*((f1*f)/((f2-f1)*(f2-f))) + oz*((f2*f)/((f1-f2)*(f1-f))) + z*((f2*f1)/((f1-f2)*(f-f1)))
  ooz = oz
  oz = z
bailout:
  |z - oz| > @bailout

default:
  title = "Inverse Quadratic Interpolation (z^3 - z^2 - z - 1)"
  maxiter = 1000

  param p1
    caption = "1st start value"
    default = 1.99999999999
    hint = "Start value for oz."
  endparam

  param p2
    caption = "2nd start value"
    default = 2.0
    hint = "Start value for ooz."
  endparam

  param bailout
    caption = "Bailout Value"
    default = 1.0e-20
    hint = "The bailout value."
  endparam
}
« Last Edit: August 15, 2020, 09:22:00 PM by FractalAlex »

Offline ottomagus

  • *
  • Fractal Phenom
  • ****
  • Posts: 56
  • Otto Magus
    • Otto Magus Digital Art
« Reply #102 on: August 17, 2020, 09:30:56 PM »
Alex. I have bookmarked the Oscar Veliz channel and will explore!

In my previous post, the Jarratt formula is right but I got the code wrong. An "accidental modification" that fortuitously turned out better than the "correct" fractal :).
Jarratt's method actually produces fairly generic results, similar to Traub-Ostrowski (in fact, for f=z^3-1 it is exactly the same as Traub-Ostrowski).

Here's a few more pics. These, like the others, are one layer UF images intended to show some of the basic qualities of the fractal.

The first two are a hybrid of two Whittaker methods:
z=z-f(z)*(2-f(z)*f''(z)/f'(z)^2 )/(2*f'(z))
and
a=f(z)*f''(z)/f'(z)^2
z=z-f(z)*(2-a+(4+2*a)/(2-a*(2-a)))/(4*f'(z))

The third is a modification of a formula attributed to Nedzhibov et al from https://file.scirp.org/pdf/AM20100300002_25669189.pdf
The modified formula is:
k=z-f(z)/f'(z)
z=k-k*f(k)/(f-a*f(k))
where a is a parameter

The last two are a Newton/Julia hybrid.

Parameters here https://onedrive.live.com/?id=10065AB7502DF479%2152982&cid=10065AB7502DF479

Offline ottomagus

  • *
  • Fractal Phenom
  • ****
  • Posts: 56
  • Otto Magus
    • Otto Magus Digital Art
« Reply #103 on: August 25, 2020, 06:37:13 PM »
Here's some UF code for a nova version of the secant method:

Code: [Select]
SecantNova {
init:
if @st==0
z1=@sec1*@nf(pixel)
z2=@sec2*@nf(pixel)
else
z1=@sec1
z2=@sec2
endif
loop:
if @v==0
fz1=z1^@pow-@c
fz2=z2^@pow-@c
elseif @v==1
fz1=z1^@pow-z1^(@pow-1)-z1^(@pow-2)-@c
fz2=z2^@pow-z2^(@pow-1)-z2^(@pow-2)-@c
elseif @v==2
fz1=z1^@pow+z1^(@pow-1)+z1^(@pow-2)-@c
fz2=z2^@pow+z2^(@pow-1)+z2^(@pow-2)-@c
elseif @v==3
fz1=z1^@pow*(z1-1.5)/(1-1.5*z1)-@c
fz2=z2^@pow*(z2-1.5)/(1-1.5*z2)-@c
elseif @v==4
fz1=z1^@pow*(z1-0.5)/(1-0.5*z1)-@c
fz2=z2^@pow*(z2-0.5)/(1-0.5*z2)-@c
elseif @v==5
fz1=z1^@pow+1/z1^@pow-@c
fz2=z2^@pow+1/z2^@pow-@c
elseif @v==6
fz1=z1^@pow-1/z1^@pow-@c
fz2=z2^@pow-1/z2^@pow-@c
elseif @v==7
fz1=((z1^2-1)/(z1+1))^@pow-@c
fz2=((z2^2-1)/(z2+1))^@pow-@c
elseif @v==8
fz1=z1^@pow-z1^(@pow-1)-z1^(@pow-2)-1/z1^@pow-@c
fz2=z2^@pow-z2^(@pow-1)-z2^(@pow-2)-1/z2^@pow-@c
elseif @v==9
fz1=z1^@pow+cotanh(z1)^@pow-@c
fz2=z2^@pow+cotanh(z2)^@pow-@c
elseif @v==10
fz1=z1^@pow-cotanh(z1)^@pow-@c
fz2=z2^@pow-cotanh(z2)^@pow-@c
elseif @v==11
fz1=((z1^@pow-1)/(@pow*z1^(@pow-1)))^@pow-@c
fz2=((z2^@pow-1)/(@pow*z2^(@pow-1)))^@pow-@c
elseif @v==12
fz1=(z1^@pow-1)/(z1^(@pow-1))-@c
fz2=(z2^@pow-1)/(z2^(@pow-1))-@c
elseif @v==13
fz1=(z1^2/(2*z1+1))^@pow-@c
fz2=(z2^2/(2*z2+1))^@pow-@c
elseif @v==14
fz1=((Z1^2-1)/(1-z1))^@pow-@c
fz2=((Z2^2-1)/(1-z2))^@pow-@c
else
fz1=((Z1^2-1.5)/(1-1.5*z1))^@pow-@c
fz2=((Z2^2-1.5)/(1-1.5*z2))^@pow-@c
endif
z=z2-@sec3*fz2*(z2-z1)/(fz2-fz1)
if @st==0
z1=z2+@np*@nf(pixel)
else
z1=z2+@nf(pixel)
endif
z2=z
fz1=fz2
if @v==0
fz2=z^@pow-@c
elseif @v==1
fz2=z^@pow-z^(@pow-1)-z^(@pow-2)-@c
elseif @v==2
fz2=z^@pow+z^(@pow-1)+z^(@pow-2)-@c
elseif @v==3
fz2=z^@pow*(z-1.5)/(1-1.5*z)-@c
elseif @v==4
fz2=z^@pow*(z-0.5)/(1-0.5*z)-@c
elseif @v==5
fz2=z^@pow+1/z^@pow-@c
elseif @v==6
fz2=z^@pow-1/z^@pow-@c
elseif @v==7
fz2=((z^2-1)/(z+1))^@pow-@c
elseif @v==8
fz2=z^@pow-z^(@pow-1)-z^(@pow-2)-1/z^@pow-@c
elseif @v==9
fz2=z^@pow+cotanh(z)^@pow-@c
elseif @v==10
fz2=z^@pow-cotanh(z)^@pow-@c
elseif @v==11
fz2=((z^@pow-1)/(@pow*z^(@pow-1)))^@pow-@c
elseif @v==12
fz2=(z^@pow-1)/(z^(@pow-1))-@c
elseif @v==13
fz2=(z^2/(2*z+1))^@pow-@c
elseif @v==14
fz2=((Z^2-1)/(1-z))^@pow-@c
else
fz2=((Z^2-1.5)/(1-1.5*z))^@pow-@c
endif
bailout:
|fz2|>@bail
default:
title="Secant Nova"
center=(0.001,0.001)
periodicity=0
param st
caption="Start Type"
enum="Pixel""Parameter"
endparam
param v
caption="Variant"
enum="z^p-c""z^p-z^(p-1)-z^(p-2)-c""z^p+z^(p-1)+z^(p-2)-c"\
"z^p*(z-2)/(1-2*z)-c""z^p*(z-0.5)/(1-0.5*z)-c"\
"z^p+1/z^p-c""z^p-1/z^p-c""((z^2-1)/(z+1))^p-c"\
"z^p-z^(p-1)-z^(p-2)-1/z^p-c""z^p+cotanh(z)^p-c"\
"z^p-cotanh(z)^p-c""((z^p-1)/(p*z^(p-1)))^p-c"\
"(z^p-1)/(z^(p-1))-c"
endparam
param sec1
caption="Secant Parameter 1"
default=(-0.7,0)
endparam
param sec2
caption="Secant Parameter 2"
default=(0.8,0)
endparam
param sec3
caption="Secant Parameter 3"
default=(1,0)
endparam
param pow
caption="Exponent"
default=(3,0)
endparam
param c
caption="Constant"
default=(1,0)
endparam
param np
caption="Nova Parameter"
default=(1,0)
visible=@st==0
endparam
func nf
caption="Nova Function"
default=ident()
endfunc
param bail
caption="Bailout"
default=0.000001
endparam
}


The "+pixel" part of the code doesn't go where I expected it to :)


... and some example images. Parameters here https://onedrive.live.com/?id=10065AB7502DF479%2153061&cid=10065AB7502DF479
Please note: These params use a plug-in version of the formula, which I only uploaded to the database today. UF users will probably need to update public formulas.


clip
Root finding via subdivision/IA

Started by marcm200 on Fractal Mathematics And New Theories

10 Replies
294 Views
Last post August 15, 2020, 06:41:45 AM
by marcm200
xx
Is it possible to speed up rendering with this method?

Started by Fluoroantimonic_Acid on Fractal Mathematics And New Theories

3 Replies
322 Views
Last post April 19, 2019, 12:59:29 AM
by gerrit
xx
Method of Early Checking Whether a Point Is In the Set

Started by Byte11 on Fractal Mathematics And New Theories

12 Replies
907 Views
Last post December 24, 2017, 12:33:10 AM
by Byte11
clip
Neat variation of Newton's method

Started by Ebanflo on Share a fractal

1 Replies
535 Views
Last post April 17, 2018, 04:34:13 PM
by mrrudewords
xx
Out of band method to contact admins?

Started by pauldelbrot on Forum Help And Support

6 Replies
478 Views
Last post September 04, 2018, 08:37:11 AM
by quaz0r