alg125.txt

Numerical Anaysis 8th Edition Burden and Faires (Maple Source)

> printf(`Are the functions P,Q,R,F,G,G1,G2 ready for input or\n`);
> printf(`have they been included in code and has the input \n`);
> printf(`file been created?  Answer Y or N.\n`);
> AA := scanf(`\n%c`)[1];
> if AA = "Y" or AA = "y" then
> printf(`Input function P(X,Y) in terms of x and y\n`);
> P := scanf(`%a`)[1];
> P := unapply(P,x,y);
> printf(`Input function Q(X,Y) in terms of x and y\n`);
> Q := scanf(`%a`)[1];
> Q := unapply(Q,x,y);
> printf(`Input function R(X,Y) in terms of x and y\n`);
> R := scanf(`%a`)[1];
> R := unapply(R,x,y);
> printf(`Input function F(X,Y) in terms of x and y\n`);
> F := scanf(`%a`)[1];
> F := unapply(F,x,y);
> printf(`Input function G(X,Y) in terms of x and y\n`);
> G := scanf(`%a`)[1];
> G := unapply(G,x);
> printf(`Input function G1(X,Y) in terms of x and y\n`);
> G1 := scanf(`%a`)[1];
> G1 := unapply(G1,x);
> #printf(`Input function G2(X) in terms of x\n`);
> #G2 := scanf(`%a`)[1];
> #G2 := unapply(G2,x);
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`for example:  A:\\DATA.DTA\n`);
> NAME := scanf(`%s`)[1];
> 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 I from 1 to NL do
> for J from 1 to 2 do
> LINE[I-1,J-1] := fscanf(INP, `%d`)[1];
> od;
> od;
> fi;
> fclose(INP);
> else
> printf(`The program will end so that the input file \n`);
> printf(`can be created and the functions defined.\n`);
> fi;
> if OK = TRUE then
> printf(`Choice of output method:\n`);
> printf(`1. Output to screen\n`);
> printf(`2. Output to text file\n`);
> printf(`Please enter 1 or 2.\n`);
> FLAG := scanf(`%d`)[1];
> if FLAG = 2 then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`for example:  A:\\OUTPUT.DTA\n`);
> NAME := scanf(`%s`)[1];
> 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 I from 1 to M do
> fprintf(OUP, ` %5d`, I);
> for J from 1 to 3 do
> fprintf(OUP, ` %4d`, NTR[I-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, LINE[KK-1,0], LINE[KK-1,1]);
> od;
> # Step 1
> if LM <> LN1 then
> for L from N1 to LM do
> GAMMA[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 I from 1 to LN1 do
> for J from 1 to N1 do
> ALPHA[I-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 I from 1 to M do
> J1 := NTR[I-1,0];
> J2 := NTR[I-1,1];
> J3 := NTR[I-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[I-1,0] := (XX[J2-1]*YY[J3-1]-YY[J2-1]*XX[J3-1])/DELTA;
> B[I-1,0] := (YY[J2-1]-YY[J3-1])/DELTA;
> C[I-1,0] := (XX[J3-1]-XX[J2-1])/DELTA;
> A[I-1,1] := (XX[J3-1]*YY[J1-1]-YY[J3-1]*XX[J1-1])/DELTA;
> B[I-1,1] := (YY[J3-1]-YY[J1-1])/DELTA;
> C[I-1,1] := (XX[J1-1]-XX[J3-1])/DELTA;
> A[I-1,2] := (XX[J1-1]*YY[J2-1]-YY[J1-1]*XX[J2-1])/DELTA;
> B[I-1,2] := (YY[J1-1]-YY[J2-1])/DELTA;
> C[I-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[I-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[I-1,I2-1];
> ZZ := B[I-1,I1-1]*B[I-1,I2-1]*QQ(1,XX,YY,DELTA,J1,J2,J3,I1,I2,A,B,C)+C[I-1,I1-1]*C[I-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*GAMMA[JJ2-1];
> fi;
> else
> if JJ2 <= LN1 then
> ALPHA[JJ2-1,N1-1] := ALPHA[JJ2-1,N1-1]-ZZ*GAMMA[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 I from 1 to M do
> for J from 1 to 3 do
> fprintf(OUP, `%4d %3d %3d %13.8f %13.8f %13.8f\n`, I, J, NTR[I-1,J-1], 
> A[I-1,J-1],B[I-1,J-1], C[I-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 I 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[I-1,I1-1];
> if LINE[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[I-1,I2-1];
> if LINE[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,I,I1,I2,A,B,C);
> XJ1 := SQ(4,XX,YY,J1,J2,I,I1,I2,A,B,C);
> XJ2 := SQ(5,XX,YY,J1,J2,I,I1,I2,A,B,C);
> XI1 := SQ(2,XX,YY,J1,J2,I,I1,I2,A,B,C);
> XI2 := SQ(3,XX,YY,J1,J2,I,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*GAMMA[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*GAMMA[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 I from 1 to LN1 do
> fprintf(OUP, `Row %4d\n`, I);
> for J from 1 to N1 do
> fprintf(OUP, ` %13.10e\n `, evalf(ALPHA[I-1,J-1]));
> od;
> od;
> fprintf(OUP, `\n`);
> # Step 20
> if LN1 > 1 then
> INN := LN1-1;
> for I from 1 to INN do
> I1 := I+1;
> for J from I1 to LN1 do
> CCC := ALPHA[J-1,I-1]/ALPHA[I-1,I-1];
> for J1 from I1 to N1 do
> ALPHA[J-1,J1-1] := ALPHA[J-1,J1-1]-CCC*ALPHA[I-1,J1-1];
> od;
> ALPHA[J-1,I-1] := 0;
> od;
> od;
> fi;
> GAMMA[LN1-1] := ALPHA[LN1-1,N1-1]/ALPHA[LN1-1,LN1-1];
> if LN1 > 1 then
> for I from 1 to INN do
> J := LN1-I;
> CCC := ALPHA[J-1,N1-1];
> J1 := J+1;
> for KK from J1 to LN1 do
> CCC := CCC-ALPHA[J-1,KK-1]*GAMMA[KK-1];
> od;
> GAMMA[J-1] := CCC/ALPHA[J-1,J-1];
> od;
> fi;
> # Step 21
> # Output GAMMA
> fprintf(OUP, `Coefficients of Basis Functions:\n`);
> for I from 1 to LM do
> LLL := LL[I-1];
> fprintf(OUP, `%3d %13.8f %d`, I, evalf(GAMMA[I-1]), I);
> for J from 1 to LLL do
> fprintf(OUP, ` %4d`, II[I-1,J-1]);
> od;
> fprintf(OUP, `\n`);
> od;
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> # Step 23
> fi;
> RETURN(0);
> end;
> alg125();