In Section 5, we described the input and output routines. The permanent material also includes routines for the solution of linear equations, the solution of systems of ordinary differential equations and operations on matrices. Further routines may be added from time to time.
routine spec eqn solve(arrayname A,b, realname det)
This routine solves the equations
A(1,1)x(1) + A(1,2)x(2) +....+ A(1,n)x(n) = b(1) A(2,1)x(1) + A(2,2)x(2) +....+ A(2,n)x(n) = b(2) ' ' ' ' ' ' ' ' ' A(n,1)x(1) + A(n,2)x(2) +....+ A(n,m)x(n) = b(n)
(i.e. Ax = b), where the coefficients A(i,j) are stored in the matrix A, and b(i) in the vector b. A is destroyed and the solution is placed in b. If during the elimination process, the equations are found to be linearly dependant, then 'det' is set to zero and the routine is left, with both A and b upset, otherwise 'det' is set to the determinant of A. Consequently 'det' should be tested after each call of the routine.
The matrix routines operate on two dimensional arrays(i.e. matrices not vectors). The dimensions of the arrays are not required as parameter as the routines automatically find these from the declarations, and check them for compatibility. The programmer may insert similar tests in his own routines by means of the functions
integer fn spec dim (arrayname A) integer fn spec bound (arrayname A, integer n)
The first gives the dimensionality of the array(1 for a vector, 2 for a matrix etc.), and the second the nth bound (upper or lower) of the array counting from left to right. For example, if A were declared by:-
array A(-5:+5, 1:p)
where p = 10 then
dim(A) would have the value 2
bound(A,1) ' ' ' ' -5
bound(A,2) ' ' ' ' +5
bound(A,3) ' ' ' ' 1
bound(A,4) ' ' ' ' 10
The routines
routine spec unit (arrayname A) routine spec null (arrayname A)
set A to be a unit matrix (checking that it is square) and a null matrix respectively. The routines
routine spec matrix add (arrayname A,B,C) routine spec matrix sub (arrayname A,B,C) routine spec matrix copy (arrayname A,B)
set A to B+C, B-C and B respectively. Although the parameters are of type arrayname, the operation of the routines is such that the same array can be substituted for more than one of the parameters. For example:-
matrix add(A,A,A)
doubles A. The same is not true of the following routines:-
routine spec matrix mult(arrayname A,B,C) routine spec matrix mult'(arrayname A,B,C) routine spec matrix trans (arrayname A,B)
These set A to B*C, B*C' and B' respectively where the ' denotes transposition. If it is required to, say, set a matrix to the product of itself and another then the call
matrix mult(A,A,B)
will fail. It is necessary to declare another array, 'dummy' say, and then use 'matrix copy' and'matrix mult':-
matrix copy (dummy,A) matrix mult(A,dummy,B)
Alternatively, a routine with parameters of type array may be defined which calls the permanent routines :-
routine MATRIX MULT(arrayname A, array B,C)
matrix mult(A,B,C)
end
In this case a call of the form
MATRIX MULT(A,A,B) or even MATRIX MULT(A,A,A)
is possible.
The routines
routine spec matrix div(arrayname A,B, realname det) routine spec invert(array name A,B, realname det)
set A to inv(B) A and inv(B) respectively. In the process B is destroyed and the value of its determinant placed in det. Should the matrix be found to be singular, 'det' is set to zero. Consequently 'det' should be tested after every call for these routines. If B is required at the end of the routine, then the techniques described above should by used. The function
real fn det (arrayname B)
sets 'det' to the determinant of B and destroys B.
There are two routines available for advancing the solution of a system of first order ordinary differential equations
dy(i)/dx = f(i)(x,y(1),y(2),....,y(n)) i=1,2,...,n
from x to x+h, using the Kutta-Merson fourth-order integration method. The system is defined by means of an auxiliary routine, which must be supplied by the user, of the form:-
routine spec aux(arrayname f, real x)
which must evaluate the derivatives f(i) in terms of y(1),y(2),..,y(n) and x and then place them in f(1),f(2), ....., f(n).
The first routine
routine spec int step(arrayname y,real x,h, integer n c
realname e, routine aux)
advances the solution by a single step of length h of the Kutta-Merson process.
[2] Reference, L.FOX (Ed. ) Numerical Solution of Ordinary and Partial Differential equations. Pergamon 1962, P.24.
The parameters are:-
The second routine
routine spec kutta merson(array name y,real x0,x1, realname e c
integer n,k, routine aux)
advances the solution, by means of a series of calls for 'intstep', from x0 to x1, keeping, if possible the estimate of the maximum truncation error less than e. An initial step length of (x1-x0)/2↑m where 2↑(m+1) >k ≥ 2↑m, is taken. If over a step the local truncation error (given by 'int step') is greater than e, then the step length is halved; if the error is less than .01e then the step length is doubled. If three successive reductions in step length give no improvement in the estimated truncation error, then e is replaced by twice the smallest error achieved, and the integration process continued. The parameters are :-