Solving Equations

Numerical Methods

David Mayerich

Scalable Tissue Imaging and Modeling (STIM) Laboratory

Department of Electrical and Computer Engineering

Cullen College of Engineering

University of Houston

David Mayerich

STIM Laboratory, University of Houston

Solving Equations

Solving Equations Algebraically

Solving Quadratic Equations

Solutions as Roots

David Mayerich

STIM Laboratory, University of Houston

Solving Equations Algebraically

  • We often require the dependent variable \(x\) that satisfies an equation

David Mayerich

STIM Laboratory, University of Houston

y=f(x)

where the function \(f\) and \(y\) are known

  • This can be accomplished by calculating the inverse function:

x=f^{-1}(y)

Implementing Solutions

David Mayerich

STIM Laboratory, University of Houston

  • The inverse function \(f^{-1}(z)=\sin^{-1}(z) + e\) could be implemented as a series:

\sin^{-1}(x) = x + \frac{x^3}{6} + \frac{3x^5}{40} + O\left(x^6\right)
f^{-1}(z) \approx x + \frac{x^3}{6} + \frac{3x^5}{40} + e

  • Consider the equation: \(y = \sin(x - e)\)

\sin^{-1}(y) = x-e
\sin^{-1}(y) + e = x
f(z)=\sin(z-e)
f^{-1}(z)=\sin^{-1}(z) + e

Computer Algebra Systems

David Mayerich

STIM Laboratory, University of Houston

  • How can we get the inverse formulation: \(f(z)=\sin(z-e) \rightarrow f^{-1}(z)=\sin^{-1}(z) + e\)?

  • Computer Algebra Systems (CAS) can get us there:

from sympy import *

y, z = symbols("y z")
solve(sin(z - exp(1)) - y, z)

[{z: asin(y) + E}, {z: -asin(y) + E + pi}]

Numerical Solutions

  • Consider functions without inverses

David Mayerich

STIM Laboratory, University of Houston

y = \sin\left(z^2 - e^z\right)

  • Wolfram Alpha cannot calculate the inverse but can generate a graph

  • This is done by expressing the equation as a root:

(the vast majority of them)

\sin\left(z^2 - e^z\right) - y = 0

  • Finding the roots of this equation for a given \(y\) gives the solution \(z\)

Root-Finding Formulas

  • Find the solution to \(y = \sin\left(z^2 - e^z\right)\) for \(y = -0.5\)

David Mayerich

STIM Laboratory, University of Houston

\(z=-0.3909\)

  • Express the special functions as series:

\sin(z)\approx z
e^z \approx 1 + z + \frac{z^2}{2}
\sin\left(z^2 - e^z\right)\approx \frac{z^2}{2} - z - 1
\frac{z^2}{2} - z - 1 - y=0
\frac{1}{2}z^2 - z - 0.5=0
z=\frac{-b\pm \sqrt{b^2 - 4ac}}{2a}
z=\frac{1\pm \sqrt{(-1)^2 - 4(0.5)(-0.5)}}{2(0.5)} = [\ -0.4142,\ 2.414\ ]
\epsilon_r\approx 5.9\%

Simple Roots

David Mayerich

STIM Laboratory, University of Houston

  • Consider an interval \([a,b]\) of a function \(f(x)\in\mathbb{R}\) where \(f(a)\) and \(f(b)\) have opposite signs:

f(a)\cdot f(b) < 0

  • If \(f(x)\) is continuous on \([a,b]\), then there exists a number \(c\in[a,b]\) such that \(f(c)=0\)

  • The value \(x=c\) is a simple root of \(f\)

Bisection Method

Binary Searches for Roots

Bisection Algorithm

Error

David Mayerich

STIM Laboratory, University of Houston

Bisection Method

David Mayerich

STIM Laboratory, University of Houston

a
c
b
y=0

Bisection Method

David Mayerich

STIM Laboratory, University of Houston

y=0
a
b
c

Bisection Method

David Mayerich

STIM Laboratory, University of Houston

y=0
a
b
c

Bisection Method

David Mayerich

STIM Laboratory, University of Houston

y=0
a
b
c

Bisection Method

David Mayerich

STIM Laboratory, University of Houston

y=0
a
b
c

Bisection Method

David Mayerich

STIM Laboratory, University of Houston

y=0
a
b
c

Bisection Method

David Mayerich

STIM Laboratory, University of Houston

y=0
f(c) \approx 0

Bisection Method Evaluation

  • Calculate the solution to \(y=\sin\left(x^2 - e^x\right)\) for \(y=-0.5\) given the interval \([-1.0, -0.1]\) out to \(4\) digits of precision

David Mayerich

STIM Laboratory, University of Houston

