# The Pegasus and related methods for solving equations

Share on:

In the previous post, we saw that a small change to the method of false position provided much faster convergence, while retaining its bracketing.

This was the Illinois method which is only one of a whole host of similar methods, some of which converge even faster.

And as a reminder, here's its definition, with a very slight change:

Given $$x_{i-1}$$ and $$x_i$$ that bracket a root and their function values $$f_{i-1}$$, $$f_i$$, first compute the secant value

$x_{i+1}=\frac{x_{i-1} f_i - x_i f_{i-1}}{f_i - f_{i-1}}.$

and let $$f_{i+1}=f(x_{i+1})$$. Then:

1. if $$f_if_{i+1}<0$$, replace $$(x_{i-1},f_{i-1})$$ with $$(x_i,f_i)$$
2. if $$f_if_{i+1}>0$$, replace $$(x_{i-1},f_{i-1})$$ with $$(x_{i-1},\gamma f_{i-1})$$ with $$\gamma=0.5$$.

In each case we replace $$(x_i,f_i)$$ by $$(x_{i+1},f_{i+1})$$.

Much research since has been investigating possible scaling values for $$\gamma$$. If $$\gamma$$ is to be constant, then it can be shown that $$\gamma=0.5$$ is optimal. But $$\gamma$$ need not be constant.

## The Pegasus method

This was defined by Dowell & Jarratt, whose form of the Illinois method we used in the last post; see their article "The 'Pegasus' method for computing the root of an equation." BIT Numerical Mathematics 12 (1972),pp503-508.

Here we use

$\gamma = \frac{f_i}{f_i+f_{i+1}}$

And here's the Julia function for defining it; basically the same as the Illinois function of the previous post, with the differences both of $$\gamma$$ and of showing the absolute difference between successive iterations:

 1function pegasus(f,a,b;num_iter = 20)
2    if f(a)*f(b) > 0
3	error("Values given are not guaranteed to bracket a root")
4    else
5	fa, fb = f(a), f(b)
6	c = b
7	for k in 1:num_iter
8	    c_old = c
9	    c = b - fb*(b-a)/(fb-fa)
10	    fc = f(c)
11	    if fb*fc < 0
12		a, fa = b, fb
13	    else
14		fa = fa*fb/(fb+fc)
15	    end
16	    b, fb = c, fc
17	    @printf("%2d: %1.60f, %1.15e\n",k,c,abs(c_old-c))
18	    #println(c,",  ",abs(c_old-c))
19	end
20    end
21end


And with the same function $$f(x)=x^5-2$$:

 1Julia> pegasus(f,BigFloat("2"),BigFloat("3"))
2
3 1: 1.032258064516129032258064516129032258064516129032258064516128, 9.677419354838710e-01
4 2: 1.058249216160286723536401978566436704093370062081736874593804, 2.599115164415769e-02
5 3: 1.095035652659330505424147240084976534578418952240763525441535, 3.678643649904378e-02
6 4: 1.131485704080638653175790037904708402277455906284414289399986, 3.645005142130815e-02
7 5: 1.147884687198048718506398066361222614002745776909137282797727, 1.639898311741007e-02
8 6: 1.148720321893174344989720370927796480725850707839146637042703, 8.356346951256265e-04
9 7: 1.148698323855563082475350143443841164951825177665336110156093, 2.199803761126251e-05
10 8: 1.148698354995843974508573437584337692570858791644071608835544, 3.114028089203322e-08
11 9: 1.148698354997035006796886004382924209539326506468986782249386, 1.191032288312567e-12
1210: 1.148698354997035006798626946777931199507739956238584846007635, 1.740942395006990e-21
1311: 1.148698354997035006798626946777927589443850889097797494571041, 3.610063889067141e-33
1412: 1.148698354997035006798626946777927589443850889097797505513712, 1.094267079108742e-53
1513: 1.148698354997035006798626946777927589443850889097797505513712, 0.000000000000000e+00
1614: 1.148698354997035006798626946777927589443850889097797505513711, 1.244603055572228e-60
1715: 1.148698354997035006798626946777927589443850889097797505513711, 0.000000000000000e+00
1816: 1.148698354997035006798626946777927589443850889097797505513712, 1.244603055572228e-60
1917: 1.148698354997035006798626946777927589443850889097797505513711, 1.244603055572228e-60
2018: 1.148698354997035006798626946777927589443850889097797505513711, 0.000000000000000e+00
2119: 1.148698354997035006798626946777927589443850889097797505513712, 1.244603055572228e-60
2220: 1.148698354997035006798626946777927589443850889097797505513711, 1.244603055572228e-60


This is indeed faster than the Illinois method, with an efficiency index of $E≈ 1.64232$\$.

For another example, compute that value of Lambert's W function $$W(100)$$; this is the solution of $$xe^x-100=0$$; Lambert's function is the inverse of $$y=xe^x$$.

 1Julia> pegasus(x->x*exp(x)-100,BigFloat("3"),BigFloat("4"),num_iter=11)
2
3 1: 3.251324125460162218273541775021703846514592591436893826064677, 7.486758745398378e-01
4 2: 3.340634428196726809051441645629843725138212599102446141079158, 8.931030273656459e-02
5 3: 3.380778785367380815168698736373250080216344738301354721341292, 4.014435717065401e-02
6 4: 3.385665875268764090984779722929318950325131118803147296896820, 4.887089901383276e-03
7 5: 3.385630033731458485133511320033990390684414498737973972888242, 3.584153730560585e-05
8 6: 3.385630140287712138407043176029520651934410763166067578057173, 1.065562536532735e-07
9 7: 3.385630140290050184887119964591950428421045978948047751759305, 2.338046480076789e-12
10 8: 3.385630140290050184888244364529728481623112988246686294800948, 1.124399937778053e-21
11 9: 3.385630140290050184888244364529726867491694170157806679271792, 1.614131418818089e-33
1210: 3.385630140290050184888244364529726867491694170157806680386175, 1.114382727073817e-54
1311: 3.385630140290050184888244364529726867491694170157806680386175, 0.000000000000000e+00


## Other similar methods

These differ only in the definition of the scaling value $$\gamma$$; several are given by J. A. Ford in "Improved Algorithms of Illinois-Type for the Numerical Solution of Nonlinear Equations" University of Essex, Department of Computer Science (1995). This article is happily available online.

Following Ford, define

$\phi_k=\frac{f_{k+1}}{f_k}$

and then:

Method $$\gamma$$ Efficiency Index
Anderson & Björck if $$f_i>f_{i+1}$$ then $$1-\phi_i$$ else 0.5 1.70998 or 1.68179
Ford method 1 $$(1-\phi_i-\phi_{i-1})/(1+\phi_i-\phi_{i-1})$$ 1.55113
Ford method 2 $$(1-\phi_i)/(1-\phi_{i-1})$$ 1.61803
Ford method 3 $$1-(\phi_i/(1-\phi_{i-1}))$$ 1.70998
Ford method 4 $$1-\phi_i-\phi_{i-1}$$ 1.68179
Ford method 5 $$(1-\phi_i)/(1+\phi_i-\phi_{i-1})$$ not given

The efficiency of the Anderson-Björck method depends on whether the sign of

$K = \left(\frac{c_2}{c_1}\right)^{\!2}-\frac{c_3}{c_1}$

is positive or negative, where

$c_k = \frac{f^{(k)}(x^*)}{k!},\;k\ge 1$

and $$x^*$$ is the solution. Note that the $$c_k$$ values are simply the coefficients of the Taylor series expansion of $$f(x)$$ about the root $$x^*$$; that is

$f(x-x^*) = c_1x+c_2x^2+c_3x^3+\cdots$