Nonlinear equations Numerical methods in R Systems of nonlinear equations Bi7740: Scientific computing Root finding Vlad Popovici popovici@iba.muni.cz Institute of Biostatistics and Analyses Masaryk University, Brno Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Outline 1 Nonlinear equations 2 Numerical methods in R 3 Systems of nonlinear equations Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Nonlinear equations scalar problem: f : R → R, find x ∈ R such that f(x) = 0 vectorial problem: f : Rn → Rn , fing x ∈ Rn such that f(x) = 0 in any case, here we consider f to be continuously differentiable everywhere in the neighborhood of the solution an interval [a, b] is a bracket for the function f if f(a)f(b) < 0 f continuous → f([a, b]) is an interval Bolzano’s thm.: if [a, b] is a bracket for f than there exists at least one x∗ ∈ [a, b] s.t. f(x∗) = 0 if f(x∗) = f (x∗) = · · · = f(m−1)(x∗) = 0 but f(m) 0 then x∗ has multiplicity m note: in Rn things are much more complicated Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Conditioning absolute condition number is for a scalar problem 1/|f (x∗)| a multiple is ill-conditioned absolute condition number for a vectorial problem is J−1 f (x∗) , where Jf is the Jacobian matrix of f, [Jf (x)]ij = ∂fi(x) ∂xj if the Jacobian is nearly singular, the problem is ill-conditioned Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Sensitivity and conditioning possible interpretations of the approximate solution: f(ˆx) − f(x∗ ) ≤ : small residual ˆx − x∗ ≤ closeness to the true solution the two criteria might not be satistfied simultaneously if the problem is well-conditioned: small residual implies accurate solution Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Convergence rate usually, the solution is find iteratively let ek = xk − x∗ be the error at the k−th iteration, where xk is the approximation and x∗ is the true solution the method converges with rate r if lim k→∞ ek+1 ek r = C, for C > 0 finite if the method is based on improving the bracketing, then ek = bk − ak if r = 1 and C < 1, the convergence is linear and a constant number of digits are "gained" per iteration r = 2 the convergence is quadratic, the number of exact digits doubles at each iteration r > 1 the converges is superlinear, increasing number of digits are gained (depends on r) Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Outline 1 Nonlinear equations 2 Numerical methods in R 3 Systems of nonlinear equations Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Bisection method Idea: refine the bracketing of the solution until the length of the interval is small enough. Assumption: there is only one solution in the interval. Algorithm 1: Bisection while (b − a) > do m ← a + b−a 2 if sign(f(a)) = sign(f(m)) then a ← m else b ← m a bm f(a) f(b) HOMEWORK Implement the above procedure in Matlab. Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Bisection, cont’d convergence is certain, but slow convergence rate is linear (r = 1 and C = 1/2) after k iterations, the length of the interval is (b − a)/2k , so achieving a tolerance requires log2 b − a iterations, idependently of f. the value of the function is not used, just the sign Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Fixed-point methods a fixed point for a function g : R → R is a value x ∈ R such that f(x) = x the fixed-point iteration xk+1 = g(xk ) is used to build a series of successive approximations to the solution for a given equation f(x) there might be several equivalent fixed-point problems x = g(x) Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations The solutions of the equation x2 − x − 2 = 0 are the fixed points of each of the following functions: g(x) = x2 − 2 g(x) = √ x + 2 g(x) = 1 + 2/x g(x) = x2 +2 2x−1 2 3 1 1 2 3 g(x)=1+2/x Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Conditions for convergence a function g : S ⊂ R → R is called Lipschitz-bounded if ∃α ∈ [0, 1] so that |f(x1) − f(x0)| ≤ α|x1 − x0|, ∀x0, x1 ∈ S in other words, if |g (x∗)| < 1, then g is Lipschitz-bounded for such functions, there exists an interval containing x∗ s.t. iteration xk+1 = g(xk ) converges to x∗ if started within that interval if |g (x∗)| > 1 the iterations diverge in general, convergence is linear smoothed iterations: xk+1 = λk xk + (1 − λk )f(xk ) with λk ∈ [0, 1] and limk→∞ λk = 0 Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Stopping criteria If either 1 |xk+1 − xk | ≤ 1|xk+1| (relative error) 2 |xk+1 − xk | ≤ 2 (absolute iteration error) 3 |f(xk+1) − f(xk )| ≤ 3 (absolute functional error) stop the iterations. Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Newton-Raphson method from Taylor series: f(x + h) ≈ f(x) + f (x)h so in a small neighborhood around x f(x) can be approximated by a linear function of h with the root −f(x)/f (x) Newton iteration: xk+1 = xk − f(xk ) f (xk ) x_kx_k+1 Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Implement the Newton-Raphson procedure in Matlab. function [r, xk] = NRiter(f, fprime, x, tol, max_iter) xk(1) = x; for k=1:max_iter xk(k+1) = ................ if break end r = xk(...); return Apply it like: >> f = @(x) (x^2 −4*sin(x)); >> fp = @(x) (2*x−4*cos(x)); >> NRiter(f, fp, 3) ans = 1.9338 Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Newton-Raphson method, cont’d convergence for a simple root is quadratic to converge, the procedure needs to start close enough to the solution, where the function f is monotonic Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Secant method (lat.: Regula falsi) secant method approximates the derivative by finite differences: xk+1 = xk − f(xk ) xk − xk−1 f(xk ) − f(xk−1) convergence is normally superlinear, with r ≈ 1.618 it must be started in a properly chosen neighborhood implement in Matlab [r,xk] = secant(f, x0, x1,...) x_k+1 x_k x_k−1 Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Interpolation methods and other approaches secant method uses linear interpolation one can use higher-degree polynomial interpolation (e.g. quadratic) and find the roots of the interpolating polynomial inverse interpolation: xk+1 = p−1 (yk ) where p is an interpolating polynomial fractional interpolation special methods for finding roots of the polynomials Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Fractional interpolation previous methods have difficulties with functions having horizontal or vertical asymptotes linear fractional interpolation uses φ(x) = x − u vx − w function, which has a vertical asymptote at x = w/v, a horizontal asymptote at y = 1/v and a zero at x = u Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Fractional interpolation, cont’d let x0, x1, x2 be three points where the function is estimates, yielding f0, f1, f2 find u, v, w for φ by solving   1 x0f0 −f0 1 x1f1 −f1 1 x2f2 −f2     u v w   =   x0 x1 x2   the iteration step swaps the values: x0 ← x1 and x1 ← x2 the new approximate solution is the zero of the linear fraction, x2 = u. This can be implemented as x2 ← x2 + (x0 − x2)(x1 − x2)(f0 − f1)f2 (x0 − x2)(f2 − f1)f0 − (x1 − x2)(f2 − f0)f1 Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Outline 1 Nonlinear equations 2 Numerical methods in R 3 Systems of nonlinear equations Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Systems of nonlinear equations much more difficult than the scalar case no simple way to ensure convergence computational overhead increases rapidly with the dimension difficult to determine the number of solutions difficult to find a good starting approximation Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Fixed-point methods in Rn g : Rn → Rn , x = g(x) = [g1(x), . . . , gn(x)] fixed-point iteration: xk+1 = g(xk ) denote ρ(Jg(x)) the spectral radius (maximum absolute eigenvalue) of the Jacobian matrix of g evaluated at x if ρ(Jg(x∗)) < 1, the fixed point iteration converges if started close enough to the solution the convergence is linear with C = ρ(Jg(x∗)) Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Newton-Raphson method in Rn xk+1 = xk − J−1 f (xk )f(xk ) no need for inversion; solve the system Jf (xk )sk = −f(xk ) for Newton step sk and iterate xk+1 = xk + sk HOMEWORK: Implement in Matlab the above procedure Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Broyden’s method uses approximations of the Jacobian the initial approximation of J can be the actual Jacobian (if available) or even I Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Algorithm 2: Broyden’s method for k = 0, 1, 2, . . . do solve Bk sk = −f(xk ) for sk xk+1 = xk + sk yk = f(xk+1) − f(xk ) Bk+1 = Bk + ((yk − Bk sk )sT k )/(sT k sk ) if xk+1 − xk ≥ 1(1 + xk+1 ) then continue if f(xk+1) < 2 then x∗ = xk+1 break else algorithm failed Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Further topics secant method is also extended to Rn (see Boryden’s method) robust Newton-like methods: enlarge the region of convergence, introduce a scalar parameter to ensure progression toward solution in Matlab: roots() for roots of a polynomial; fzero() for scalar case; fsolve() for the Rn case. Vlad Bi7740: Scientific computing Nonlinear equations Numerical methods in R Systems of nonlinear equations Matlab exercise help fsolve try the examples from the documentation solve the system:    x2 + 2y2 − 5x + 7y = 40 3x2 − y2 + 4x + 2y = 28 with initial approximation [2, 3] Vlad Bi7740: Scientific computing