\sin\left(x^2-e^x\right) = -0.5
-0.5500
0.2290
-0.3250
-0.07851
-0.4375
0.06122
-0.3813
-0.01206
-0.4094
0.0237015
-0.3954
0.005665
-0.3884
-0.003190
-0.3919
0.001224
-0.3902
-0.0009235
-0.3911
0.0002126
-0.3907
-0.0002925
-0.3909
-0.00004002
-1.000
1.091
-0.1000
-0.2801
-0.5500
0.2290
-0.1000
-0.2801
-0.5500
0.2290
-0.3250
-0.07851
-0.4375
0.06122
-0.3250
-0.07851
-0.4375
0.06122
-0.3813
-0.01206
-0.4094
0.0237015
-0.3813
-0.01206
-0.3954
0.005665
-0.3813
-0.01206
-0.3954
0.005665
-0.3884
-0.003190
-0.3919
0.001224
-0.3884
-0.003190
-0.3919
0.001224
-0.3902
-0.0009235
-0.3911
0.0002126
-0.3902
-0.0009235
-0.3911
0.0002126
-0.3907
-0.0002925
\sin\left(x^2-e^x\right)+0.5 = 0
\mathbf{a}
\mathbf{b}
\mathbf{f(a)}
\mathbf{f(b)}
\mathbf{c}
\mathbf{f(c)}

Bisection Method Algorithm

David Mayerich

STIM Laboratory, University of Houston

// implement the bisection method where:
// f is a function pointer that takes a parameter x and returns a value
// a and b are brackets bounding a root
// N is the maximum number of iterations
// return an approximation to the root bounded by a and b
double bisection(double (*f)(double), double a, double b, int N){
  float c;

  return c; // return the last approximation to the root
}

Bisection Method Algorithm

David Mayerich

STIM Laboratory, University of Houston

// implement the bisection method where:
// f is a function pointer that takes a parameter x and returns a value
// a and b are brackets bounding a root
// N is the maximum number of iterations
// return an approximation to the root bounded by a and b
double bisection(double (*f)(double), double a, double b, int N){
  float c;
  for(int i = 0; i < N; i++){ // iterate N times

  }
  return c; // return the last approximation to the root
}

Bisection Method Algorithm

David Mayerich

STIM Laboratory, University of Houston

// implement the bisection method where:
// f is a function pointer that takes a parameter x and returns a value
// a and b are brackets bounding a root
// N is the maximum number of iterations
// return an approximation to the root bounded by a and b
double bisection(double (*f)(double), double a, double b, int N){
  float c;
  for(int i = 0; i < N; i++){ // iterate N times
    c = (a + b) / 2.0; // find midpoint to approximate f(c) = 0
    if(f(c) == 0) return c; // if f(c) actually IS zero, return it (rare)
    if(f(a) * f(c) > 0) a = c; // if f(c) has the same sign as f(a), move a
    else b = c; // otherwise move b
  }
  return c; // return the last approximation to the root
}

Bisection Method - Example

  • Calculate the solution to \(y=\sin\left(x^2 - e^x\right)\) for \(y=-0.5\) given the interval \([-1.9, 1.7]\)

David Mayerich

STIM Laboratory, University of Houston

double func(double z) {
	return = sin(z*z - exp(z));
}

Bisection Method - Example

  • Calculate the solution to \(y=\sin\left(x^2 - e^x\right)\) for \(y=-0.5\) given the interval \([-1.9, 1.7]\)

David Mayerich

STIM Laboratory, University of Houston

double func(double z) {
	return = sin(z*z - exp(z));
}

int main(int argc, char** argv) {
	int N = 100;
    
    double result = bisection(func, -1.9, 1.7, N);

}

Bisection Method Theorem

If the bisection method is applied to a continuous function \(f\) on an interval \([a, b]\), where \(f(a)f(b)<0\), then after \(n\) steps an approximate root will be computed with an error at most \(\epsilon = \frac{b-a}{2^n}\).

David Mayerich

