Root Bracketing For Two Nonlinear Equations
By
Namir C Shammas
Finding roots of nonlinear equations has attracted a lot of attention among mathematicians and numerical analysts. Today’s scientific and graphics calculators have built in solvers that use root-seeking algorithms to solve user-defined non-linear equations. When it comes to solving single nonlinear equations you get to choose from a myriad of algorithms. They range from root-bracketing methods like the Bisection and the False Position algorithms, to the jaguar of root-seeking algorithms, namely the Newton-Raphson method. In as much as a programmer feels spoiled, like a kid in a candy store, by so many methods that solve the roots of single nonlinear equations, the reverse is true when it comes to solving multiple nonlinear equations. In this case, it’s slim picking! Looking at numerical analysis textbooks, a programmer has no route except for using Newton’s method. And if that was not bad enough, numerical analysis books comment that Newton’s method may not work as well for multiple nonlinear equations as it did for their single counterparts! Why? Because the derivatives of the nonlinear equations may lead the iteration process through a wild goose chase! This situation is a nightmare akin to going to a restaurant and being given a single-menu item. To add insult to injury, the restaurant folks inform you that the eatery may even be out of that single menu item!
Being a numerical analysis enthusiast (or fanatic if you listen to my brother) I often wondered about the availability of methods like Bisection and False Position for multiple nonlinear equations. The number-crunching gurus have declared in their books that such methods are difficult to implement for multiple nonlinear equations. “Can’t we at least try”, I often asked myself? Searching through the Internet for such methods led me nowhere. I decided to give it a stab myself. The results are versions of the Bisection and False Position algorithms that work with two nonlinear equations. This article presents these algorithms. As a bonus, I will also include an extended version of the Newton method which works for multiple nonlinear equations. As far as this article is concerned, I will focus on examples that solve sets of two nonlinear equations. So, fasten you seat built, cause her we go!
The Bisection method represents the simplest and albeit the slowest method for finding a root of a nonlinear function. The method requires you to define an interval of X, call it Xa and Xb, where F(Xroot) = 0 and F(Xa) and F(Xb) are non-zero values that have opposite signs. The method systematically selects the value in the middle of that range to close in on the root. Each iteration replaces either Xa or Xb with the median value and eventually narrows the range for the root within the pre-assigned tolerance value. While the Bisection method is no speed demon, it always finds the root. The Bisection method does not use the derivatives of the function F(x), as does Newton’s method, and is therefore immune to the bad side effects caused by derivative values approaching zero. Such low value for the derivates can lead the iteration on a wild goose chase.
The False Position method is similar to the Bisection, except it uses the values of the functions at F(Xa) and F(Xb) as weights to decide how to locate the next guess. Unlike the Bisection, the False Position does not narrow the initial interval for the root. Instead, the method converges when the new guess for the root and the value of Xa or Xb it replaces are within the tolerance limit. This is an important distinction to keep in mind when coding the False Position.
The task of defining a Bisection method for two nonlinear equations took me through different approaches. In the first approach, I sought to define two ranges—of for the values of X and the other for the values of Y to solve F(X, Y) = 0 and G(X, Y) = 0. The main question was how to select the ranges for X and Y regarding the values of functions F and G? I though of a range defined by (Xa,Ya) and (Xb,Yb) such that F(Xa,Ya) and F(Xb,Yb) have opposite signs and, likewise, G(Xa,Ya) and G(Xb,Yb) have opposite signs. The next question was how do the signs of F(Xa,Ya) and G(Xa,Ya) relate to each other? Do they have to maintain the same signs or opposite signs? The same question goes for the values of F(Xb,Yb) and G(Xb,Yb). After testing a few methods I began to see why the numerical analysis text books have written off the Bisection method for multiple nonlinear equations.
After several years, I eventually gave up on the above approach and took a new and simpler one. I focused on the F(X,Y) and G(X,Y) are both zero (the Z vertical plane). Each of these functions yields a line, a curve, or even a spherical shape. In other words, when the functions F and G are both zero, there are two functions, call them Yfx and Ygx that relate the values of X and Y. Where these two functions intersect you will find the root or roots of F(X, Y) and G(X, Y). The function Yfx(X) calculates the value of Y for a given X when F(X, Y) is zero. Likewise, the function Ygx(X) calculates the value of Y for a given X when G(X, Y) is zero. The definition of these two spin-off functions allows me to define a Bisection method for two nonlinear equations that works as follows:
Figure 1 shows the curves for the functions Yfx(X) and Ygx(X) and the root-containing interval (Xa, Xb).
Figure 1. The intersection of the two nonlinear functions.
Once I obtained a Bisection method for two nonlinear equations, defining a False Position version was easy. The False position calculates the new guess as Xm = (DeltaB * Xa – DeltaA * Xb) / (DeltaA – DeltaB). The methods compares the value of Xm with the value of Xa or Xb (whichever Xm replaces) to determine if an acceptable root value was reached.
Figure 2 shows the three equations involved in refining the guesses for the roots of functions F(X, Y) and G(X, Y). In Figure 2, I use F, G, Fx, Fy, Gx, and Gy to abbreviate F(X, Y), G(X, Y), dF(X, Y)/dX, dF(X, Y)/dY, dG(X, Y)/dX, and dG(X, Y)/dY, respectively. Equations 2 and 3 in Figure 2 show the guess refinements for X and Y, respectively. Newton’s method performs well if the initial guesses are close to the actual roots. Depending on the initial guess values and the function derivatives, Newton’s method may converge on different roots.
Figure 2. The three equations used in refining the guesses for the roots of functions F(X, Y) and G(X, Y).
As a bonus in this article, I included a method that extends Newton’s method in solving multiple nonlinear equations. The method is based on my DDJ article in June 2002 (see the Algorithm Alley column), but is a simpler version. I replace the derivative dF(X, Y)/dX (abbreviated by F’(X, Y)) with the expression F’(X, Y) – (F’’(X, Y) * F(X, Y)) / (2 * F’(X, Y)), where F’’(X, Y) is the second derivative of F(X, Y) with respect to X. I used the similar expressions to extend the derivative of F(X, Y) with respect to Y, and the derivatives of G(X, Y) with respect to X and Y. By adding an extra term after the value of the first derivative, I make Newton’s method accelerate in converging to the roots. This method works for any number of nonlinear equations.
Let’s shift our focus to some C++ code that implements the Bisection, False Position, Newton, and Extended Newton methods. Listing 1 shows the source code for the Root.cpp file that contains the following groups of code:
The source code defines three sets of functions. Each set has four C++ functions that represent F(X, Y), G(X, Y), Yfx(X), and Ygx(X). The functions YFx(X) and Ygx(X) relate the values of X to Y when F(X, Y) and G(X, Y) are zero, respectively. Each function set has an integer appended to the names of the C++ functions. The first set of C++ functions model equations that have quadratic terms for both X and Y. The second set is an example of equations that contain exponential functions. The third set is an example of nonlinear functions that require iterative solution (at least function F(X, Y) does in this article). The code for function fx3 uses a simple implementation of Newton’s method to speed up convergence to calculate Y for a give value of X. I implemented function fx3 such that it returns the current iteration value, if the number of iterations exceeds the maximum limit.
The function fDeriv1X and fDeriv2X calculate the first and second derivates of a function with respect to X, respectively. Likewise, the functions fDeriv1Y and fDeriv2Y calculate the first and second derivates of a function with respect to X, respectively. All these C++ functions work with functions that take arguments for X and Y.
The class CRoot models a solver that implements the four methods we are dealing with. The class declares the following protected data members:
The class declares a constructor that has default parameters fTolerX and fTolerY that assign tolerances for the X and Y root values, respectively. The constructor also performs initialization for the data members of the class. The class also declares three access member functions. The member function getIters returns the number of iterations. The member functions getX and getY recall the values of the roots for X and Y, stored in data members m_fx and m_fy, respectively.
The member function Bisection implements the Bisection method that I described earlier in this article. The member function has the following parameters:
The member function stores the approximated roots in data members m_fx and m_fy. The member function returns a Boolean value to indicate success or failure. The member function monitors the number of iteration just in case the arguments for parameters fx and fy represent functions that do not meet the basic assumptions.
The member function FalsePosition implements the False Position method that I mentioned earlier article. The member function has the same parameter list as the member function Bisection.
The member function BisectionIterative implements a special version of the Bisection method where either or both functions Yfx(X) and Ygx(X) require iterations to calculate their results. The member function has the following parameters:
I did not implement a version of the FalsePosition member function that handles iterations in calculating Yfx(X) and Ygx(X). I leave that as an exercise for you.
The member function Newton implements the Newton’s method using the equations in Figure 2. The member function has the following parameters:
The member function ExNewton implements the extended version of Newton’s method as described earlier in this article. The member function has the same parameter list as member function Newton.
All the root-seeking member functions initialize the data members m_fx and m_fy with very high values. If the iteration fails, each member function returns false without updating the values in data members m_fx and m_fy. If the iterations succeed, each member function assigns the root approximations to the data members m_fx and m_fy and then returns true.
The function main tests the class CRoot with the three sets of C++ functions. The function main performs the following general tasks:
In the case of the first two tested mathematical functions, the function main performs the following tasks:
In the case of the third function set, the function main performs the following tasks:
Figure 3 shows the contents of the output file RootBrak.txt which mirrors the console output.
In solving the first set of mathematical functions, the function main uses the range (0, 3) for X with the Bisection and False Position methods. The Bisection method takes 15 iterations compared to just 8 for the False Position. In testing the Newton and extended Newton methods, the function main uses the initial guesses of 1.5 for each X and Y. Both methods return results that are better than their root-bracketing counterparts and in just 5 iterations! In this example the extended Newton method did not show improvement over the standard Newton method.
In solving the second set of mathematical functions, the function main uses the range (0, 3) for X with the Bisection and False Position methods. The Bisection method takes 15 iterations compared to 13 for the False Position. In testing the Newton and extended Newton methods, the function main uses the initial guesses of 1.5 for each X and Y. Surprisingly, Newton’s methods takes 17 as opposed to just 5 iterations for the extended Newton method! Here, the latter method shows superiority.
In solving the third set of mathematical functions, the function main uses the range (1, 5) for X with the Bisection method. The Bisection method takes 20 iterations. In testing the Newton and extended Newton methods, the function main uses the initial guesses of 2.5 for each X and Y. Newton’s methods takes 6 iterations as opposed to just 4 iterations for the extended Newton method! Here, again the latter method shows superiority.
This article presented versions of the Bisection and False Position algorithms that work with pairs of nonlinear equations. The article also presented an extended Newton method that works with multiple nonlinear equations. The Bisection and False Position methods, while slower than the Newton method and extended Newton method, are more certain to find a root if at least one exists in the specified range of values. I hope that my efforts push algorithms for solving nonlinear equations just a bit beyond current limits. I also hope that they can inspire other numerical analysis enthusiast into implementing versions of Bisection and False Position for more than two nonlinear functions.
Figure 3. The output.
Solving f(x,y) = x^2 + y^2 - (2.1^2 + 3.1^2)
g(x,y) = x^2 - y^2 - (2.1^2 - 3.1^2)
Exact roots are x = 2.1 and y = 3.1
Using Bisection
Initial range for x is (0, 3)
Number of iterations 15
Roots are at x = 2.09995 and y = 3.09999
Using False Position
Initial range for x is (0, 3)
Number of iterations 8
Roots are at x = 2.09999 and y = 3.1
Using Newton's Method
Initial guesses for x = 1.5 and for y = 1.5
Number of iterations 5
Roots are at x = 3.1 and y = 2.1
Using Extended Newton's Method
Initial guesses for x = 1.5 and for y = 1.5
Number of iterations 5
Roots are at x = 3.1 and y = 2.1
Solving f(x,y) = e^x - y^2 - (e^2 - 4^2)
g(x,y) = x^3 - e^y - (2^3 - e^4)
Exact roots are x = 2 and y = 4
Using Bisection
Initial range for x is (0, 3)
Number of iterations 15
Roots are at x = 1.99997 and y = 4.00001
Using False Position
Initial range for x is (0, 3)
Number of iterations 14
Roots are at x = 1.9999 and y = 3.99991
Using Newton's Method
Initial guesses for x = 1.5 and for y = 1.5
Number of iterations 17
Roots are at x = -3.59715 and y = -2.93911
Using Extended Newton's Method
Initial guesses for x = 1.5 and for y = 1.5
Number of iterations 5
Roots are at x = 2 and y = 4
Solving f(x,y) = x^2 - y^2 * e^x -(3^2 - 5^2 * e^3)
g(x,y) = e^x - y - (e^3 - 5)
Exact roots are x = 3 and y = 5
Using Bisection
Initial range for x is (1, 5)
Number of iterations 20
Roots are at x = 3 and y = 5
Using Newton's Method
Initial guesses for x = 2.5 and for y = 2.5
Number of iterations 6
Roots are at x = 3 and y = 5
Using Extended Newton's Method
Initial guesses for x = 2.5 and for y = 2.5
Number of iterations 4
Roots are at x = 3 and y = 5
Here is the listing for the source code. The C++ program was tested using Borland Builder C++ 6.
/* Program that calculates the roots of two nonlinear equations using special versions of the bisection and false position methods Author: Namir Clement Shammas. Date: September 5, 2002 */ #pragma hdrstop #include#include #include void pressAnyKey(char* sMsg) { char c[2]; cout << sMsg; cin.getline(c, 1); } //-------------- FIRST SET OF FUNCTIONS ----------------------- double f1(double x, double y) { // f(x,y) = x^2 - y^2 - (2.1^2 - 3.1^2) return x * x - y * y - (3.1 * 3.1 - 2.1 * 2.1); } double g1(double x, double y) { // g(x,y) = x^2 + y^2 - (2.1^2 - 3.1^2) return x * x + y * y - (3.1 * 3.1 + 2.1 * 2.1); } double fx1(double x) { // f(x,y) = x^2 - y^2 - (2.1^2 - 3.1^2) return sqrt(3.1 * 3.1 - 2.1 * 2.1 + x * x); } double gx1(double x) { // g(x,y) = x^2 + y^2 - (2.1^2 - 3.1^2) return sqrt(3.1 * 3.1 + 2.1 * 2.1 - x * x); } //-------------- SECOND SET OF FUNCTIONS ----------------------- double f2(double x, double y) { // f(x,y) = e^x - y^2 - (e^2 - 4^2) return exp(x) - y * y - (exp(2) - 16); } double g2(double x, double y) { // g(x,y) = x^3 - e^y - (2^3 - e^4) return x * x * x - exp(y) - (8 - exp(4)); } double fx2(double x) { // f(x,y) = e^x - y^2 - (e^2 - 4^2) return sqrt(exp(x) - (exp(2) - 16)); } double gx2(double x) { // g(x,y) = x^3 - e^y - (2^3 - e^4) return log(x * x * x - (8 - exp(4))); } //-------------- THIRD SET OF FUNCTIONS ----------------------- double f3(double x, double y) { return x * x - y * y * exp(x) - (9 - 25 * exp(3)); } double g3(double x, double y) { // g(x,y) = e^x - y - (e^3 - 5) return exp(x) - y - (exp(3) - 5); } double fx3(double x, double y, int nMaxIter) { // f(x,y) = x^2 - y^2 * e^x -(3^2 - 5^2 * e^3) // this function is iterative double h, fDiff, fTolerance = 1.0e-5; int nIter = 0; // do { h = fabs(y) > 1 ? 0.01 * y : 0.01; fDiff = 2 * h * f3(x, y) / (f3(x, y+h) - f3(x, y-h)); y -= fDiff; if (++nIter > nMaxIter) return y; } while (fabs(fDiff) > fTolerance); return y; } double gx3(double x, double, int) { // g(x,y) = e^x - y - (e^3 - 5) return exp(x) - (exp(3) - 5); } //----------------------- Functions to calculate Derivates ---------- double fDeriv1X(double x, double y, double (*f)(double, double)) { double h = fabs(x) > 1 ? 0.01 * x : 0.01; return ((*f)(x+h, y) - (*f)(x-h, y)) / (2 * h); } double fDeriv1Y(double x, double y, double (*f)(double, double)) { double h = fabs(y) > 1 ? 0.01 * y : 0.01; return ((*f)(x, y+h) - (*f)(x, y-h)) / (2 * h); } double fDeriv2X(double x, double y, double (*f)(double, double)) { double h = fabs(x) > 1 ? 0.01 * x : 0.01; return ((*f)(x+h, y) - 2 * (*f)(x, y) + (*f)(x-h, y)) / (h * h); } double fDeriv2Y(double x, double y, double (*f)(double, double)) { double h = fabs(y) > 1 ? 0.01 * y : 0.01; return ((*f)(x, y+h) - 2 * (*f)(x, y) + (*f)(x, y-h)) / (h * h); } // ---------------------- Root-Seeking Class ---------------- class CRoot { public: CRoot(double fTolerX = 1.0e-4, double fTolerY = 1.0e-4); int getIters() { return m_nIter; } double getX() { return m_fx; } double getY() { return m_fy; } bool Bisection(double fXlo, double fXhi, double (*fx)(double x), double (*gx)(double x)); bool FalsePosition(double fXlo, double fXhi, double (*fx)(double x), double (*gx)(double x)); bool BisectionIterative(double fXlo, double fXhi, double fYguessFx, double fYguessGx, int nMaxIter, double (*fx)(double x, double y, int nMaxIter), double (*gx)(double x, double y, int nMaxIter)); bool Newton(double x, double y, int MaxIter, double (*f)(double x, double y), double (*g)(double x, double y)); bool ExNewton(double x, double y, int MaxIter, double (*f)(double x, double y), double (*g)(double x, double y)); protected: int m_nMaxIter; int m_nIter; double m_fx; double m_fy; double m_fTolerX; double m_fTolerY; }; CRoot::CRoot(double fTolerX, double fTolerY) { m_nMaxIter = 10000; m_fTolerX = fTolerX; m_fTolerY = fTolerY; m_nIter = 0; m_fx = 1.e+100; m_fy = 1.e+100; } bool CRoot::Bisection(double fXlo, double fXhi, double (*fx)(double x), double (*gx)(double x)) { double fXm; double fYflo, fYfhi, fYfm; double fYglo, fYghi, fYgm; double fDeltaLo, fDeltaHi, fDeltaM; bool bGoOn; m_fx = 1.e+100; m_fy = 1.e+100; fYflo = (*fx)(fXlo); fYfhi = (*fx)(fXhi); fYglo = (*gx)(fXlo); fYghi = (*gx)(fXhi); fDeltaLo = fYflo - fYfhi; fDeltaHi = fYglo - fYghi; m_nIter = 0; do { fXm = (fXlo + fXhi) / 2; fYfm = (*fx)(fXm); fYgm = (*gx)(fXm); fDeltaM = fYfm - fYgm; if (fDeltaLo * fDeltaM > 0) { fXlo = fXm; fYflo = fYfm; fYglo = fYgm; fDeltaLo = fDeltaM; } else { fXhi = fXm; fYfhi = fYfm; fYghi = fYgm; fDeltaHi = fDeltaM; } bGoOn = fabs(fXlo - fXhi) > m_fTolerX || fabs(fYflo - fYfhi) > m_fTolerY || fabs(fYglo - fYghi) > m_fTolerY; if (++m_nIter > m_nMaxIter) return false; } while (bGoOn); m_fx = fXm; m_fy = (fYflo + fYfhi) / 2; return true; } bool CRoot::BisectionIterative(double fXlo, double fXhi, double fYguessFx, double fYguessGx, int nMaxIter, double (*fx)(double x, double y, int nMaxIter), double (*gx)(double x, double y, int nMaxIter)) { double fXm; double fYflo, fYfhi, fYfm; double fYglo, fYghi, fYgm; double fDeltaLo, fDeltaHi, fDeltaM; bool bGoOn; m_fx = 1.e+100; m_fy = 1.e+100; fYflo = (*fx)(fXlo, fYguessFx, nMaxIter); fYfhi = (*fx)(fXhi, fYguessFx, nMaxIter); fYglo = (*gx)(fXlo, fYguessGx, nMaxIter); fYghi = (*gx)(fXhi, fYguessGx, nMaxIter); fDeltaLo = fYflo - fYfhi; fDeltaHi = fYglo - fYghi; m_nIter = 0; do { fXm = (fXlo + fXhi) / 2; fYfm = (*fx)(fXm, fYguessFx, nMaxIter); fYgm = (*gx)(fXm, fYguessGx, nMaxIter); fDeltaM = fYfm - fYgm; if (fDeltaLo * fDeltaM > 0) { fXlo = fXm; fYflo = fYfm; fYglo = fYgm; fDeltaLo = fDeltaM; } else { fXhi = fXm; fYfhi = fYfm; fYghi = fYgm; fDeltaHi = fDeltaM; } bGoOn = fabs(fXlo - fXhi) > m_fTolerX || fabs(fYflo - fYfhi) > m_fTolerY || fabs(fYglo - fYghi) > m_fTolerY; if (++m_nIter > m_nMaxIter) return false; } while (bGoOn); m_fx = fXm; m_fy = (fYflo + fYfhi) / 2; return true; } bool CRoot::FalsePosition(double fXlo, double fXhi, double (*fx)(double x), double (*gx)(double x)) { double fXm; double fYflo, fYfhi, fYfm; double fYglo, fYghi, fYgm; double fDeltaLo, fDeltaHi, fDeltaM; double fDiffX, fDiffY; bool bGoOn; m_fx = 1.e+100; m_fy = 1.e+100; fYflo = (*fx)(fXlo); fYfhi = (*fx)(fXhi); fYglo = (*gx)(fXlo); fYghi = (*gx)(fXhi); fDeltaLo = fYflo - fYfhi; fDeltaHi = fYglo - fYghi; m_nIter = 0; do { fXm = (fDeltaHi * fXlo - fDeltaLo * fXhi) / (fDeltaHi - fDeltaLo); fYfm = (*fx)(fXm); fYgm = (*gx)(fXm); fDeltaM = fYfm - fYgm; if (fDeltaLo * fDeltaM > 0) { fDiffX = fXlo - fXm; fDiffY = fYflo - fYfm; fXlo = fXm; fYflo = fYfm; fYglo = fYgm; fDeltaLo = fDeltaM; } else { fDiffX = fXhi - fXm; fDiffY = fYfhi - fYfm; fXhi = fXm; fYfhi = fYfm; fYghi = fYgm; fDeltaHi = fDeltaM; } bGoOn = fabs(fDiffX) > m_fTolerX || fabs(fDiffY) > m_fTolerY; if (++m_nIter > m_nMaxIter) return false; } while (bGoOn); m_fx = fXm; m_fy = fYfm; return true; } bool CRoot::Newton(double x, double y, int nMaxIter, double (*f)(double x, double y), double (*g)(double x, double y)) { double f0, g0, fx, fy, gx, gy; double fDiffX, fDiffY; double J; m_fx = 1.e+100; m_fy = 1.e+100; m_nIter = 0; do { f0 = (*f)(x, y); g0 = (*g)(x, y); fx = fDeriv1X(x, y, f); gx = fDeriv1X(x, y, g); fy = fDeriv1Y(x, y, f); gy = fDeriv1Y(x, y, g); J = fx * gy - gx * fy; fDiffX = (f0 * gy - g0 * fy) / J; fDiffY = (g0 * fx - f0 * gx) / J; x -= fDiffX; y -= fDiffY; if (++m_nIter > nMaxIter) return false; } while (fabs(fDiffX) > m_fTolerX || fabs(fDiffY) > m_fTolerY); m_fx = x; m_fy = y; return true; } bool CRoot::ExNewton(double x, double y, int nMaxIter, double (*f)(double x, double y), double (*g)(double x, double y)) { double f0, g0, fx, fy, gx, gy; double fDiffX, fDiffY; double J; m_fx = 1.e+100; m_fy = 1.e+100; m_nIter = 0; do { f0 = (*f)(x, y); g0 = (*g)(x, y); fx = fDeriv1X(x, y, f) - fDeriv2X(x, y, f) / 2 * f0 / fDeriv1X(x, y, f); gx = fDeriv1X(x, y, g) - fDeriv2X(x, y, g) / 2 * g0 / fDeriv1X(x, y, g); fy = fDeriv1Y(x, y, f) - fDeriv2Y(x, y, f) / 2 * f0 / fDeriv1Y(x, y, f); gy = fDeriv1Y(x, y, g) - fDeriv2Y(x, y, g) / 2 * g0 / fDeriv1Y(x, y, g); J = fx * gy - gx * fy; fDiffX = (f0 * gy - g0 * fy) / J; fDiffY = (g0 * fx - f0 * gx) / J; x -= fDiffX; y -= fDiffY; if (++m_nIter > nMaxIter) return false; } while (fabs(fDiffX) > m_fTolerX || fabs(fDiffY) > m_fTolerY); m_fx = x; m_fy = y; return true; } int main(int argc, char* argv[]) { CRoot objRoot; ofstream fout("RootBrak.txt"); double fXlo, fXhi, fX, fY; fXlo = 0; fXhi = 3; cout << "Solving f(x,y) = x^2 + y^2 - (2.1^2 + 3.1^2) " << endl << " g(x,y) = x^2 - y^2 - (2.1^2 - 3.1^2)" << endl << "Exact roots are x = 2.1 and y = 3.1" << endl; fout << "Solving f(x,y) = x^2 + y^2 - (2.1^2 + 3.1^2) " << endl << " g(x,y) = x^2 - y^2 - (2.1^2 - 3.1^2)" << endl << "Exact roots are x = 2.1 and y = 3.1" << endl; cout << "Using Bisection" << endl; fout << "Using Bisection" << endl; cout << "Initial range for x is (" << fXlo << ", " << fXhi << ") " << endl; fout << "Initial range for x is (" << fXlo << ", " << fXhi << ") " << endl; if (objRoot.Bisection(fXlo, fXhi, fx1, gx1)) { cout << "Number of iterations " << objRoot.getIters() << endl; cout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; fout << "Number of iterations " << objRoot.getIters() << endl; fout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; } else { cout << "Method failed!" << endl; fout << "Method failed!" << endl; } cout << "Using False Position" << endl; fout << "Using False Position" << endl; cout << "Initial range for x is (" << fXlo << ", " << fXhi << ") " << endl; fout << "Initial range for x is (" << fXlo << ", " << fXhi << ") " << endl; if (objRoot.FalsePosition(fXlo, fXhi, fx1, gx1)) { cout << "Number of iterations " << objRoot.getIters() << endl; cout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; fout << "Number of iterations " << objRoot.getIters() << endl; fout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; } else { cout << "Method failed!" << endl; fout << "Method failed!" << endl; } fX = 1.5; fY = 1.5; cout << "Using Newton's Method" << endl; fout << "Using Newton's Method" << endl; cout << "Initial guesses for x = " << fX << " and for y = " << fY << endl; fout << "Initial guesses for x = " << fX << " and for y = " << fY << endl; if (objRoot.Newton(fX, fY, 50, f1, g1)) { cout << "Number of iterations " << objRoot.getIters() << endl; cout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; fout << "Number of iterations " << objRoot.getIters() << endl; fout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; } else { cout << "Method failed!" << endl; fout << "Method failed!" << endl; } fX = 1.5; fY = 1.5; cout << "Using Extended Newton's Method" << endl; fout << "Using Extended Newton's Method" << endl; cout << "Initial guesses for x = " << fX << " and for y = " << fY << endl; fout << "Initial guesses for x = " << fX << " and for y = " << fY << endl; if (objRoot.ExNewton(fX, fY, 50, f1, g1)) { cout << "Number of iterations " << objRoot.getIters() << endl; cout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; fout << "Number of iterations " << objRoot.getIters() << endl; fout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; } else { cout << "Method failed!" << endl; fout << "Method failed!" << endl; } cout << endl; pressAnyKey("Press Enter to resume ..."); cout << endl; fout << endl; fXlo = 0; fXhi = 3; cout << "Solving f(x,y) = e^x - y^2 - (e^2 - 4^2) " << endl << " g(x,y) = x^3 - e^y - (2^3 - e^4)" << endl << "Exact roots are x = 2 and y = 4" << endl; fout << "Solving f(x,y) = e^x - y^2 - (e^2 - 4^2) " << endl << " g(x,y) = x^3 - e^y - (2^3 - e^4)" << endl << "Exact roots are x = 2 and y = 4" << endl; cout << "Using Bisection" << endl; fout << "Using Bisection" << endl; cout << "Initial range for x is (" << fXlo << ", " << fXhi << ") " << endl; fout << "Initial range for x is (" << fXlo << ", " << fXhi << ") " << endl; if (objRoot.Bisection(fXlo, fXhi, fx2, gx2)) { cout << "Number of iterations " << objRoot.getIters() << endl; cout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; fout << "Number of iterations " << objRoot.getIters() << endl; fout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; } else { cout << "Method failed!" << endl; fout << "Method failed!" << endl; } cout << "Using False Position" << endl; fout << "Using False Position" << endl; cout << "Initial range for x is (" << fXlo << ", " << fXhi << ") " << endl; fout << "Initial range for x is (" << fXlo << ", " << fXhi << ") " << endl; if (objRoot.FalsePosition(fXlo, fXhi, fx2, gx2)) { cout << "Number of iterations " << objRoot.getIters() << endl; cout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; fout << "Number of iterations " << objRoot.getIters() << endl; fout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; } else { cout << "Method failed!" << endl; fout << "Method failed!" << endl; } fX = 1.5; fY = 1.5; cout << "Using Newton's Method" << endl; fout << "Using Newton's Method" << endl; cout << "Initial guesses for x = " << fX << " and for y = " << fY << endl; fout << "Initial guesses for x = " << fX << " and for y = " << fY << endl; if (objRoot.Newton(fX, fY, 50, f2, g2)) { cout << "Number of iterations " << objRoot.getIters() << endl; cout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; fout << "Number of iterations " << objRoot.getIters() << endl; fout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; } else { cout << "Method failed!" << endl; fout << "Method failed!" << endl; } fX = 1.5; fY = 1.5; cout << "Using Extended Newton's Method" << endl; fout << "Using Extended Newton's Method" << endl; cout << "Initial guesses for x = " << fX << " and for y = " << fY << endl; fout << "Initial guesses for x = " << fX << " and for y = " << fY << endl; if (objRoot.ExNewton(fX, fY, 50, f2, g2)) { cout << "Number of iterations " << objRoot.getIters() << endl; cout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; fout << "Number of iterations " << objRoot.getIters() << endl; fout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; } else { cout << "Method failed!" << endl; fout << "Method failed!" << endl; } cout << endl; pressAnyKey("Press Enter to resume ..."); cout << endl; fout << endl; fXlo = 1; fXhi = 5; double fYguessFx = 8; double fYguessGx = 0; int nMaxIter = 100; cout << "Solving f(x,y) = x^2 - y^2 * e^x -(3^2 - 5^2 * e^3) " << endl << " g(x,y) = e^x - y - (e^3 - 5)" << endl << "Exact roots are x = 3 and y = 5" << endl; fout << "Solving f(x,y) = x^2 - y^2 * e^x -(3^2 - 5^2 * e^3) " << endl << " g(x,y) = e^x - y - (e^3 - 5)" << endl << "Exact roots are x = 3 and y = 5" << endl; cout << "Using Bisection" << endl; fout << "Using Bisection" << endl; cout << "Initial range for x is (" << fXlo << ", " << fXhi << ") " << endl; fout << "Initial range for x is (" << fXlo << ", " << fXhi << ") " << endl; if (objRoot.BisectionIterative(fXlo, fXhi, fYguessFx, fYguessGx, nMaxIter, fx3, gx3)) { cout << "Number of iterations " << objRoot.getIters() << endl; cout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; fout << "Number of iterations " << objRoot.getIters() << endl; fout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; } else { cout << "Method failed!" << endl; fout << "Method failed!" << endl; } fX = 2.5; fY = 2.5; cout << "Using Newton's Method" << endl; fout << "Using Newton's Method" << endl; cout << "Initial guesses for x = " << fX << " and for y = " << fY << endl; fout << "Initial guesses for x = " << fX << " and for y = " << fY << endl; if (objRoot.Newton(fX, fY, 50, f3, g3)) { cout << "Number of iterations " << objRoot.getIters() << endl; cout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; fout << "Number of iterations " << objRoot.getIters() << endl; fout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; } else { cout << "Method failed!" << endl; fout << "Method failed!" << endl; } fX = 2.5; fY = 2.5; cout << "Using Extended Newton's Method" << endl; fout << "Using Extended Newton's Method" << endl; cout << "Initial guesses for x = " << fX << " and for y = " << fY << endl; fout << "Initial guesses for x = " << fX << " and for y = " << fY << endl; if (objRoot.ExNewton(fX, fY, 50, f3, g3)) { cout << "Number of iterations " << objRoot.getIters() << endl; cout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; fout << "Number of iterations " << objRoot.getIters() << endl; fout << "Roots are at x = " << objRoot.getX() << " and y = " << objRoot.getY() << endl; } else { cout << "Method failed!" << endl; fout << "Method failed!" << endl; } cout << endl; pressAnyKey("Press Enter to end the program"); fout.close(); return 0; } //--------------------------------------------------------------------
Copyright (c) Namir Shammas. All rights reserved.