# Questions on Numerical Solution Methods

Why do we bother with Newton Iterations when there are better way to solve things.

This is one of those leaps of faith that we ask you to make in the Education Business, and then we wonder why you didn't absorb it. The fact is that in a large range of applications, Newton (or some close relation) is the best way to solve things. Mathematica uses a Newton iteration when you invoke FindRoot. I use a generalized Newton iteration on a daily basis to solve systems of several thousand coupled non-linear equations. No Newton is not the universal best way to solve equations. You will have to learn several methods, if you get heavily involved in the equation solution business, but this is a good start to your education in the field.

How can a computer generate a random number when there is nothing random about the way it operates?

I'm not the best guy to ask this question, because I can't give you a complete answer. There are large numbers of papers in journals on the subject if you are interested, and source code available for "random number" generators. The simple answer is that the don't generate random numbers. They generate "pseudo-random" numbers. These are sequences of numbers that proceed from some starting point (called the seed) and via a well determined algorithm produce a series of numbers that for most purposes look random. Start with the same seed, generate the same sequence of numbers. You can introduce a certain amount of true randomness in the whole process with your selection of the seed. Pull out the number of milliseconds that have elapsed during the day and scale by an appropriate number.

Could you give us some tips on sorting arrays alphabetically?

Aside from declaring the right variables CHARACTER, very little changes from a sort of numbers. You can continue to use the same methods with the same .LT. or .GT. (maybe .LE. or .GE.) tests. The secret is understanding how Fortran views the order or character strings. The best way to understand this ordering this is to play with the names in "grades.in" once you have your program working. Generally 'AA'<'Aa'<'Aaa'<'Baa' <'Bba'<'Zba'<'aAa'<'aaa'<'zaa'.

Can you rate how expensive (time consuming) the various sort techniques are?

You can't predict time consumption without knowing something about the data. As a general rule bubble comes in last. See the next question.

Is there ever a case when it is very bad to use a particular sorting type?

Yes, if you have a large number (million, billion, ...) of numbers or names that are randomly ordered, then you should not use bubble, selection, or insertion sort. Use the NETLIB routine or one of its relatives that take of order n log(n) computer operations rather than n2 operations. If you have a list that was recently sorted, and a couple of new elements are tagged on the end, but out of order, then you do not want to do a selection sort, and probably shouldn't use a bubble or the NETLIB sort in grades.f. Something related to an insertion sort (maybe with a halving search that I'll show you later) is a very good idea.

Why doesn't Fortran have a built in Sort routine?

It is neither common enough nor a simple enough to rate as an intrinsic function. You will find other more frequently used procedures that fall into this category. That is why NETLIB was developed. Having a subroutine that you can plug into your program is nearly as easy as having an intrinsic function, and generally a lot more flexible.

Can you always trust NETLIB files?

No, but they almost always warn you if a package is inadequately tested. In general these subroutines are as well tested as any software that you will find. Packages like LINPACK and BLAS are at the zero bug level.

How would you find a subprogram in NETLIB if you didn't know its name.

In NETSCAPE select the line that searches the NETLIB data base by topic. You can also look for a file with a name beginning "index". Also check the file "linpack.doc" in ~jhm for a description of the LINPACK systematic naming convention, giving another approach to finding some subroutines.

When is Gauss Elimination Used?

The simple answer is: When you are going to solve a set of linear equations without returning for another solution with a different right hand side, but the same coefficients for the unknowns. Your question is actually deeper than you may suspect. At some point, depending on the number and location of non-zero coefficients, and the total number of equations, Gauss Elimination becomes too expensive in both computer time and memory requirements. Then you must shift to iterative solution methods. I will not be teaching these methods. However, it is worth knowing that your first introduction to iterative solution of linear equations will probably be the Gauss-Seidel method. This method is relatively simple to teach, but not used frequently in modern applications.

What is the purpose of having an additional array that carries along the values of an original array?

I assume you are referring to my use of arrays "a" and "at" in some of my examples. The reason is that when you pass an array "a" to a LINPACK subroutine as part of a solution of a set of linear equations, the subroutine overwrites the original values in "a". The extra array "at" is needed only if you have further use for those original coefficients. I needed them for a cross-check of results.

How Deep will we go into Differential Equations?

We will look at a simple first order system of differential equations. No fancy theorems. We will follow the steps to convert the differential equations to algebraic equations and proceed with a solution. If time permits, I may discuss a second order ODE.

Could you list methods (Secant, Adams, ...) in order of Accuracy.

Decouple Secant from Adams, Runge-Kutta, etc. Secant like Newton is used to solve non-linear algebraic equations (yes those algebraic equations could be generated by an implicit Adam's-Moulton, but that's a second level of solution detail). For a given number of iterations, Secant usually produces a less accurate answer than Newton. On the ODE side, Explicit and Implicit Euler formally have the same first order (lowest) accuracy. Crank-Nicholson has the next higher (second order) accuracy. Runge-Kutta, Adams-Bashford, Adams-Moulton, and Implicit Adams-Moulton are all tied for the highest formal accuracy (fourth order) of those that we have studied. In practice Runge-Kutta is usually a little more accurate than the others for a given step size. This accuracy has meaning only as long as the methods stay numerically stable. The ordering of methods from first to last to give up as step size is increased is usually (but not always): Explicit Euler or Adams-Bashford give up first ; Adams-Moulton holds on until a little larger step size; Implicit Adams-Moulton and Runge-Kutta last further; and Implicit Euler stays stable to any time step size (results may be very inaccurate, but not fatal to the floating point unit)

What is the difference between the Secant and Newton method.

Just the way the derivative term is evaluated. In Newton you evaluate the derivative analytically (standard calculus). With Secant you are doing a numerical approximation to the derivative.

Why does it seem like ODE's are only solved numerically with a computer. Why can't they be solved analytically like Mathematica can do?

The number of equations that Mathematica can solve analytically is actually relatively limited. When you get to the ones involved in really interesting physical models of complex systems, neither Mathematica nor the best Mathematicians can give you a useful analytic solution.

Could it sometimes be better to fit a curve and then solve it analytically rather than numerically?

You've gotten to the heart of many numerical methods. The act of fitting a curve is usually an approximation so that you are not producing a true analytic solution to the original differential equation. The use of fit curves generate algebraic equations that can be solved relatively with computers (LINPACK, Newton, ...). This idea is in a sense at the heart of Finite Element methods. You can take quite a few courses on that subject.

If we could only remember a couple sentences about ODE solutions what would they be?

Use existing solution packages whenever possible from places like NETLIB (ODEPACK). However, test them carefully for limits of accuracy and stability. One good class of tests is to make a series of runs, each with different integration steps.

## Up one level / Home

Maintained by John Mahaffy : jhm@cac.psu.edu