STIM Laboratory, University of Houston

  • Each iteration computes the root with one additional binary bit of precision
     

  • How many iterations are required to guarantee an error less than \(\epsilon'\)?

\frac{b-a}{2^n} < \epsilon'
b-a < 2^n \epsilon'
\log(b-a) < n \log 2 + \log \epsilon'
\log(b-a) - \log \epsilon' < n \log 2
n > \frac{\log(b-a) - \log \epsilon'}{\log 2}
n = \left\lceil\frac{\log(b-a) - \log \epsilon'}{\log 2}\right\rceil

Bisection Method Characteristics

  • Very robust - guaranteed to converge to a root as long as \(f(a)f(b)<0\)
     

  • Does not work for multiple roots when \(f(a)f(b) \geq 0\)
     

  • Binary search reduces error by \(O(2^n)\) after \(n\) iterations

    • This is equivalent to \(1\) bit per iteration

    • \(\approx 23\) iterations for float, \(\approx 53\) iterations for double

David Mayerich

STIM Laboratory, University of Houston

Method of False Position

Combining Function Size and Magnitude

Method of False Position Algorithm

Problems

David Mayerich

STIM Laboratory, University of Houston

Method of False Position

David Mayerich

STIM Laboratory, University of Houston

a
c
b
y=0

Method of False Position

David Mayerich

STIM Laboratory, University of Houston

a
c
b
y=0

Method of False Position

David Mayerich

STIM Laboratory, University of Houston

a
c
b
y=0

Method of False Position

  • Instead of using the midpoint for \(c\), use the x-intercept of \(\overline{f(a)f(b)}\)

  • Given the bracket \([a, b]\), we calculate two points: \(\left( a, f(a)\right)\) and \(\left( b, f(b)\right)\)

  • The line between them is specified by:

David Mayerich

STIM Laboratory, University of Houston

y=mx+b
y(x)=\frac{f(a)-f(b)}{b-a} (x-a) + f(a)

  • Setting \(y(x)=0\) and simplifying yields the root of the line at:

x=\frac{a f(b) - b f(a)}{f(b)-f(a)}

Method of False Position Algorithm

David Mayerich

STIM Laboratory, University of Houston

// implement the method of false position where:
// f is a function pointer that takes a parameter x and returns a value
// a and b are brackets bounding a root
// N is the maximum number of iterations
// return an approximation to the root bounded by a and b
double mfp(double (*f)(double), double a, double b, int N){
  double c;
  
  return c; // return the last approximation to the root
}

Method of False Position Algorithm

David Mayerich

STIM Laboratory, University of Houston

// implement the method of false position where:
// f is a function pointer that takes a parameter x and returns a value
// a and b are brackets bounding a root
// N is the maximum number of iterations
// return an approximation to the root bounded by a and b
double mfp(double (*f)(double), double a, double b, int N){
  double c;
  for(int i = 0; i < N; i++){ // iterate N times

  }
  return c; // return the last approximation to the root
}

Method of False Position Algorithm

David Mayerich

STIM Laboratory, University of Houston

// implement the method of false position where:
// f is a function pointer that takes a parameter x and returns a value
// a and b are brackets bounding a root
// N is the maximum number of iterations
// return an approximation to the root bounded by a and b
double mfp(double (*f)(double), double a, double b, int N){
  double c;
  for(int i = 0; i < N; i++){ // iterate N times
    double fa = f(a); // precompute f(a) and f(b)
    double fb = f(b);
    c = (a * fb – b * fa) / (fb – fa);// WAS (a + b) / 2
    if(f(c) == 0) return c; // if f(c) actually IS zero, return it (rare)
    if(f(a) * f(c) > 0) a = c; // if f(c) has the same sign as f(a), move a
    else b = c; // otherwise move b
  }
  return c; // return the last approximation to the root
}

Bracketing Convergence

Convergence Analysis

Bisection and False Position

David Mayerich

STIM Laboratory, University of Houston

Bisection Method Behavior

  • What we know given the bracket \([a, b]\) on a continuous function \(f(x)\) if \(a \leq b\):

    • At each iteration \(n\), the update for \(a_n\) will always increase:
                             \(a\rightarrow a_0 \leq a_1 \leq a_2 \leq \cdots \leq a_n\)

    • At each iteration \(n\), the update for \(b_n\) will always decrease:
                             \(b\rightarrow b_0 \geq b_1 \geq b_2 \geq \cdots \geq b_n\)

    • As \(n\) increases, we expect both sides of the bracket to converge to the root \(r\):

       

    • From the Bisection Method Theorem, we know that:

David Mayerich

STIM Laboratory, University of Houston

\lim_{n\rightarrow \infty} a_n = \lim_{n\rightarrow \infty} b_n = r
b_n - a_n = \frac{b_0 - a_0}{2^n}

Bisection Convergence Rate

  • For any iteration \(n\geq 1\), the Bisection Method calculates a new value \(c_n\):

David Mayerich

STIM Laboratory, University of Houston

|r - c_n| \leq \frac{b_n - a_n}{2} = \frac{b - a}{2^n}
|r - c_n| \leq \left( \frac{1}{2} \right)^n |b-a|

  • The convergence rate specifies how quickly the error reduces across iterations:

\lim_{n\rightarrow\infty} \frac{|r - c_n|}{|r-c_{n-1}|} = \frac{\left(\frac{1}{2}\right)^{n}|b-a|}{\left(\frac{1}{2} \right)^{n-1}|b-a|}
= \frac{\left(\frac{1}{2}\right)^{n}}{\left(\frac{1}{2} \right)^{n-1}}
= \frac{1}{2}