Lesson 08

Root-Finding

INTRODUCTION

In this lesson we will explore methods for finding numerical solution(s) to the equation $f(x) = 0$. Note that this is identical to finding solutions to an equation $g(x) = A$, since we can define a new function $f(x) = g(x) - A$ so that we are again looking for the value of $x$ such that $f(x) = 0$. The problem of computing the range of a projectile subject to linear drag is an example of how one often ends up with a transcendental equation rather than an analytic solution.

Projectile Motion The motion of a projectile is known to have parabolic shape. This is true only in first approximation, when the air resistance is neglected. In this case the equations of motion are written as:

$\ddot{x} = 0\\ \ddot{y} = -g$

for which we can easily solve for the analytic solution

$x(t) = x_0+v_{x,0}t\\ y(t) = y_0+v_{y,0}t-\frac{1}{2}gt^2$

Solving for the time in the $x$ equation and substituting into the $y$ equation yields a solution for the distance at which the projectile hits the ground $R$:

$R=\frac{v_0^2}{g}\sin(2\theta)$

where $v_0$ is the initial velocity of the projectile and $\theta$ the angle of the cannon with respect to the horizontal direction. Air resistance, or drag, will change this answer. Will it make the range larger or smaller?

Air Resistance is a retarding force. For relatively small velocities, air resistance is proportional (and opposed to) the velocity vector and the equations of motion modify into:

$\ddot{x} = -k\dot{x}\\ \ddot{y} = -g-k\dot{y}$

where $k$ is a constant that depends on the size of the projectile and the air viscosity. These equations have analytic solutions:

$x(t) = x_0+\frac{v_{x,0}}{k}(1-e^{-kt})\\ y(t) = y_0-\frac{gt}{k}+\frac{kv_{y,0}+g}{k^2}(1-e^{-kt})$

To find the time at which the projectile hits the ground, one can set $y(t)=0$, yielding a landing time $T$ given by the implicit equation:

$T=\frac{kv_{y,0}+g}{gk}(1-e^{-kT})$

which is a trascendental equation without an analytic solution. We need to find a way to solve numerically the equation:

$f(T) = T-\frac{kv_{y,0}+g}{gk}(1-e^{-kT})=0$

A Program Function provides a convenient way to repeatedly calculate a function without having to retype the equation in your code. 

def root(time):
   z = time - (k*vy0 + g)/(g*k) * (1.0 - exp(-k*time))
   return z;

You can then compute the value of f(T) in your program by calling this function with a given value of T:

guess = root(12.0)

Note however that you would have to set the values of the constants within this function as they are treated as local variables - they are only known inside the function.

Find That Root

How would you go about finding the value for $T$ that satisfies the above equation? A simple, but highly inefficient method is to start with a guess, and take small steps in time until the sign of $f(T)$ changed. The root, or value of $T$ where $f(T)=0$ must therefore be inbetween the last two values.

Challenge: What might you use for your initial guess of the landing time? How would you determine which direction to start changing - smaller or larger time?

The problem with this method is that the accuracy of your answer is determined by your step size, but if you choose a small step size you are forced to take many steps (unless your initial guess was really good).

The Bisection Method is a simple way of decreasing your step size in order to get arbitrary accuracy. It is a simple, but robust method. The process starts by identifying two points $x=a$ and $x=b$ that bound the root you are looking for. If $f(a)$ and $f(b)$ have opposite signs, then there must be a (or multiple) value(s) of $x$ inbetween $a$ and $b$ such that $f(x) = 0$. What might you choose for bounding values in our projectile problem?

The idea of the bisection method is to now bisect the space defined by our bounds by finding the midpoint, $c=(a+b)/2$ and evaluating the function at this midpoint. If $f(c)$ has the same sign as $f(a)$, then we know there is a root between $c$ and $b$, and if $f(c)$ has the opposite sign as $f(a)$, then we know there is a root between $a$ and $c$. You then use $c$ to replace the appropropriate bounding value and repeat the procedure with the new bounds. You thus cut your error in half with each step. This routine is very robust but converges slowly (linearly).

Exercise: write a routine that finds the landing time and distance for $v_{0}=100$ m/s, $\theta=45^\circ$, and $k=0.01$.

Newton's Method converges faster than the bisection method (reaching a comparable error in fewer steps), but requires that you know the derivative of your function. If $x_0$ is the root we are searching for and $x_1$ is our initial guess for the root, an approximate expression for the derivative would be

$f'(x_1) = \frac{f(x_1) - f(x_0)}{x_1 - x_0}$

We know $f(x_0)=0$, so we can turn this equation around to solve for $x_0$. However, since this equation is just an approximation, this won't be the exact value but rather a better approximation that we will call $x_2$.

$x_2 = x_1 - \frac{f(x_1)}{f'(x_1)}$

We can then iterate this equation over and over again, replacing our old guess $x_1$ with our new guess $x_2$ and recomputing the above equation. This method converges quadratically.

Challenge: What can go wrong? Can you think of a function and an initial guess that would make Newton's method diverge from the true root of the function?

Exercise: write a routine that uses Newton's method to find the landing time and distance for $v_{0}=100$ m/s, $\theta=45^\circ$, and $k=0.01$.

Assignment: Submit your bisection code that finds the range of a projectile for $v_{0}=100$ m/s, $\theta=45^\circ$, and $k=0.01$ with an accuracy of one part in 106.