CUBIC SPLINE
RAYLEIGH-RITZ ALGORITHM
*
* To approximate the solution to the boundary-value problem
*
* -D(P(X)Y')/DX + Q(X)Y = F(X), 0 <= X <= 1, Y(0)=Y(1)=0
*
* With a sum of cubic splines:
*
* INPUT: Integer n
*
* OUTPUT: Coefficients C(0),...,C(n+1) of the basis functions
*
* To change problems do the following:
* 1. Change functions P, Q and F in the subprograms named P,Q,F
* 2. Change dimensions in main program to values given by
* dimension A(N+2,N+3),X(N+2),C(N+2),CO(N+2,4,4),
* DCO(N+2,4,3),AF(N+1),BF(N+1),CF(N+1),DF(N+1),
* AP(N+1),BP(N+1),CP(N+1),DP(N+1),
* AQ(N+1),BQ(N+1),CQ(N+1),DQ(N+1)
* 3. Change dimensions in subroutine COEF to
* dimension AA(N+2),BB(N+2),CC(N+2),DD(N+2),
* H(N+2),XA(N+2),XL(N+2),XU(N+2),XZ(N+2)
*
* GENERAL OUTLINE
*
* 1. Nodes labelled X(I)=(I-1)*H, 1 <= I <= N+2, where
* H=1/(N+1) so that zero subscripts are avoided
* 2. The functions PHI(I) and PHI'(I) are shifted so that
* PHI(1) and PHI'(1) are centered at X(1), PHI(2) and PHI'(2)
* are centered at X(2), . . . , PHI(N+2) and
* PHI'(N+2) are centered at (X(N+2)---for example,
* PHI(3) = S((X-X(3))/H)
* = S(X/H + 2)
* 3. The functions PHI(I) are represented in terms of their
* coefficients in the following way:
* (PHI(I))(X) = CO(I,K,1) + CO(I,K,2)*X + CO(I,K,3)*X**2
* CO(I,K,4) *X**3
* for X(J) <= X <= X(J+1) where
* K=1 IF J=I-2, K=2 IF J=I-1, K=3 IF J=I, K=4 IF J=I+1
* since PHI(I) is nonzero only between X(I-2) and X(I+2)
* unless I = 1, 2, N+1 or N+2
* (see subroutine PHICO)
* 4. The derivative of PHI(I) denoted PHI'(I) is represented
* as in 3. By its coefficients DCO(I,K,L), L = 1, 2, 3
* (See subroutine DPHICO).
* 5. The functions P,Q and F are represented by their cubic
* spline interpolants using clamped boundary conditions
* (see Algorithm 3.5). Thus, for X(I) <= X <= X(I+1) we
* use AF(I)+BF(I)*(X-X[I])+CF(I)*(X-X[I])^2+DF(I)*(X-X[I])^3
* to represent F(X). Similarly, AP,BP,CP,DP are used for P
* and AQ,BQ,CQ,DQ are used for Q. (see subroutine COEF).
* 6. The integrands in STEPS 6 and 9 are replaced by products
* of cubic polynomial approximations on each subinterval of
* length H and the integrals of the resulting polynomials
* are computed exactly. (see subroutine XINT).
LINEAR
FINITE-DIFFERENCE ALGORITHM
*
* To approximate the solution of the boundary-value problem
*
* Y'' = P(X)Y' + Q(X)Y + R(X), A<=X<=B, Y(A) = ALPHA, Y(B) = BETA:
*
* INPUT: Endpoints A, B; boundary conditions ALPHA, BETA;
* integer N.
*
* OUTPUT: Approximations W(I) to Y(X(I)) for each I=0,1,...,N+1.
LINEAR SHOOTING ALGORITHM
*
* To approximate the solution of the boundary-value problem
*
* -Y'' + P(X)Y' + Q(X)Y + R(X) = 0, A<=X<=B, Y(A)=ALPHA, Y(B)=BETA:
*
*
* INPUT: Endpoints A,B; boundary conditions ALPHA, BETA; number of
* subintervals N.
*
* OUTPUT: Approximations W(1,I) to Y(X(I)); W(2,I) to Y'(X(I))
* for each I=0,1,...,N.
NONLINEAR
FINITE-DIFFERENCE ALGORITHM
*
* To approximate the solution to the nonlinear boundary-value problem
* Y'' = F(X,Y,Y'), A<=X<=B, Y(A) = ALPHA, Y(B) = BETA:
*
* INPUT: Endpoints A,B; boundary conditions ALPHA, BETA;
* integer N; tolerance TOL; maximum number of iterations M.
*
* OUTPUT: Approximations W(I) TO Y(X(I)) for each I=0,1,...,N+1
* or a message that the maximum number of iterations was
* exceeded.
NONLINEAR
SHOOTING ALGORITHM
*
* To approximate the solution of the nonlinear boundary-value problem
*
* Y'' = F(X,Y,Y'), A<=X<=B, Y(A) = ALPHA, Y(B) = BETA:
*
*
* INPUT: Endpoints A,B; boundary conditions ALPHA, BETA; number of
* subintervals N; tolerance TOL; maximum number of iterations M.
*
* OUTPUT: Approximations W(1,I) TO Y(X(I)); W(2,I) TO Y'(X(I))
* for each I=0,1,...,N or a message that the maximum
* number of iterations was exceeded.
PIECEWISE
LINEAR RAYLEIGH-RITZ ALGORITHM
*
* To approximate the solution of the boundary-value problem
*
* -D(P(X)Y')/DX + Q(X)Y = F(X), 0 <= X <= 1,
* Y(0) = Y(1) = 0,
*
* with a piecewise linear function:
*
* INPUT: integer N; mesh points X(0) = 0 < X(1) < ...
* < X(N) < X(N+1) = 1
*
* OUTPUT: coefficients C(1),...,C(N) of the basis functions.