The Pegasus and related methods for solving equations
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:
- if \(f_if_{i+1}<0\), replace \((x_{i-1},f_{i-1})\) with \((x_i,f_i)\)
- 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 \]