this is mth:newton.hlp as of 26 september 1986, a help file for mth:newton, a program which carries out newton iteration on a system of up to 10 equations. newton was written by john burkardt staff programmer mathematics department university of pittsburgh pittsburgh, pa ********************************************************************* newton can solve two distinct but related groups of problems. the first type of problem involves n equations in n unknowns. you supply a single approximate guess, and newton iteration is used to find the solution. the second type of problem involves n-1 equations in n unknowns. this is known as the 'continuation process'. in this case, the underdetermined system now defines a curve of points. you supply an approximate point on the curve, and then use newton iteration to get onto the curve. now you have the program take a step in the tangent direction to find another approximate point. newton iteration is performed on this point to get back on the curve. by repeated applications of this process, you can find a sequence of points along the curve. note that in this case, the program adds an extra equation, so that the 'number of equations' must be increased by 1. also, because of the crude nature of this algorithm, there may be some problems. there is no stepsize algorithm in the code. that means you must choose a stepsize h to use, which may be too large, causing newton to fail. also, the form of the extra equation may not work well with some problems. again, you will most likely notice this when newton does not converge well. ********************************************************************* to use newton, you must first write two subroutines which define your problem. the first of these must be named funk and must have the following form subroutine funk(fx,neqn,x) dimension fx(neqn) dimension x(neqn) from time to time, the program will call this subroutine, and will give you the value of x. it expects you to assign values to fx and return. if you are using the continuation option, do not bother assigning a value to fx(neqn). that is done for you. the second subroutine evaluates the jacobian matrix of partial derivatives of the function defined in funk. the subroutine must have the following form subroutine jacob(fprime,nrow,neqn,x) dimension fprime(nrow,neqn) dimension x(neqn) from time to time, the program will call jacob, and supply the value of x. you must evaluate the entries of the jacobian where fprime(i,j)=d f(i)/d x(j). note that if you decide to use the finite difference jacobian option, you do not have to perform any calculations in jacob, but you still must supply a subroutine of that name. if you are using the continuation option, do not bother assigning values to the last row of fprime. that will be done for you. ********************************************************************* running newton once you have written your subroutines, and supposing that they are both stored in one file named 'myprog.for', you can run newton by typing the following command execute myprog.for, newton.for[??????,??????] the program will first ask you whether you want to use the standard method, or the continuation option. then it will ask you for the number of equations in your system. then it will want the current starting x. after that, the program will ask you again and again what the next command is, and will carry it out. the following commands are available to you - b = begin newton iteration c = continue current newton iteration d = set damping option e = set number of equations (aborts current newton iteration) f = set finite difference option g = check jacobian versus finite difference jacobian (aborts current newton iteration) h = set continuation stepsizes i = set starting value of predicted next continuation point. j = print jacobian k = define how often jacobian is evaluated l = list commands m = set maximum number of newton iterations n = set the norm to use o = change amount of output p = print x, f(x) (and tan(x) for continuation problems) s = stop t = set error tolerances v = print values of control quantities w = print work statistics x = set current x approximation y = set target for continuation (aborts current newton iteration) the meaning of these commands is startup commands e = means that you are declaring the number of equations. since this question is asked at the beginning, you probably should not need to use this command thereafter. a clever programmer, however, could run several problems in one run of different sizes using the size as a flag in the function and jacobian routines. x = means that you input a point x from which all newton iterations will begin. the improved x produced by newton iteration can be used to proceed via the 'c' command, but any 'b' command causes the program to begin from this reference starting point. you can change x at any time. y = means that in your continuation process, you would like the program to search for a special point. this special point is defined by the value of a given index. if you were tracing a circle of radius 9, you might want to find the point (3,0). to do this, you would request the point with index 1 equal to 3, or with index 2 equal to 0. note however that the program only checks for these points by comparing the two most recent points it computed and seeing if the given value lies between them. thus, on a circle, the two points might be (2.98, -0.1) and (2.97, +0.2). in that case, the program will notice that the point with x(2)=0 lies between these two points, but it would not realize that a point with x(1)=3 lies between them. so be careful how you specify this target point. print out commands p = causes the program to print the value of the starting point and the function f(x). j = causes the program to print the value of the jacobian matrix at the current point. if no newton steps have been taken, this is the starting point. if newton steps have been taken, the jacobian is displayed for the point of last evaluation, which is the point just before the current iterate, unless the modified newton method has been used. commands for the continuation step h = set continuation stepsize. if the program has a point x, it will compute the tangent to the curve at x, tan. to get the next approximate point on the curve, it will set xn new = x + h*tan. thus, for small h, the new point should be close to the curve, and for large h, probably far off. since this point is the starting point of a newton iteration, far away points caused by large h values may cause newton to fail. you may set h negative to reverse direction. i = take continuation step. this command causes the program to set the starting point xn new = x + h*tangent. it also means that the program will inform you of the extra equation, which will have the form fx(n)=x(ihold)-xhold. the program picks values for ihold and xhold each time. ihold is called the continuation parameter. the effect of this equation is to demand that the newton iterates all fix the ihold-th index of x. once you have issued the i command, it is up to you to use the b and c commands to get a point on the curve. thus a continuation process will involve using the commands b, i, b, i, b etc. of course you can insert other output commands or change variables as you go along. commands that define the newton iteration d = set the damping option. normally, newton tries to set xn new=xn old+correction. sometimes, xn new is not an improvement on xn old. in these cases, if damping is not used, then the program immediately returns to you. however, it is a fact that if 'damp' is a small enough number, xn new=xn old+damp*correction will be an improvement on xn old. (this assumes that we are using the standard newton method, and are not lagging the jacobian). if you specify damping, then the program will attempt to save the day by searching for a small enough damp to reduce the function norm of xn new. in fact, it sets damp to 1 until a failure occurs, and then repeatedly halves damp until the point is acceptable. we actually only allow 5 halvings before giving up. f = set the finite difference option. the program assumes that you have actually written formulas for the jacobian in your jacob subroutine. if you haven't, you can ask the program to approximate the jacobian using finite differences. even if you have, you can use this option to compare the results of exact and approximate jacobians. the finite difference options include backward, centered, and forward differences to approximate d f(i)/d x(j). k = define how often the jacobian is to be evaluated. although newton's method generally evaluates the jacobian at every step, it is cheaper (especially for large problems) to evaluate the jacobian only every so often. every so often is defined to be 1 for the original newton's method, but you can reset that to a higher value to see what happens. m = set maximum number of steps before return. to keep newton from running wild, the program allows it to take no more than m steps before it asks you for a new command, which might be simply 'c' for continue if you're happy with the progress of the iteration. right now, the number of steps is 10. change this value for your own taste, but don't make it too large. n = set the norm. the error tests for newton's method are based on the norms of the function, the point x and the step size (xn new-xn old). there are three norms available, the maximum norm, the sum of the absolute values, and the square root of the sum of squares. o = define how much output you want in newton iteration. 0 means that the newton iteration will only report whether it failed, succeeded, or hasn't reached a conclusion yet. 1 means that the norms of x, f(x) and the step will be printed for every step. 2 means that the actual vectors x, f(x) and the step will be printed for every step. t = set error tolerances. you can tighten or loosen the error tests this way. two numbers, abserr and relerr, are used. the error tests are of the following form where xnorm1, fnorm1 and snorm1 are the norms of the current x, fx and step, and xnorm0, fnorm0 and snorm0 are the same for the previous step. acceptance tests- step 1: accept if fnorm1.le.abserr step n: accept if fnorm1.le.abserr and snorm1.le.relerr*(xnorm1+abserr) rejection tests step 1: reject if fnorm1.gt.fnorm0+abserr step n: reject if fnorm1.gt.fnorm0+abserr or if snorm1.gt.snorm0+abserr execution commands b = begin newton iteration. you must give this command at least once before you can give the 'c' command. the newton iteration begins from the starting point you gave and takes up to the maximum number of steps. c = continue newton iteration. not to be given unless the b command has been given at least once. the newton iteration picks up from the last iterate of the iteration that had paused. other commands g = check jacobian versus finite difference jacobian. you should only issue this command when you are not in the middle of a newton iteration. the code will compare your jacobian with the finite difference jacobian. you can see the errors caused by discretization, and you may also be able to catch mistakes in your jacobian routine. l = list the commands. a short list of the commands is given, in case you forgot. v = print control quantities. the program will print the error tolerances, the damping option, the jacobian option, the norm used, the maximum number of newton steps allowed, and the number of newton steps taken before the jacobian is re-evaluated. if the continuation option is being used, you will also see the current stepsize and the current parameter. w = print work statistics. you can see how many times the function and jacobian were evaluated, how many times the jacobian was factored and solved. s = stop. *********************************************************************