> restart:
> # Finite Element Algorithm 12.5
> #
> # 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, 
> #       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) GAMMA1 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.
> global P,Q,R,F,G,G1,G2,G3,G4,G5,RR,FFF,GG1,GG2,GG3,GG4,GG5,QQ,SQ,S:
> # The following code is commented out.  It could be used to include
> # the functions P, Q, R, F, G, G1 as part of the code.
> #P := proc(X,Y) local p:
> #p := 1:
> #RETURN(p):
> #Q := proc(X,Y) local q:
> #q := 1:
> #RETURN(q):
> #R := proc(X,Y) local r:
> #r := 0:
> #RETURN(r):
> #F := proc(X,Y) local f:
> #f := 0:
> #RETURN(f):
> #G := proc(X,Y) local g:
> #g :=  4:
> #RETURN(g):
> #end:
> #G1 := proc(X,Y) local g1:
> #g1 :=  0:
> #RETURN(g1):
> #end:
> # Function G2 is defined in code for this example instead of input
> G2 := proc(X,Y) local T, Z1, g2:
> T := 1.0E-05:
> Z1 := 0:
> if 0.2-T <= X and X <= 0.4+T and abs(Y-0.2) <= T then
> Z1 := X:
> fi:
> if 0.5-T <= X and X <= (0.6+T) and abs(Y-0.1) <= T then
> Z1 := X:
> fi:
> if -T <= Y and Y <= 0.1+T and abs(X-0.6) <= T then
> Z1 := Y:
> fi:
> if -T <= X and X <= 0.2+T and abs(Y+X-0.4) <= T then
> Z1 := (X+Y)/sqrt(2):
> fi:
> if 0.4 -T <= X and X <= 0.5+T and abs(Y+X-0.6) <= T then
> Z1 := (X+Y)/sqrt(2):
> fi:
> g2 := Z1:
> RETURN(g2):
> end:
> RR := proc(X,Y,A,B,C,I3,I1,I2) local rr:
> rr := R(X,Y)*(A[I3-1,I1-1]+B[I3-1,I1-1]*X+C[I3-1,I1-1]*Y)*(A[I3-1,I2-1]+B[I3-1,I2-1]*X+C[I3-1,I2-1]*Y):
> RETURN(rr):
> end:
> FFF := proc(X,Y,A,B,C,I3,I1) local fff:
> fff := F(X,Y)*(A[I3-1,I1-1]+B[I3-1,I1-1]*X+C[I3-1,I1-1]*Y):
> RETURN(fff):
> end:
> GG1 := proc(X,Y,A,B,C,I3,I1,I2) local gg1:
> gg1 := G1(X,Y)*(A[I3-1,I1-1]+B[I3-1,I1-1]*X+C[I3-1,I1-1]*Y)*(A[I3-1,I2-1]+B[I3-1,I2-1]*X+C[I3-1,I2-1]*Y):
> RETURN(gg1):
> end:
> GG2 := proc(X,Y,A,B,C,I3,I1) local gg2:
> gg2 := G2(X,Y)*(A[I3-1,I1-1]+B[I3-1,I1-1]*X+C[I3-1,I1-1]*Y):
> RETURN(gg2):
> end:
> GG3 := proc(X,Y,A,B,C,I3,I2) local gg3:
> gg3 := G2(X,Y)*(A[I3-1,I2-1]+B[I3-1,I2-1]*X+C[I3-1,I2-1]*Y):
> RETURN(gg3):
> end:
> GG4 := proc(X,Y,A,B,C,I3,I1) local gg4:
> gg4 := G1(X,Y)*(A[I3-1,I1-1]+B[I3-1,I1-1]*X+C[I3-1,I1-1]*Y)*(A[I3-1,I1-1]+B[I3-1,I1-1]*X+C[I3-1,I1-1]*Y):
> RETURN(gg4):
> end:
> GG5 := proc(X,Y,A,B,C,I3,I2) local gg5:
> gg5 := G1(X,Y)*(A[I3-1,I2-1]+B[I3-1,I2-1]*X+C[I3-1,I2-1]*Y)*(A[I3-1,I2-1]+B[I3-1,I2-1]*X+C[I3-1,I2-1]*Y):
> RETURN(gg5):
> end:
> QQ := proc(L,XX,YY,DELTA,J1,J2,J3,I1,I2,A,B,C)  local X, Y, I3, T1, T2, T3, QQQ, qq: global S:
> X[0] := XX[J1-1]:
> Y[0] := YY[J1-1]:
> X[1] := XX[J2-1]:
> Y[1] := YY[J2-1]:
> X[2] := XX[J3-1]:
> Y[2] := YY[J3-1]:
> X[3] := 0.5*(X[0]+X[1]):
> Y[3] := 0.5*(Y[0]+Y[1]):
> X[4] := 0.5*(X[0]+X[2]):
> Y[4] := 0.5*(Y[0]+Y[2]):
> X[5] := 0.5*(X[1]+X[2]):
> Y[5] := 0.5*(Y[1]+Y[2]):
> X[6] := (X[0]+X[1]+X[2])/3.0:
> Y[6] := (Y[0]+Y[1]+Y[2])/3.0:
> if L = 1 then
> for I3 from 1 to 7 do
> S[I3-1] := P(X[I3-1],Y[I3-1]):
> od:    
> fi:
> if L = 2 then
> for I3 from 1 to 7 do
> S[I3-1] := Q(X[I3-1],Y[I3-1]):
> od:    
> fi:
> if L = 3 then
> for I3 from 1 to 7 do
> S[I3-1] := RR(X[I3-1],Y[I3-1],A,B,C,I3,I1,I2):
> od:    
> fi:
> if L = 4 then
> for I3 from 1 to 7 do
> S[I3-1] := FFF(X[I3-1],Y[I3-1],A,B,C,I3,I1):
> od:    
> fi:
> T1 := 0:
> for I3 from 1 to 3 do
> T1 := T1+S[I3-1]:
> od:
> T2 := 0:
> for I3 from 4 to 6 do
> T2 := T2+S[I3-1]:
> od:
> T3 := S[6]:
> QQQ := 0.5*(T1/20+2*T2/15+9*T3/20)*abs(DELTA):
> qq := QQQ:
> RETURN(qq):
> end:
> SQ := proc(L,XX,YY,J1,J2,I3,I1,I2,A,B,C)  local X1, Y1, X2, Y2, H, T1, T2, T3, K, X, Q3, Q1, Q2, SSQ, sq: global S:
> X1 := XX[J1-1]:
> Y1 := YY[J1-1]:
> X2 := XX[J2-1]:
> Y2 := YY[J2-1]:
> H := 0.01:
> T1 := X2-X1:
> T2 := Y2-Y1:
> T3 := sqrt(T1*T1+T2*T2):
> for K from 1 to 101 do
> X := (K-1)*H:
> if L = 1 then
> S[K-1] := T3*GG1(T1*X+X1,T2*X+Y1,A,B,C,I3,I1,I2):
> fi:
> if L = 2 then
> S[K-1] := T3*GG2(T1*X+X1,T2*X+Y1,A,B,C,I3,I1):
> fi:
> if L = 3 then
> S[K-1] := T3*GG3(T1*X+X1,T2*X+Y1,A,B,C,I3,I2):
> fi:
> if L = 4 then
> S[K-1] := T3*GG4(T1*X+X1,T2*X+Y1,A,B,C,I3,I1):
> fi:
> if L = 5 then
> S[K-1] := T3*GG5(T1*X+X1,T2*X+Y1,A,B,C,I3,I2):
> fi:
> od:
> Q3 := S[0]+S[100]:
> Q1 := 0:
> Q2 := S[99]:
> for K from 1 to 49 do
> Q1 := Q1+S[2*K]:
> Q2 := Q2+S[2*K-1]:
> od:
> SSQ := (Q3+2*Q1+4*Q2)*H/3:
> sq := SSQ:
> RETURN(sq):
> end:
> print(`This is the Finite Element Method.`):
> OK := FALSE:
> print(`The input will be from a text file in the following form:`):
> print(`K  N  M  n  m  nl`):
> print(`Each of the above is an integer -`):
> print(`separate with at least one blank.`):
> print(`Follow with the input for each node in the form:`):
> print(`x-coord., y-coord., number of triangles in which the`):
> print(` node is a vertex.`):
> print(`Continue with the triangle number and vertex number for`):
> print(`each triangle in which the node is a vertex.`):
> print(`Separate each entry with at least one blank.`):
> print(`After all nodes have been entered follow with information`):
> print(`on the lines over which line integrals must be computed.`):
> print(`The format of this data will be the node number of the`):
> print(`starting node, followed by the node number of the ending`):
> print(`node for each line, taken in the positive direction.`):
> print(`There should be 2 * nl such entries, each an integer`):
> print(`separated by a blank.`):
> print(`Functions can be input or coded as procedures.`):
> print(`The example has G2 as a procedure and the rest`):
> print(`as input.`):
> print(`Are the functions P,Q,R,F,G,G1,G2 ready for input or`):
> print(`have they been included in code and has the input `):
> print(`file been created?  Answer 1 for yes or 2 for no.`):
> AA := scanf(`%d`)[1]:  print(`Input is `): print(AA):
> if AA = 1 then
> print(`Input function P(X,Y) in terms of x and y`):
> P := scanf(`%a`)[1]: print(`P(x,y) = `): print(P):
> P := unapply(P,x,y):
> print(`Input function Q(X,Y) in terms of x and y`):
> Q := scanf(`%a`)[1]: print(`Q(x,y) = `): print(Q):
> Q := unapply(Q,x,y):
> print(`Input function R(X,Y) in terms of x and y`):
> R := scanf(`%a`)[1]: print(`R(x,y) = `): print(R):
> R := unapply(R,x,y):
> print(`Input function F(X,Y) in terms of x and y`):
> F := scanf(`%a`)[1]: print(`F(x,y) = `): print(F):
> F := unapply(F,x,y):
> print(`Input function G(X,Y) in terms of x and y`):
> G := scanf(`%a`)[1]: print(`G(x,y) = `): print(G):
> G := unapply(G,x):
> print(`Input function G1(X,Y) in terms of x and y`):
> G1 := scanf(`%a`)[1]: print(`G1(x,y) = `): print(G1):
> G1 := unapply(G1,x):
> #print(`Input function G2(X) in terms of x`):
> #G2 := scanf(`%a`)[1]: print(`G2(x,y) = `): print(G2):
> #G2 := unapply(G2,x):
> print(`Input the file name in the form - drive:\\name.ext`):
> print(`for example:  A:\\DATA.DTA`):
> NAME := scanf(`%s`)[1]:  print(`Input file is `): print(NAME):
> INP := fopen(NAME,READ,TEXT):
> OK := TRUE:
> K := fscanf(INP, `%d`)[1]:
> N := fscanf(INP, `%d`)[1]:
> M := fscanf(INP, `%d`)[1]:
> LN1 := fscanf(INP, `%d`)[1]:
> LM := fscanf(INP, `%d`)[1]:
> NL := fscanf(INP, `%d`)[1]:
> for KK from 1 to LM do
> XX[KK-1] := fscanf(INP, `%f`)[1]:
> YY[KK-1] := fscanf(INP, `%f`)[1]:
> LLL := fscanf(INP, `%d`)[1]:
> for J from 1 to LLL do
> II[KK-1,J-1] := fscanf(INP, `%d`)[1]:
> NV[KK-1,J-1] := fscanf(INP, `%d`)[1]:
> od:
> LL[KK-1] := LLL:
> for J from 1 to LLL do
> N1 := II[KK-1,J-1]:
> N2 := NV[KK-1,J-1]:
> NTR[N1-1,N2-1] := KK:
> od:
> od:
> if NL > 0 then
> for I3 from 1 to NL do
> for J from 1 to 2 do
> LINE1[I3-1,J-1] := fscanf(INP, `%d`)[1]:
> od:
> od:
> fi:
> fclose(INP):
> else
> print(`The program will end so that the input file `):
> print(`can be created and the functions defined.`):
> fi:
> if OK = TRUE then
> print(`Choice of output method:`):
> print(`1. Output to screen`):
> print(`2. Output to text file`):
> print(`Please enter 1 or 2.`):
> FLAG := scanf(`%d`)[1]: print(`Input is `): print(FLAG):
> if FLAG = 2 then
> print(`Input the file name in the form - drive:\\name.ext`):
> print(`for example:  A:\\OUTPUT.DTA`):
> NAME := scanf(`%s`)[1]: print(`Output file is `): print(NAME):
> OUP := fopen(NAME,WRITE,TEXT):
> else
> OUP := default:
> fi:
> fprintf(OUP, `FINITE ELEMENT METHOD\n\n\n`):
> K1 := K+1:
> N1 := LN1+1:
> fprintf(OUP, `Vertices and Nodes of Triangles\n`):
> fprintf(OUP, `Triangle-node number for vertex 1 to 3\n`):
> for I3 from 1 to M do
> fprintf(OUP, ` %5d`, I3):
> for J from 1 to 3 do
> fprintf(OUP, ` %4d`, NTR[I3-1,J-1]):
> od:
> fprintf(OUP, `\n`):
> od:
> fprintf(OUP, `x and y coordinates of nodes\n`):
> for KK from 1 to LM do
> fprintf(OUP, `%5d %11.8f %11.8f\n`, KK, XX[KK-1], YY[KK-1]):
> od:
> fprintf(OUP, `Lines of the Domain\n`):
> for KK from 1 to NL do
> fprintf(OUP, `%5d %4d %4d\n`, KK, LINE1[KK-1,0], LINE1[KK-1,1]):
> od:
> # Step 1
> if LM <> LN1 then
> for L from N1 to LM do
> GAMMA1[L-1] := G(XX[L-1],YY[L-1]):
> od:
> fi:
> # Step 2 - initialization of ALPHA
> # Note that ALHPA[I,LN1+1] = BETA[I]
> for I3 from 1 to LN1 do
> for J from 1 to N1 do
> ALPHA[I3-1,J-1] := 0:
> od:
> # Steps 3, 4 and 6 - 12 are within the next loop
> # For each triangle I let node J1 be vertex 1, node J2 be vertex 2,
> # and node J3 be vertex 3
> # Step 3
> od:
> for I3 from 1 to M do
> J1 := NTR[I3-1,0]:
> J2 := NTR[I3-1,1]:
> J3 := NTR[I3-1,2]:
> DELTA := 
> XX[J2-1]*YY[J3-1]-XX[J3-1]*YY[J2-1]-XX[J1-1]*(YY[J3-1]-YY[J2-1])+YY[J1-1]*
> (XX[J3-1]-XX[J2-1]):
> A[I3-1,0] := (XX[J2-1]*YY[J3-1]-YY[J2-1]*XX[J3-1])/DELTA:
> B[I3-1,0] := (YY[J2-1]-YY[J3-1])/DELTA:
> C[I3-1,0] := (XX[J3-1]-XX[J2-1])/DELTA:
> A[I3-1,1] := (XX[J3-1]*YY[J1-1]-YY[J3-1]*XX[J1-1])/DELTA:
> B[I3-1,1] := (YY[J3-1]-YY[J1-1])/DELTA:
> C[I3-1,1] := (XX[J1-1]-XX[J3-1])/DELTA:
> A[I3-1,2] := (XX[J1-1]*YY[J2-1]-YY[J1-1]*XX[J2-1])/DELTA:
> B[I3-1,2] := (YY[J1-1]-YY[J2-1])/DELTA:
> C[I3-1,2] := (XX[J2-1]-XX[J1-1])/DELTA:
> # Step 4
> # I1 = J for Step 4 and I1 = K for Step 7
> for I1 from 1 to 3 do
> # Step 8
> JJ1 := NTR[I3-1,I1-1]:
> # I2 = K for Step 4 and I2 = J for Step 9
> for I2 from 1 to I1 do
> # Steps 4 and 10
> JJ2 := NTR[I3-1,I2-1]:
> ZZ := B[I3-1,I1-1]*B[I3-1,I2-1]*QQ(1,XX,YY,DELTA,J1,J2,J3,I1,I2,A,B,C)+C[I3-1,I1-1]*C[I3-1,I2-1]*QQ(2,XX,YY,DELTA,J1,J2,J3,I1,I2,A,B,C)-QQ(3,XX,YY,DELTA,J1,J2,J3,I1,I2,A,B,C):
> # Steps 11 and 12
> if JJ1 <= LN1 then
> if JJ2 <= LN1 then
> ALPHA[JJ1-1,JJ2-1] := ALPHA[JJ1-1,JJ2-1]+ZZ:
> if JJ1 <> JJ2 then
> ALPHA[JJ2-1,JJ1-1] := ALPHA[JJ2-1,JJ1-1]+ZZ:
> fi:
> else
> ALPHA[JJ1-1,N1-1] := ALPHA[JJ1-1,N1-1]-ZZ*GAMMA1[JJ2-1]:
> fi:
> else
> if JJ2 <= LN1 then
> ALPHA[JJ2-1,N1-1] := ALPHA[JJ2-1,N1-1]-ZZ*GAMMA1[JJ1-1]:
> fi:
> fi:
> od:
> HH := -QQ(4,XX,YY,DELTA,J1,J2,J3,I1,I2,A,B,C):
> if JJ1 <= LN1 then
> ALPHA[JJ1-1,N1-1] := ALPHA[JJ1-1,N1-1]+HH:
> fi:
> od:
> od:
> # Output the basis functions
> fprintf(OUP, `Basis Functions\n`):
> fprintf(OUP, `Triangle - Vertex - Node - Function\n`):
> for I3 from 1 to M do
> for J from 1 to 3 do
> fprintf(OUP, `%4d %3d %3d %13.8f %13.8f %13.8f\n`, I3, J, NTR[I3-1,J-1], 
> A[I3-1,J-1],B[I3-1,J-1], C[I3-1,J-1]):
> od:
> od:
> # Step 5
> # For each line segment JI = 1, ..., NL and for each triangle I,
> # I = K1, ..., N which may contain line JI search all 3
> # vertices for possible correspondences.
> # Step 5 and Steps 13 - 19
> if NL <> 0 and N <> K then
> for JI from 1 to NL do
> for I3 from K1 to N do
> for I1 from 1 to 3 do
> # I1 = J in Step 5 and I1 = K in Step 14
> # Step 15
> J1 := NTR[I3-1,I1-1]:
> if LINE1[JI-1,0] = J1 then
> for I2 from 1 to 3 do
> # I2 = K in Step 5 and I2 = J in Step 16
> # Step 17
> J2 := NTR[I3-1,I2-1]:
> if LINE1[JI-1,1] = J2 then
> # There is a correspondence of vertix I1 in triangle I with node J1
> # as the start of line JI and vertex I2 with node J2 as the
> # end of line JI
> # Step 5
> XJ := SQ(1,XX,YY,J1,J2,I3,I1,I2,A,B,C):
> XJ1 := SQ(4,XX,YY,J1,J2,I3,I1,I2,A,B,C):
> XJ2 := SQ(5,XX,YY,J1,J2,I3,I1,I2,A,B,C):
> XI1 := SQ(2,XX,YY,J1,J2,I3,I1,I2,A,B,C):
> XI2 := SQ(3,XX,YY,J1,J2,I3,I1,I2,A,B,C):
> # Steps 8 and 19
> if J1 <= LN1 then
> if J2 <= LN1 then
> ALPHA[J1-1,J1-1] := ALPHA[J1-1,J1-1]+XJ1:
> ALPHA[J1-1,J2-1] := ALPHA[J1-1,J2-1]+XJ:
> ALPHA[J2-1,J2-1] := ALPHA[J2-1,J2-1]+XJ2:
> ALPHA[J2-1,J1-1] := ALPHA[J2-1,J1-1]+XJ:
> ALPHA[J1-1,N1-1] := ALPHA[J1-1,N1-1]+XI1:
> ALPHA[J2-1,N1-1] := ALPHA[J2-1,N1-1]+XI2:
> else
> ALPHA[J1-1,N1-1] := ALPHA[J1-1,N1-1]-XJ*GAMMA1[J2-1]+XI1:
> ALPHA[J1-1,J1-1] := ALPHA[J1-1,J1-1]+XJ1:
> fi:
> else
> if J2 <= LN1 then
> ALPHA[J2-1,N1-1] := ALPHA[J2-1,N1-1]-XJ*GAMMA1[J1-1]+XI2:
> ALPHA[J2-1,J2-1] := ALPHA[J2-1,J2-1]+XJ2:
> fi:
> fi:
> fi:
> od:
> fi:
> od:
> od:
> od:
> fi:
> # Output ALPHA
> fprintf(OUP, `Matrix ALPHA follows:\n`):
> for I3 from 1 to LN1 do
> fprintf(OUP, `Row %4d\n`, I3):
> for J from 1 to N1 do
> fprintf(OUP, ` %13.10e\n `, evalf(ALPHA[I3-1,J-1])):
> od:
> od:
> fprintf(OUP, `\n`):
> # Step 20
> if LN1 > 1 then
> INN := LN1-1:
> for I3 from 1 to INN do
> I1 := I3+1:
> for J from I1 to LN1 do
> CCC := ALPHA[J-1,I3-1]/ALPHA[I3-1,I3-1]:
> for J1 from I1 to N1 do
> ALPHA[J-1,J1-1] := ALPHA[J-1,J1-1]-CCC*ALPHA[I3-1,J1-1]:
> od:
> ALPHA[J-1,I3-1] := 0:
> od:
> od:
> fi:
> GAMMA1[LN1-1] := ALPHA[LN1-1,N1-1]/ALPHA[LN1-1,LN1-1]:
> if LN1 > 1 then
> for I3 from 1 to INN do
> J := LN1-I3:
> CCC := ALPHA[J-1,N1-1]:
> J1 := J+1:
> for KK from J1 to LN1 do
> CCC := CCC-ALPHA[J-1,KK-1]*GAMMA1[KK-1]:
> od:
> GAMMA1[J-1] := CCC/ALPHA[J-1,J-1]:
> od:
> fi:
> # Step 21
> # Output GAMMA1
> fprintf(OUP, `Coefficients of Basis Functions:\n`):
> for I3 from 1 to LM do
> LLL := LL[I3-1]:
> fprintf(OUP, `%3d %13.8f %d`, I3, evalf(GAMMA1[I3-1]), I3):
> for J from 1 to LLL do
> fprintf(OUP, ` %4d`, II[I3-1,J-1]):
> od:
> fprintf(OUP, `\n`):
> od:
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,` created successfully`):
> fi:
> # Step 23
> fi: