111,99 €
Mehr erfahren.
- Herausgeber: Wiley-VCH
- Kategorie: Fachliteratur
- Sprache: Englisch
This third book in a suite of four practical guides is an engineer's companion to using numerical methods for the solution of complex mathematical problems. The required software is provided by way of the freeware mathematical library BzzMath that is developed and maintained by the authors. The present volume focuses on optimization and nonlinear systems solution. The book describes numerical methods, innovative techniques and strategies that are all implemented in a well-established, freeware library. Each of these handy guides enables the reader to use and implement standard numerical tools for their work, explaining the theory behind the various functions and problem solvers, and showcasing applications in diverse scientific and engineering fields. Numerous examples, sample codes, programs and applications are proposed and discussed. The book teaches engineers and scientists how to use the latest and most powerful numerical methods for their daily work.
Sie lesen das E-Book in den Legimi-Apps auf:
Seitenzahl: 534
Veröffentlichungsjahr: 2013
Ähnliche
Contents
Cover
Related Titles
Title Page
Copyright
Preface
Chapter 1: Function Root-Finding
1.1 Introduction
1.2 Substitution Algorithms
1.3 Bolzano's Algorithm
1.4 Function Approximation
1.5 Use of a Multiprocessor Machine with a Known Interval of Uncertainty
1.6 Search for an Interval of Uncertainty
1.7 Stop Criteria
1.8 Classes for Function Root-Finding
1.9 Case Studies
1.10 Tests for BzzFunctionRoot and BzzFunctionRootMP Classes
1.11 Some Caveats
Chapter 2: One-Dimensional Optimization
2.1 Introduction
2.2 Measuring the Efficiency of the Search for the Minimum
2.3 Comparison Methods
2.4 Parabolic Interpolation
2.5 Cubic Interpolation
2.6 Gradient-Based Methods
2.7 Combination of Algorithms in a General Program
2.8 Parallel Computations
2.9 Search for the Interval of Uncertainty
2.10 Stop Criteria
2.11 Classes for One-Dimensional Minimization
2.12 Case Studies
Chapter 3: Unconstrained Optimization
3.1 Introduction
3.2 Heuristic Methods
3.3 Gradient-Based Methods
3.4 Conjugate Direction Methods
3.5 Newton's Method
3.6 Modified Newton Methods
3.7 Quasi-Newton Methods
3.8 Narrow Valley Effect
3.9 Stop Criteria
3.10 BzzMath Classes for Unconstrained Multidimensional Minimization
3.11 Case Studies
3.12 Tests
Chapter 4: Large-Scale Unconstrained Optimization
4.1 Introduction
4.2 Collecting a Sparse Symmetric Matrix
4.3 Ordering the Hessian Rows and Columns
4.4 Quadratic Functions
4.5 Hessian Evaluation
4.6 Newton's Method
4.7 Inexact Newton Methods
4.8 Practical Preconditioners
4.9 openMP Parallelization
4.10 Class for Large-Scale Unconstrained Minimization
Chapter 5: Robust Unconstrained Minimization
5.1 Introduction
5.2 One-Dimensional Minimization
5.3 Classes for One-Dimensional Robust Minimization
5.4 Examples in One-Dimensional Space
5.5 Examples in Multidimensional Space
5.6 Two-Dimensional Minimization
5.7 Classes for Robust Two-Dimensional Minimization
5.8 Examples for BzzMinimizationTwoVeryRobust Class
5.9 Multidimensional Robust Minimization
5.10 Class for Robust Multidimensional Minimization
Chapter 6: Robust Function Root-Finding
6.1 Introduction
6.2 Class and Examples
Chapter 7: Nonlinear Systems
7.1 Introduction
7.2 Comparing Nonlinear Systems to Other Iterative Problems
7.3 Convergence Test
7.4 Substitution Methods
7.5 Minimization Methods
7.6 Jacobian Evaluation
7.7 Newton's Method
7.8 Gauss–Newton Method
7.9 Modified Newton Methods
7.10 Newton's Method and Parallel Computations
7.11 Quasi-Newton Methods
7.12 Quasi-Newton Methods and Parallel Computing
7.13 Stop Criteria
7.14 Classes for Nonlinear System Solution with Dense Matrices
7.15 Tests for the BzzNonLinearSystem Class
7.16 Sparse and Large-Scale Systems
7.17 Large Linear System Solution with Iterative Methods
7.18 Classes for Nonlinear System Solution with Sparse Matrices
7.19 Continuation Methods
7.20 Solution of Certain Equations with Respect to Certain Variables
7.21 Case Studies
7.22 Special Cases
7.23 Some Caveats
Chapter 8: Underdimensioned Nonlinear Systems
8.1 Introduction
8.2 Underdimensioned Linear Systems
8.3 Class for Underdimensioned Nonlinear System Solution
Chapter 9: Constrained Minimization
9.1 Introduction
9.2 Equality Constraints
9.3 Equality and Inequality Constraints
9.4 Lagrangian Dual Problem
Chapter 10: Linear Programming
10.1 Introduction
10.2 Basic Attic Method Concepts
10.3 Attic Method
10.4 Differences between the Attic Method and Traditional Approaches
10.5 Explosion in the Number of Iterations
10.6 Degeneracy
10.7 Duality
10.8 General Considerations
Chapter 11: Quadratic Programming
11.1 Introduction
11.2 KKT Conditions for a QP Problem
11.3 Equality-Constrained QP
11.4 Equality- and Inequality-Constrained Problems
11.5 Class for QP
11.6 Projection or Reduced Direction Search Methods for Bound-Constrained Problems
11.7 Equality, Inequality, and Bound Constraints
11.8 Tests
Chapter 12: Constrained Minimization: Penalty and Barrier Functions
12.1 Introduction
12.2 Penalty Function Methods
12.3 Barrier Function Methods
12.4 Mixed Penalty–Barrier Function Methods
Chapter 13: Constrained Minimization: Active Set Methods
13.1 Introduction
13.2 Class for Constrained Minimization
13.3 Successive Linear Programming
13.4 Projection Methods
13.5 Reduced Direction Search Methods
13.6 Projection or Reduced Direction Search Methods for Bound-Constrained Problems
13.7 Successive Quadratic Programming or Projected Lagrangian Method
13.8 Narrow Valley Effect
13.9 The Nonlinear Constraints Effect
13.10 Tests
Chapter 14: Parametric Continuation in Optimization and Process Control
14.1 Introduction
14.2 Algebraic Constraints
References
Appendix A: Copyrights
Index
Related Titles
Buzzi-Ferraris, G., Manenti, F.
Fundamentals and Linear Algebra for the Chemical Engineer
Solving Numerical Problems
2010
Print ISBN: 978-3-527-32552-8
Buzzi-Ferraris, G., Manenti, F.
Interpolation and Regression Models for the Chemical Engineer
Solving Numerical Problems
2010
Print ISBN: 978-3-527-32652-5
Rice, R. G., Do, D. D.
Applied Mathematics And Modeling For Chemical Engineers
Second edition
2012
ISBN: 978-1-118-02472-0
Velten, K.
Mathematical Modeling and Simulation
Introduction for Scientists and Engineers
2009
Print ISBN: 978-3-527-40758-3, also available in electronic formats
Katoh, S., Yoshida, F.
Biochemical Engineering
A Textbook for Engineers, Chemists and Biologists
2009
ISBN: 978-3-527-32536-8
Sewell, G.
Computational Methods of Linear Algebra, Second Edition
Second Edition
2005
Print ISBN: 978-0-471-73579-3; also available in electronic formats
All books published by Wiley-VCH are carefully produced. Nevertheless, authors, editors, and publisher do not warrant the information contained in these books, including this book, to be free of errors. Readers are advised to keep in mind that statements, data, illustrations, procedural details or other items may inadvertently be inaccurate.
Library of Congress Card No.: applied for
British Library Cataloguing-in-Publication Data
A catalogue record for this book is available from the British Library.
Bibliographic information published by the Deutsche Nationalbibliothek
The Deutsche Nationalbibliothek lists this publication in the Deutsche Nationalbibliografie; detailed bibliographic data are available on the Internet at <http://dnb.d-nb.de>.
© 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Boschstr. 12, 69469 Weinheim, Germany
All rights reserved (including those of translation into other languages). No part of this book may be reproduced in any form – by photoprinting, microfilm, or any other means – nor transmitted or translated into a machine language without written permission from the publishers. Registered names, trademarks, etc. used in this book, even when not specifically marked as such, are not to be considered unprotected by law.
Print ISBN: 978-3-527-33274-8
ePDF ISBN: 978-3-527-66717-8
ePub ISBN: 978-3-527-66716-1
Mobi ISBN: 978-3-527-66715-4
oBook ISBN: 978-3-527-66714-7
Preface
This book is aimed at students and professionals needing to numerically solve scientific problems involving nonlinear systems and optimization.
We assume our readers have the basic familiarity with numerical methods, modeling, and optimization topics that any undergraduate student in scientific or engineering disciplines should have. We also recommend at least a basic knowledge of C++ programming.
Readers that do not have any of the above should first refer to the companion books in this series:
Guido Buzzi-Ferraris (1994), Scientific C++–Building Numerical Libraries, the Object-Oriented Way, 2nd ed., Addison-Wesley, Cambridge University Press, 479 pp, ISBN: 0-201-63192-X.Guido Buzzi-Ferraris and Flavio Manenti (2010), Fundamentals and Linear Algebra for the Chemical Engineer: Solving Numerical Problems, Wiley-VCH Verlag GmbH, Weinheim, 360 pp, ISBN: 978-3-527-32552-8These books explain and apply the fundamentals of numerical methods in C++.
Although many books on nonlinear systems and optimization approach these topics from a theoretical viewpoint only, we wanted to explain the theoretical aspects in an informal way, by offering an applied point of view of this scientific discipline. In fact, this volume focuses on the solution of concrete problems and includes many examples, applications, code samples, programming, and overall programs, to give readers not only the methodology to tackle their specific problems but also the structure to implement an appropriate program and ad hoc algorithms to solve it.
The book describes numerical methods, high-performance algorithms, specific devices, innovative techniques and strategies, all of which are implemented in a well-established numerical library: the BzzMath library, developed by Prof. Guido Buzzi-Ferraris at the Politecnico di Milano and downloadable from www.chem.polimi.it/homes/gbuzzi.
This gives readers the invaluable opportunity to use and implement their code with a numerical library that involves some of the most appealing algorithms in the solution of differential equations, algebraic systems, optimal problems, data regressions for linear and nonlinear cases, boundary value problems, linear programming, and so on.
Unfortunately, unlike many other books that only cover theory, all these numerical contents cannot be explained in a single book because of their application to concrete problems and the need for specific code examples. We, therefore, decided to split the numerical analysis topics into several distinct areas, each one covered by an ad hoc book by the same authors and adopting the same philosophy:
Vol. I: Buzzi-Ferraris and Manenti (2010), Fundamentals and Linear Algebra for the Chemical Engineer: Solving Numerical Problems, Wiley-VCH Verlag GmbH, Weinheim, Germany.Vol. II: Buzzi-Ferraris and Manenti (2010), Interpolation and Regression Models for the Chemical Engineer: Solving Numerical Problems, Wiley-VCH Verlag GmbH, Weinheim, Germany.Vol. III: Buzzi-Ferraris and Manenti (2014), Nonlinear Systems and Optimization for the Chemical Engineer: Solving Numerical Problems, Wiley-VCH Verlag GmbH, Weinheim, Germany.Vol. IV: Buzzi-Ferraris and Manenti (in press), Differential and Differential-Algebraic Systems for the Chemical Engineer: Solving Numerical Problems, Wiley-VCH Verlag GmbH, Weinheim, Germany.Vol. V: Buzzi-Ferraris and Manenti, The Attic Method in Mixed Integer and Linear Programming for the Chemical Engineer: Solving Numerical Problems, in progress.This book proposes algorithms and methods to solve nonlinear systems and optimization, whereas the companion books cover linear algebra and linear systems, data analysis and regressions, differential and differential-algebraic systems, and linear programming, respectively. After having introduced the theoretical content, all explain their application in detail and provide optimized C++ code samples to solve general problems. This allows readers to use the proposed programs to make tackling their specific numerical issues easier using BzzMath library.
The BzzMath library can be used in any scientific field in which there is a need to solve numerical problems. Its primary use is in engineering, but it can also be used in statistics, medicine, economics, physics, management, environmental sciences, biosciences, and so on.
Outline of This Book
This book deals with the solution of nonlinear systems and optimization problems in continuous space. Analogous to its companion books, it proposes a series of robust and high-performance algorithms implemented in the BzzMath library to tackle these multifaceted and notoriously difficult issues.
Chapter 1 discusses methods for handling the function root-finding problem. Algorithms are proposed in the renewed forms to exploit the multiprocessor machines. The use of parallel computing is also common to the chapters that follow.
One-dimensional optimization is broached in Chapter 2. Comparison methods, parabolic and cubic interpolations, gradient-based methods, and hybrid combinations are presented and implemented in chemical engineering examples.
Unconstrained optimization methods are discussed in Chapter 3. Heuristic methods, gradient methods and the conjugate direction methods are introduced together with Newton's method and modified Newton and quasi-Newton methods. Convergence and stop criteria are discussed, implemented in generalized classes, and used to optimize the design and operation of batch and fixed-bed reactors.
Chapter 4 has been devoted to large-scale unconstrained optimization problems, where problems related to the management of matrix sparsity and the ordering of rows and columns are broached. Hessian evaluation, Newton and inexact Newton methods are discussed.
To tackle tough problems, such as narrow valleys or discontinuities have to be solved, robust algorithms are presented in Chapter 5. A typical case relates to the estimation of kinetic parameters. One-, two-, and multidimensional optimization cases are handled separately and parallel computing plays a key role in improving both the robustness and the efficiency.
Algorithm robustness, using the openMP directives for shared memory parallel computing, is also investigated and implemented for function root-finding in Chapter 6.
Another important topic is the solution of nonlinear algebraic systems, which demands very robust algorithms. Chapter 7 illustrates the numerical methods for square systems in their sequential and parallel implementation. In addition to this, methods and techniques are proposed by separating the small and medium dimension problems, which are considered dense for large-scale systems, where the management of matrix sparsity is crucial. Many practical examples are provided.
Chapter 8 deals with underdimensioned nonlinear systems. It proposes a stable Gauss factorization for their solution and compares the novel method to the most common factorizations. Sparse underdimensioned nonlinear systems also have a dedicated class of algorithms.
Constrained optimization is broached starting from Chapter 9. The constraints are split into three categories: bounds, equality constraints, and inequality constraints. The relationship between primal and dual problems is discussed in further depth.
Chapter 10 introduces the Attic method for linear programming, and to which we dedicate a complete volume (Vol. 5 – Buzzi-Ferraris and Manenti, in progress) because of the huge significance of the method. The Attic method is compared to the existing Simplex and Interior Point methods, emphasizing its efficiency and robustness.
Quadratic programming is discussed in Chapter 11. The quadratic programming methods are implemented by handling the numerical evaluation of Hessian in a novel way. Reduced gradient and gradient projection are described as conventional methods; they are then compared to the novel proposed approach based on the Karush–Kuhn–Tucker direction projection method. It represents the basic core for the development of successive quadratic programming (SQP) methods.
Penalty function methods, barrier function methods, and hybrid penalty–barrier function methods are discussed in Chapter 12; they are implemented and adopted in several examples.
Chapter 13 illustrates the problem of constrained optimization by introducing the active set methods. Successive linear programming (SLP), projection, reduced direction search, SQP methods are described, implemented, and adopted to solve several practical examples of constrained linear/nonlinear optimization, including the solution of the Maratos effect.
Chapter 14 introduces parametric continuation in optimization and process control, which is of significant practical interest in chemical and process engineering. This topic is broached in Vol. 4 – Buzzi-Ferraris and Manenti (in press).
Notation
These books contain icons not only to highlight some important features and concepts but also to underscore that there is potential for serious errors in programming or in selecting the appropriate numerical methods.
New concepts or new ideas. As they may be difficult to understand, they require a change of point of view.
Description and remarks regarding important concepts; smart and interesting ideas.
Positive aspects, benefits, and advantages of algorithms, methods, and techniques in solving a specific problem.
Negative aspects and disadvantages of algorithms, methods, and techniques used to solve a specific problem.
Some aspects are intentionally neglected.
Caveat, risk of making sneaky mistakes, and spread errors.
Description of BzzMath library classes or functions.
Definitions and properties.
Conditioning status of the mathematical formulation.
Algorithm stability.
The algorithm efficiency assessment.
The problem, method, … is obsolete.
Example folders collected in WileyVol3.zip or in BzzMath7.zip files available at www.chem.polimi.it/homes/gbuzzi.
BzzMath Library Style
In order to facilitate both implementation and program reading, it was necessary to diversify the style of the identifiers.
C++ is a case-sensitive language and thus distinguishes between capital letters and small ones. Moreover, C++ identifiers are unlimited in the number of chars for their name, unlike FORTRAN77 identifiers. It is thus possible, and we feel indispensable, to use these prerogatives by giving every variable, object, constant, function, and so on, an identifier that allows us to immediately recognize what we are looking at.
Programmers typically use two different styles to characterize an identifier that consists of two words. One possibility is to separate the word by means of an underscore, that is, dynamic_viscosity. The other possibility is to begin the second word with a capital letter, that is, dynamicViscosity.
The style adopted in BzzMath library is described hereinafter:
Constants: The identifier should have more than two capital letters. If several words have to be used, they must be separated by an underscore.Another style-based decision was to standardize the bracket positions at the beginning and at the end of each block, to make C++ programs easier to read.
In this case also, programmers adopt two alternatives: some put the first bracket on the same row where the block starts, some others put it on the following line with the same indenting of the bracket that closes the block.
The former case results in the following style:
whereas the latter case results in
This latter alternative is adopted in the BzzMath library.
A third important style-based decision concerned the criterion to pass variables of a function either by value or by reference. In the BzzMath library, we adopt the following criterion:
If the variable is standard and the function keeps it unchanged, it is passed by value.If the variable is an object and the function keeps it unchanged, it is passed by reference and, if possible, as const type.If the variable (either standard or object) is to be modified by the function, its pointer must be provided.The object C only is modified in the following statements:Basic Requirements for Using the BzzMath Library
BzzMath library, release 7.0, is designed for a Microsoft Windows environment only. The library is released for the following compilers: Visual C++ 6 (1998), C++ 9 (2008), C++ 2010, C++ 2012, and INTEL 11.0.
openMP directives for parallel computing are available for the compilers C++ 2008, C++ 2010, C++ 2012, and INTEL 11.0.
All the compilers work on a 32-bit machine except for C++ 2010 and C++ 2012, which work on a 64-bit machine.
Moreover, FORTRAN users can either adopt all the classes belonging to the BzzMath library using opportune interfaces or directly use pieces of C++ codes in FORTRAN, by means of the so-called mixed language (see Appendix A of Vol. 2 – Buzzi-Ferraris and Manenti, 2010b).
The previous BzzMath library release, 6.0, can be installed on other operating systems, such as UNIX/LINUX with compatibility extended to gcc and icc (INTEL 11) LINUX compilers. The last version of release 6.0 is updated to the 20 May 2011 and will not undergo any further development. Moreover, the new release 7.0 has been extended quite significantly, particularly for classes dedicated to optimization and the exploitation of openMP directives.
How to Install Examples Collected in This Book
Download and unzip WileyVol3.zip from Buzzi-Ferraris's homepage (www.chem.polimi.it/homes/gbuzzi).
A Few Steps to Install Bzzmath Library
The following general tasks must be accomplished to use the BzzMath library on a computer.
Download BzzMath7.zip from Buzzi-Ferraris's homepage (www.chem.polimi.it/homes/gbuzzi).Unzip the file BzzMath7.zip in a convenient directory (e.g., in C:\NumericalLibraries\). This directory will be referred to as DIRECTORY in the following. This unzip creates the subdirectory BzzMath including other five subdirectories:Note that when a new directory is added to the PATH environment variable a; character must be included before specifying the new directory.
After having changed the PATH environment variable, you must restart the computer. At the next machine start, you can use BzzMathexe programs and/or the BzzMathTutorial.pps file, placed into the directory DIRECTORY\BzzMath\BzzMathTutorial.
Including BzzMath Library in a Calculation Program
Although the previous paragraph describes an operation that should need to be performed only once, the following operations are required whenever a new project is opened and at least one BzzMath library object is used:
The appropriate compiler library BzzMath.lib must be added to the project.The appropriate compiler must be defined by choosing one of the following alternatives:Moreover, when at least one BzzMath library object is used, it is always necessary to introduce the statementsat the beginning of the program, just below the BZZ_COMPILER selection. For example, using the INTEL 11 with openMP, the user must enter the following statements:
1
Function Root-Finding
Examples of this chapter can be found in the Vol3_Chapter1 directory in the WileyVol3.zip file available at www.chem.polimi.it/homes/gbuzzi.
1.1 Introduction
This chapter deals with the problem of finding the value that zeroes a function in the one-dimensional space:
(1.1)
We will, for instance, find the value of that zeroes the function
(1.2)
in the interval .
Only real problems involving the variable are considered.
Later, we describe the main iterative methods used to numerically solve this problem: these are referred to as iterative because they provide a series of values such that the value of for a sufficiently large approximates the solution with the required precision.
A series converges to the solution with the order if
(1.3)
where C and are appropriately selected.
The function could prove very complicated and it is not necessary to have it as an explicit algebraic function. For example, the problem might be to find either the value that zeroes the integral
(1.4)
with or the value for which a combination of the elements of the vector is zero during the integration of a differential system
(1.5)
When the function is particularly simple, however, certain special algorithms can be used.
Only general algorithms are proposed here and the potential peculiarities of the function are not exploited. This means, for instance, that special methods designed to find the roots of a -degree polynomial are not considered.
It is important to distinguish between the following three situations:
When zeroing a function, there is a considerable difference between the first two situations and the third one.
This chapter introduces the algorithms that solve the first two families of problems. The third family of problems, however, demands more robust algorithms (see Chapter 6).
When the interval of uncertainty for a function is known or the function is monotone, special algorithms can be used to guarantee the following appealing features:
It is also necessary to distinguish between two different situations in this discussion:
The computer is either a monoprocessor or parallel computing is unsuitable or cannot be exploited. Specifically, it is unsuitable when the function is very simple and it cannot be exploited when root-finding is part of a calculation in which parallel computations have already been used.It is possible – and helpful – to exploit parallel computing for function root-finding.The algorithms used to find roots of functions can be grouped into three families:
Analyzing simple problems and algorithms like this is particularly instructive as it illustrates the various ways of solving numerical problems: by using manual methods and classical analysis (without round-off errors) or, alternatively, a computer program (in which round-off errors play an important role) (Vol. 1 – Buzzi-Ferraris and Manenti, 2010a).
Anyone working on numerical problems may benefit from writing their own general root-finding function.
In the following, we will refer to the function evaluated in () simply as . We will also refer to the derivative as .
Before proceeding, however, we feel it is useful to reiterate an important statement from Vol. 1 (Buzzi-Ferraris and Manenti, 2010a).
When a problem is numerically solved on a computer, the round-off errors may become so relevant that even certain modifications to traditional classical analysis thought appear necessary.
Some methods for the collocation of a new point given a series of values demand function monotonicity.
A function, which is theoretically monotone in classical analysis (without round-off errors), can be constant in numerical analysis with respect to two points, if they are too close.
For instance, given the function
(1.6)
and two points and , the function value is the same in double precision because of the round-off errors.
The resolution, , represents the minimum distance required to ensure that a theoretically monotone function is also numerically monotone.
It is interesting to estimate the order of magnitude of for the root-finding problem. Let us suppose that the function can be expanded in Taylor series in the neighborhood of the solution, :
(1.7)
The term is not negligible with respect to when
(1.8)
In these relations, is the macheps of the calculation.
From (1.8), we obtain
(1.9)
The root of a function can be commonly found with a precision in the order of the macheps.
The previous discussion is valid only if ; otherwise, the function intersects the axis of abscissa horizontally and the points must be at least at a distance
(1.10)
When at the solution, the precision of a numerical method is proportional to rather than to . In double precision on 32 bit computers, the error is on the 6–7th significant digit rather than the 14–15th.
1.2 Substitution Algorithms
These algorithms transform the problem to produce expressions such as
(1.11)
so as to exploit the iterative formulation:
(1.12)
For example, given the function
(1.13)
we have
(1.14)
(1.15)
(1.16)
(1.17)
Starting from an initial point , the following series are obtained in correspondence with the functions (1.14)–(1.17), respectively:
: the series diverges.: the series diverges.: the series diverges.: the series properly converges to the problem solution.This approach was often adopted when the calculations were manual as no general programs were available to solve the problem in its original form.
However, this kind of procedure also has numerous disadvantages.
The iterative procedure can also diverge with very simple functions; for example, it can diverge with a linear equation.
If the iterative formula is linear
(1.18)
the method converges only if
(1.19)
It is useful to plot the two trends to illustrate the reason for this limitation:
(1.20)
(1.21)
and proceed with the iterations. Note that the method diverges if equation (1.21) has a slope and (Figure 1.1).
Figure 1.1 Convergence conditions of substitution algorithms.
When the procedure converges, the convergence is usually slow.
For example, the problem mentioned above
(1.22)
is solved using Newton's method, which has quadratic convergence, yielding the following series results: .
If the substitution method converges, the convergence speed can be improved using the Aitken method (Buzzi-Ferraris and Manenti, 2010a).
Given a linearly convergent series of , the Aitken formula generates a new and more efficient series :
(1.23)
For example, given the series
(1.24)
linearly converging to the value , the series from Table 1.1 is obtained.
A laborious preliminary study is required to convert the problem into a convergent form.
The substitution method is only useful from an educational point of view.
When such calculations were performed manually, a lot of time had to be devoted to search for a special method for solving a specific problem. We encountered several disadvantages to retaining this approach when developing computer programs:
Table 1.1 Series for the iterative method.
2.000 0002.050 0002.074 2442.097 0642.085 4482.095 0742.090 5022.094 6572.094 5462.092 7572.094 5732.094 5522.094 2002.094 5552.094 3962.094 5532.094 4832.094 5522.094 5212.094 5382.094 5462.094 5492.094 5512.094 5521.3 Bolzano's Algorithm
The following conditions must be satisfied to use Bolzano's method.
If the previous conditions are both satisfied, Bolzano's theorem has demonstrated the existence of at least one solution within the interval .
From a numerical point of view, if a single processor is available, the Bolzano method is equivalent to the bisection method.
The algorithm is very simple: the function is evaluated at the central point:
(1.25)
if has the same sign of , replaces the original boundary ; otherwise replaces . This procedure allows the original interval to be iteratively halved, preserving the two necessary conditions as above.
Let us denote the width of the initial interval . As the interval is halved at each iteration, after iterations, we have
(1.26)
From (1.26), it is possible either to evaluate a priori the interval of uncertainty after a predetermined amount of iterations or, given the final interval, to know a priori the number of iterations required to reach it.
Let us select a point different from the center of the interval. If we are lucky, the solution is in the smallest of the two new subintervals; otherwise, the solution is in the largest interval. Consequently, the maximum interval of uncertainty is larger than Bolzano's interval of uncertainty at each iteration.
In the case of a single processor, Bolzano's method minimizes the maximum interval of uncertainty.
If several processors are available simultaneously, Bolzano's method can be generalized: the function should be simultaneously evaluated in an evenly spaced number of points equal to the number of available processors and within the interval of uncertainty.
Equally, in the general form available for several processors, Bolzano's method retains the advantages of the corresponding method for single-processor machines.
However, it also has one major deficiency:
Bolzano's algorithms do not exploit the peculiarities of the function. While this is an advantage for complicated functions (e.g., nonmonotone or multimodal), it may prove a handicap in other situations where the function is smooth.
For example, if we apply Bolzano's method to the linear function
(1.27)
or to any other function, the same number of iterations is required.
1.4 Function Approximation
The key point on which this family of algorithms is based is the approximation of the function through a simpler formulation that makes searching for the root easier; the function root is then adopted as the initial guess for the successive iteration.
This kind of algorithm can be implemented using two different strategies.
The first strategy is the only one of any use when the interval of uncertainty is unknown. If the function is monotone, it is always possible to find an interval of uncertainty and, consequently, to use ad hoc techniques to find its root.
Algorithms for function root-finding with a known interval of uncertainty are discussed below, whereas the other situation is considered in Section 1.6.
1.4.1 Newton's Method
If the function and its derivatives are continuous, the same function can be expanded in Taylor series in the neighborhood of :
(1.28)
If we stop the Taylor series expansion at the first-degree term, we obtain the equation of the tangent to the curve in :
(1.29)
The point , where the tangent intersects the axis of abscissas, can be obtained by imposing the condition
(1.30)
thus
(1.31)
Given an initial guess and evaluating the function and its first derivative in correspondence with this point, it is possible to get a value for the new point using (1.31).
For example, to evaluate the square root of 2 using Newton's method, the root of the function
(1.32)
needs to be found. Relation (1.31) becomes
(1.33)
By choosing as the initial guess, the series , , , and is obtained, and the final value is already a good approximation of the solution.
If the point is close to the solution and the function can be expanded using a Taylor series
(1.34)
Dividing by
(1.35)
If the method converges
(1.36)
Newton's method has quadratic convergence.
This method too, however, has certain disadvantages.
Consider, for example, the function
(1.37)
which is equal to zero in . Starting from , the series is obtained and the method diverges.
The divergence of Newton's method with a monotone function is reported in Figure 1.2.
The advantage of Newton's method is in the fast convergence of the cases in which it effectively converges. Convergence is quadratic in these cases.
Newton's method is never used in the BzzMath library classes dedicated to root-finding.
1.4.2 The Secant Method
In relation (1.31), by approximating the derivative by the line joining the points and , the following relation takes place:
(1.38)
Two points and are needed to initialize the algorithm.
Let us consider again the previous example of evaluating the square root of 2, which corresponds to zero in equation (1.32).
Let and be the initial guesses. The series is , , , , and .
The secant method has a convergence speed raised to the power of 1.618. It is slower than Newton's method, but the first derivative does not need to be evaluated. When the computational effort involved in evaluating the derivative is in the order of the computational time required to calculate the function, the secant method performs better numerically than Newton's method.
Figure 1.2 Divergence of Newton's method with monotone function.
The secant method has the same pros and cons as Newton's method, except for the need to provide the analytical expression of the first derivative of the function.
The secant method in its pure iterative form is never used in the classes of the BzzMath library dedicated to root-finding.
1.4.3 Regula Falsi Method
The regula falsi algorithm is very similar to the previous one. The difference is in the support points adopted to linearize the function: the last two values at each iteration are used in the secant method, whereas the boundaries of the interval of uncertainty are adopted in the regula falsi method.
The advantage this has over the previous criteria is that the new iteration falls within the interval of uncertainty and thus method convergence is ensured.
The new prevision replaces the bound where the function assumes the same sign.
Its main disadvantage is that its convergence speed is slower than the secant method.
To calculate the square root of 2 using this method
(1.39)
and with and as initial guesses, we obtain the series , , , , , , , , and .
This algorithm too is important only from an educational point of view, as it is the simplest algorithm in this family, ensuring convergence when the interval of uncertainty is known.
The efficiency of the regula falsi method can be improved by using the two best values from the previous iterations rather than the boundary values for the linear approximation. This device is implemented by checking the interval of uncertainty.
This device has the advantage of a better convergence speed, similar to the one of secant method. Many programs adopt this strategy as a basic algorithm (often combined with Bolzano's method to guarantee convergence even with complex and/or nonmonotone functions).
Other authors attempt to improve the efficiency of the regula falsi method with another device. Typically the method would not work satisfactorily when the function has very different absolute values at the boundaries of the interval. If we select these boundaries as support points for linearization, we can guarantee convergence to the solution, but the interval reduction might be very small. The device consists of the selection of a smaller value (i.e., dividing by 2) in correspondence with the boundary where the function has its maximum absolute value. A reduction of this kind in the ordinate can also be proposed in successive iterations.
The regula falsi method using the two best values of the previous iterations rather than the boundary values for the linear approximation is implemented in the classes of the BzzMath library dedicated to root-finding.
1.4.4 Muller's Method or Parabolic Interpolation
This method approximates the function with a parabola. Contrary to the previous methods, three points are needed to estimate the parameters of the parabola.
The algorithm has the following pros and cons.
As the method approximates the function using a parabola, it is necessary to choose a solution from the two possible predictions.
Except in the case of polynomial root-finding and complex roots, there are more efficient and easier-to-implement algorithms.
The parabolic interpolation method is never used in the BzzMath library classes dedicated to root-finding.
1.4.5 Hyperbolic Interpolation Method
The function can be approximated using the simplest rational function: a hyperbole:
(1.40)
Given three support points, the three parameters , , and can be calculated by solving the linear system:
(1.41)
Unlike Muller's method, this algorithm has a single solution:
(1.42)
Moreover, an interpolation with a rational function is usually preferable to using a polynomial function (Vol. 2 – Buzzi-Ferraris and Manenti, 2010b).
The hyperbolic interpolation method is used when the root-finding classes from the BzzMath library have performed three iterations.
1.4.6 Inverse Polynomial Interpolation Method
Muller's method has to select the new point between the two solutions of the interpolating parabola. This difficulty is enhanced by interpolating the function with a higher degree polynomial.
However, it is possible to overcome it by adopting an inverse interpolation. In other words, a polynomial in the dependent variable can be developed.
Given support points, it is possible to calculate the polynomial parameters, approximating with respect to :
(1.43)
and estimate the polynomial prediction in .
Since the polynomial is useful to predict the point where , the Neville method, which does not require any parameter estimation, can be adopted (see Buzzi-Ferraris and Manenti, 2010b).
The inverse polynomial interpolation method is never used in the BzzMath library classes dedicated to root-finding.
1.4.7 Inverse Rational Interpolation Method
An exact polynomial interpolation can deteriorate as the polynomial degree increases (Buzzi-Ferraris and Manenti, 2010b), whereas this problem does not arise with rational functions.
An inverse interpolation can be effectively exploited using a rational function rather than a polynomial and the Bulirsch–Stoer algorithm rather than the Neville method.
It is worth remarking that this method cannot be implemented on its own in a program, rather certain additional (and less efficient) methods must be integrated with it to ensure convergence. These additional algorithms should be used whenever the inverse rational function interpolation encounters any difficulty.
Inverse interpolation with rational functions is very efficient and is the ideal basis for the development of a general function root-finding program, even though it is slightly more complex than the other algorithms.
In the BzzInterpolation class, the Neville function allows an interpolation with polynomials to be performed, whereas BulirschStoer performs an interpolation with rational functions.
These functions also provide an estimation of the improvement obtained by increasing the order of the interpolating polynomial.
The results are
The solution is 3.999963311389726e -002. The error in the prediction obtained using all the five points is 0.16% with Neville and 0.000013% with Bulirsch–Stoer.
The inverse rational interpolation method is used as the basic method when more than three iterations of the BzzMath classes for function root-finding have been performed.
1.5 Use of a Multiprocessor Machine with a Known Interval of Uncertainty
If an interval of uncertainty is known and some processors are available to find the root, Bolzano's method can effectively be coupled with a function approximation method.
There are various strategies when more than a single processor, , are available. The following is the one adopted in the BzzMath library.
Suppose we have available processors. The first one is used to calculate the function in the point selected by the efficient method being used, whereas the remaining processors evaluate the function in the points opportunely placed in the two new subintervals, created by the splitting of the interval of uncertainty.
1.6 Search for an Interval of Uncertainty
The search for an interval of uncertainty is very simple if it is known a priori that the function is monotone. Specific techniques for the most general case are explained in Chapter 6.
Heuristic techniques or function approximations (i.e., the secant method) can be exploited when the function is monotone.
One widely adopted technique is to move in the steepest descent direction with a step , which progressively increases (e.g., with ), until an interval of uncertainty is found.
If necessary, the prediction provided by the secant method can also be exploited on the condition that said prediction is checked: it should be in the order of the current value of .
When the interval of uncertainty is identified, a dedicated algorithm capable of exploiting it can be adopted.
Suppose we have available processors: it is useful to place points at a distance that gradually increases up to the detection of a proper interval of uncertainty.
It is worth remarking that, if the function is monotone, it is possible to find a rough interval of uncertainty and, once detected, the solution can be obtained using very efficient methods, although the interval of uncertainty is wide.
1.7 Stop Criteria
In all problems solved using iterative methods, it is important to define a criterion to stop the evaluations either when a reasonable approximation of the solution is achieved or when convergence is impossible.
There are various possibilities where the specific problem of function root-finding is concerned and it is advisable for us to analyze some of them here.
(1.45)
(1.46)
(1.47)
(1.48)
(1.49)
(1.50)
(1.51)
(1.52)
(1.53)
(1.54)
(1.55)
(1.56)
(1.57)
(1.58)
1.8 Classes for Function Root-Finding
In the BzzMath library, the following classes are implemented for function root-finding:
The BzzFunctionRoot class is designed to search for the root of a function when an interval of uncertainty is given and when the user wants to find one solution within that interval.
The function must be continuous within the interval. The BzzFunctionRoot class will not benefit from several processors.
Conversely, if multiprocessor computers and openMP directives are used, the BzzFunctionRootMP class should be adopted under the same above-mentioned condition. It is worthwhile remarking that if we use objects from the BzzFunctionRootMP class, we must also use a compiler that supports openMP directives and activate them.
The BzzFunctionRootRobust class is designed to search for root functions in the following cases:
When the interval of uncertainty is unknown and the function is not monotone. When the function has many roots and all of them need to be found.When the function is not evaluable in certain regions of the domain.When the function does not change sign in the neighborhood of certain roots.Only the BzzFunctionRoot and BzzFunctionRootMP classes are described in this chapter, whereas the BzzFunctionRootRobust class is explained in Chapter 6.
The BzzFunctionRoot class has several constructors. The simplest ones are
The first is the default one and no data are needed. The second requires the boundaries and of the interval of uncertainty.
If the values of the function at the two boundaries have the same sign, the calculation is stopped immediately and a warning message is printed out.
The third constructor requires the function values at the interval boundaries, and .
As the object is defined, it is possible to use the operator () to initialize it either with different values for the interval boundaries and/or with a different function.
The search for the solution is performed by the overloaded operator (). Two different versions of the overloaded operator () are implemented:
For example
Using this operator, the search proceeds until the required accuracy is obtained. On the other hand, using the operator
the root-finding process is stopped when the number of calculations is equal to niter, which is an int that limits the number of calculations. In this instance, it is worth stressing the following feature.
If the function is iteratively called, the object accounts for the previous calculations and executes the calculations as if an overall single call were carried out.
In other words, invoking or six times is the same from a numerical point of view. This property allows the user to manage the object and check the root-finding of the same while the search is proceeding without compromising efficiency. Moreover, it should be stressed that the object manages its own data, which are masked to the user inside the operator () and therefore there is no risk of modifying them unintentionally, for example, by using several objects simultaneously. As the object is aware of previous calculations, it does not execute any new calculations when it has solved the problem with the required accuracy even though it is requested to.
It is possible to adjust the accuracy of the root-finding process using the following functions: SetTolRel, SetTolAbs, and SetTolY. The former two functions allow the default value to be modified in (1.58). The default values are and , respectively. The function SetTolY allows the value for test (1.56) to be modified; the default value is .
For example, the tolerances of the object z are changed in the following code sample:
If the function is iteratively called, the tolerances can be modified anywhere. For example, accuracy can be progressively increased. The object preserves its history in this case too and exploits it as much as possible also.
In the operator (), the return is a char and indicates the reason for which the search was stopped. It can assume the following values:
-2: the values of the function at the two boundaries have the same sign.-1: the program is stopped by the user through the global variable bzzStop = 1.0: the program encountered certain numerical problems or niter is less than 1.1: the number of iterations is equal to niter in this specific call of the function.2: the test on the interval of uncertainty (1.58) is satisfactory:3: the absolute value of the function is smaller than yTolAbs.Moreover, it is possible to use the BzzPrint and BzzMessage functions, for each class implemented (Buzzi-Ferraris and Manenti, 2010a), as they summarize all the above information. For example
The function root-finding is performed by exploiting one of the following algorithms depending on the situation: the inverse rational interpolation with the Bulirsch–Stoer prevision; the regula falsi method using the two best values of the previous iterations; the hyperbolic method with three points; Bolzano's method. The former is the basic method. The regula falsi method is adopted either when only two points are available or when there are problems with the Bulirsch–Stoer prediction, that is, when the error between the predictions does not decrease with the degree of rational functions. The hyperbolic method is adopted when only three points are available. Bolzano's methods are adopted when the first three methods do not improve the objective function significantly; it is also used when the interval of uncertainty does not decrease satisfactorily in the iterations of the three previous methods. A check is performed for each new prediction to ensure that it falls within the search interval and is also significantly distant from the points already analyzed.
The BzzFunctionRootMP class has the same constructors and functions as the BzzFunctionRoot class. The difference between them lies in the fact that the former requires more processors to be available in a shared memory architecture. Consequently, anyone availing of it must use a compiler that accepts openMP libraries and the right version of the BzzMath library. Moreover, the selected compiler must be pointed out by means of an ad hocdefine command (see the readme.pdf file in the tutorial).
If we are using the INTEL C++ 11 compiler for an MS Windows environment, for example, we need to start the program as follows:
1.9 Case Studies
This section illustrates a set of case studies in which root-finding plays an important role in chemical engineering including the calculation of the volume of a nonideal gas, bubble point, and zero-crossing. However, these scenarios also crop up in several other areas. For instance, the calculation of the volume of a nonideal gas is a typical problem in fluid dynamics, whereas the zero-crossing problem is very common in all disciplines involving differential and differential-algebraic systems as convolutions models, such as the optimal control for electrical and electronic purposes.
1.9.1 Calculation of the Volume of a Nonideal Gas
The calculation of the volume of a nonideal gas is a root-finding problem. Let us consider, for example, Rice's problem (Rice, 1993) where the Beattie–Bridgeman equation of state joins the temperature , the pressure , and the molar volume through the formula
(1.60)
where is the gas constant and
(1.61)
(1.62)
(1.63)
(1.64)
1.9.2 Calculation of the Bubble Point of Vapor–Liquid Equilibrium
The calculation of vapor–liquid equilibrium is a classic problem usually included in the introductions to books on numerical analysis. It is, in fact, conceptually and mathematically very simple, and since it does not require any deep knowledge in the chemical engineering, it is also very useful to a broader readership especially as a dedicated code has not been used for this specific problem.
Chemical and industrial engineering both provide a series of examples involving nonideal vapor–liquid equilibrium calculation. For example, the equations governing a flash drum separator with molar fractions , ,…, , flow rate , and enthalpy are
component mass balance: overall mass balance: equilibrium: stoichiometry: enthalpy balance:where , , and are the vapor, liquid, and inlet flow rates; , , and are their corresponding enthalpies; and are the molar fractions; is the duty provided to the flash; and are the equilibrium constants. In the nonideal case, . The system consists of equations, where is the number of components. If we assign , , , , and and the functions , , and are known, the system is solved in its equations to get the variables , , , , and .
Nevertheless, there are several variants of the problem that may be of interest. It is possible to assign the value for some of the previous variables and consider an equal number of the remaining ones as unknown.
Some combinations of given-unknown variables are not valid since they lead to ill-posed physical problems and therefore to unfeasible or unreliable numerical solutions.
Two important problems are the bubble and dew point calculations, obtained by zeroing and , respectively, while the duty is unknown.
