CRANK-NICOLSON ALGORITHM
*
* To approximate the solution of the parabolic partial-differential
* equation subject to the boundary conditions
* u(0,t) = u(l,t) = 0, 0 < t < T = max t
* and the initial conditions
* u(x,0) = F(x), 0 <= x <= l:
*
* INPUT: endpoint l; maximum time T; constant ALPHA; integers m, N:
*
* OUTPUT: approximations W(I,J) to u(x(I),t(J)) for each
* I = 1,..., m-1 and J = 1,..., N.
Finite Element Algorithm
*
* To approximate the solution to an elliptic partial-differential
* equation subject to Dirichlet, mixed, or Neumann boundary
* conditions:
*
* Input: see STEP 0
*
* Output: description of triangles, nodes, line integrals, basis
* functions, linear system to be solved, and the
* coefficients of the basis functions
*
*
* Step 0
* General outline
*
* 1. Triangles numbered: 1 to K for triangles with no edges on
* Script-S-1 or Script-S-2, K+1 to N for triangles with
* edges on script-S-2, N+1 to M for remaining triangles.
* Note: K=0 implies that no triangle is interior to D.
* Note: M=N implies that all triangles have edges on
* script-S-2.
*
* 2. Nodes numbered: 1 to LN for interior nodes and nodes on
* script-S-2, LN+1 to LM for nodes on script-S-1.
* Note: LM and LN represent lower case m and n resp.
* Note: LN=LM implies that script-S-1 contains no nodes.
* Note: If a node is on both script-S-1 and script-S-2, then
* it should be treated as if it were only on script-S-1.
*
* 3. NL=number of line segments on script-S-2
* line(I,J) is an NL by 2 array where
* line(I,1)= first node on line I and
* line(I,2)= second node on line I taken
* in positive direction along line I
*
* 4. For the node labelled KK,KK=1,...,LM we have:
* A) coordinates XX(KK),YY(KK)
* B) number of triangles in which KK is a vertex= LL(KK)
* C) II(KK,J) labels the triangles KK is in and
* NV(KK,J) labels which vertex node KK is for
* each J=1,...,LL(KK)
*
* 5. NTR is an M by 3 array where
* NTR(I,1)=node number of vertex 1 in triangle I
* NTR(I,2)=node number of vertex 2 in triangle I
* NTR(I,3)=node number of vertex 3 in triangle I
*
* 6. Function subprograms:
* A) P,Q,R,F,G,G1,G2 are the functions specified by
* the particular differential equation
* B) RR is the integrand
* R*N(J)*N(K) on triangle I in step 4
* C) FFF is the integrand F*N(J) on triangle I in step 4
* D) GG1=G1*N(J)*N(K)
* GG2=G2*N(J)
* GG3=G2*N(K)
* GG4=G1*N(J)*N(J)
* GG5=G1*N(K)*N(K)
* integrands in step 5
* E) QQ(FF) computes the double integral of any
* integrand FF over a triangle with vertices given by
* nodes J1,J2,J3 - the method is an O(H**2) approximation
* for triangles
* F) SQ(PP) computes the line integral of any integrand PP
* along the line from (XX(J1),YY(J1)) to (XX(J2),YY(J2))
* by using a parameterization given by:
* X=XX(J1)+(XX(J2)-XX(J1))*T
* Y=YY(J1)+(YY(J2)-YY(J1))*T
* for 0 <= t <= 1
* and applying Simpson's composite method with H=.01
*
* 7. Arrays:
* A) A,B,C are M by 3 arrays where the basis function N
* for the Ith triangle, Jth vertex is
* N(X,Y)=A(I,J)+B(I,J)*X+C(I,J)*Y
* for J=1,2,3 and i=1,2,...,M
* B) XX,YY are LM by 1 arrays to hold coordinates of nodes
* C) line,LL,II,NV,NTR have been explained above
* D) Gamma and Alpha are clear
*
* 8. Note that A,B,C,XX,YY,I,I1,I2,J1,J2,J3,Delta are reserved
* storage so that in any subprogram we know that
* triangle I has vertices (XX(J1),YY(J1)),(XX(J2),YY(J2)),
* (XX(J3),YY(J3)). That is, vertex 1 is node J1, vertex 2 is
* node J2, vertex 3 is node J3 unless a line integral is
* involved in which case the line integral goes from node J1
* to node J2 in triangle I or unless vertex I1 is node J1
* and vertex I2 is node J2 - the basis functions involve
* A(I,I1)+B(I,I1)*X+C(I,I1)**Y for vertex I1 in triangle I
* and A(I,I2)+B(I,I2)*X+C(I,I2)*Y for vertex I2 in triangle I
* delta is 1/2 the area of triangle I
*
* To change problems:
* 1) change function subprograms P,Q,R,F,G,G1,G2
* 2) change data input for K,N,M,LN,LM,NL.
* 3) change data input for XX,YY,LLL,II,NV.
* 4) change data input for line.
* 5) change definition statements to read :
* A(M,3),B(M,3),C(M,3),XX(LM),YY(LM)
* definition LINE(NL,2),LL(LM),II(LM,MAX LL(LM)),
* NV(LM,MAX LL(LM))
* definition NTR(M,3),GAMMA(LM),ALPHA(LN,LN+1), use
* the appropriate numbers for the variables M, LM,
* NL.
HEAT
EQUATION BACKWARD-DIFFERENCE ALGORITHM
* To approximate the solution to the parabolic partial-differential
* equation subject to the boundary conditions
* u(0,t) = u(l,t) = 0, 0 < t < T = max t,
* and the initial conditions
* u(x,0) = F(x), 0 <= x <= l:
*
* INPUT: endpoint l; maximum time T; constant ALPHA; integers m, N.
*
* OUTPUT: approximations W(I,J) to u(x(I),t(J)) for each
* I = 1, ..., m-1 and J = 1, ..., N.
POISSON
EQUATION FINITE-DIFFERENCE ALGORITHM
* To approximate the solution to the Poisson equation
* DEL(u) = F(x,y), a <= x <= b, c <= y <= d,
* SUBJECT TO BOUNDARY CONDITIONS:
* u(x,y) = G(x,y),
* if x = a or x = b for c <= y <= d,
* if y = c or y = d for a <= x <= b
*
* INPUT: endpoints a, b, c, d; integers m, n; tolerance TOL;
* maximum number of iterations M
*
* OUTPUT: approximations W(I,J) to u(X(I),Y(J)) for each
* I = 1,..., n-1 and J=1,..., m-1 or a message that the
* maximum number of iterations was exceeded.
WAVE EQUATION
FINITE-DIFFERENCE ALGORITHM
*
* To approximate the solution to the wave equation:
* subject to the boundary conditions
* u(0,t) = u(l,t) = 0, 0 < t < T = max t
* and the initial conditions
* u(x,0) = F(x) and Du(x,0)/Dt = G(x), 0 <= x <= l:
*
* INPUT: endpoint l; maximum time T; constant ALPHA; integers m, N.
*
* OUTPUT: approximations W(I,J) to u(x(I),t(J)) for each I = 0, ..., m
* and J=0,...,N.