On the other hand, finding the roots of an equation containing a combination of powers, exponential functions, and trigonometic functions is much more difficult.
In addition, since machine accuracy is always limited, a method which uses a large number of arithmetic operations always runs the risk of introducing and propogating so many inaccuracies, that the final answer may be meaningless. (Math 166, Numerical Analysis, addresses some of the issues regarding accuracy in numerical methods.)
One approach commonly taken to solve numerical problems is by iteration, in which approximate answers are repeatedly obtained using some method. The iteration stops when two successive approximate answers are so close that they can be considered (for the sake of the problem) to be identical.
We can demonstrate this approach using Newton's method for finding roots (as found in Calculus books).
It is good to do some preliminary analysis first. We look at the function and try to determine where the roots lie. Since this function is a cubic and not a quadratic, we cannot use the quadratic formula.
We can first try to graph the function. To do this, it helps to learn where the maximum and minimum points are. Thus we need to take the derivative.
We compute the derivative of f(x), f'(x) = 3x2 - 6x. Setting this equal to zero, we find: 3x2 - 6x = 3x(x-2) = 0 which implies that x=0,2.
Next we find the values of f(x) at x=0,2. We learn that f(0) = 1 and f(2) = 8 - 12 + 1 = -3.
Looking at the horizontal asymptotes, if x goes to infinity, f(x) goes to positive infinity, and if x goes to negative infinity, f(x) goes to negative infinity.
Since we have discovered points where the function values alternate in sign, this analysis indicates that the roots must lie between (1) negative infinity (negative function value) and 0 (positive function value), (2) 0 and 2 (negative function value), and (3) 2 and positive infinity (positive function value).
STEP 1. Make an educated guess for the first approximation to the solution x.
For our purposes, let us choose x = 0.5
STEP 2. Apply Newton's Formula.
Remember that Newton's Formula is based on the slope of a straight line. The actual recursive formula is:
xi = xi-1 - f(xi-1)/f'(xi-1)
For our case we have
xi = xi-1 - (xi-13 -3xi-12+1) /(3xi-12 - 6xi-1)STEP 3. Repeat Step 2 until the desired accuracy is achieved.
In our example,
x0 = .5When two successive answers are essentially the same, we can stop Newton's iteration and assume that this common answer is our desired answer.
x1 = .6667
x2 = .6528
x3 = .6527
x4 = .6527
We will therefore code the program using two functions, one for f(x) and one for the derivative f'(x).
We also need to have a loop to repeat the iterative process with some sort of exit conditiion when the old approximation and the new approximation are "close enough".
program newton
implicit none
double precision :: x0, x1, epsilon, f, dfdx
epsilon = 1.0d-14
x0 = 0.5d00
do
x1=x0-(f(x0)/dfdx(x0))
write(6,*) "Root = ", x1
if (abs (x1-x0) < epsilon ) exit
x0=x1
end do
write(6,*) "Root = ", x1
end
!
double precision function f(x)
implicit none
double precision :: x
f = x**3 - 3*(x**2) + 1
return
end function f
!
double precision function dfdx(x)
implicit none
double precision :: x
dfdx = 3*x**2 - 6*x
return
end function dfdx
Problem 2: We might not get a root. The algorithm might either oscillate in approximate values or actually diverge!
For example, consider f(x) = x1/3.
The root is 0.
(Step 1.) Choose x0 = 1.
(Step 2.) xi = xi-1 -
xi-11/3/[(1/3)xi-1-2/3]
= xi-1 - 3xi-1 = -2xi-1
Thus we have:
x0 = 1This oscillates, diverging toward both positive infinity and negative infinity, even though the root is, simply, 0.
x1 = -2
x2 = 4
x3 = -8
This page is maintained by Dennis C. Smolarski, S.J.
dsmolarski@math.scu.edu
© Copyright 2000 Dennis C. Smolarski, SJ, All rights reserved.
Last changed: 16 February 2000. Minor correction 29 February 2000.