COMMENT THESE ROUTINES WERE SCANNED IN 1998 FROM LASER 9700 OUTPUT;
COMMENT
COPYRIGHT (C) 1981
COLUMBIA UNIVERSITY BIOENGINEERING INSTITUTE
&
VASOS-PETER JOHN PANAGIOTOPOULOS II
TECHNICAL STANFORD-SAIL-ALGOL SUBROUTINE PACKAGE ONE SBRPC1
This project was funded in part (since April, 1979) by
NSF Rat Skull Grant-Professor Melvin L. Moss-Fall 1979 only
NIH HD14371 & HL16851 Professor Richard Skalak
CU-CCA NCA-ID fund and other funds for educational computing
Samani International Enterprises Siengg Computational Research Fund
Vasos-Feter John Panagiotopoulos II
I wish to acknowledge the unlimited and most kind assistance of
Professor Gautam Dasgupta.
References utilised:
Numerical Analysis, Scheid, Schaum/McGrawHill
Applied Numerical Methods, Carnahan, Luther&Wilkes,Wiley
Introductory Computer Methods and Numerical Analysis,Pennington
Advanced Calculus for Applications,Hildebrand
Handbook of Mathematical Functions,Abramowitz&Stegun
The following Columbia courses:
ChE3100X79,ChE4510Y81,ChE4520YS1
ENTRY Factorial,Lagrange!lnterpolation,Average,Norm,Quadratic,matsum,
Summation!Of!Function,secant!root,Newton!Raphson,Regula!Falsi,
Power!Method,RUNGEK,realrd,read!arai,File!Read!Arai,Matrix!Print,
FIBONACCI,gauss!jordan,Simpson,Recirculation,Max!Array,Min!Array,
InstaPlot,cc!trapper,sink,Structural!Section,Monte!Carlo,
Binary!Monte!Carlo,Gradient!Search,Gamma,Legendre!P,
Weighted!Gauss!Legendre!Integration,
C!Conjugate,C!Add,C!Subtract,C!Abs,C!Multiply,C!Divide,C!ExP,CvCmplx
;
BEGIN "SBRPCl"
Require "sai:abbrev.sai" source!file;comment sail abbreviations package.;
Require "omni:disprm.sai" source!file;comment sail omnigraph package.;
Require "sai:routin.hdr" source!file;comment sail routine package.;
Require "sai:ttyio.hdr" source!file;comment sail i/o package.;
Define Tan(x)= {sin(x)/cos(x)},
Degrees (radians)={radians*180/4/atan(1)},
Round!Off(x)= {2*(x%2)},
Even!Odd(x)= {2*abs(x/2-x%2)},
Kroenecker(i,j)= {if i=j then 1 else 0},
Heaviside(x)= {if x>O then 1 else 0},
Dirac(x)= {if x=O then 1 else 0},
Permutation(i.J,k)={if i=j or j=k or i=k then 0 else
if (i=1 and j=2 and k=3) or(i=2 and j=3 and k=1)
or (i=3 and j=l and k=2)
then 1 else -1},
ARCOS(x)={4*ATAN(1)+ATAN((1+x^2)^.5/X)},
ARSIN(x)={-ARCOS((1-x^2)^.5)},
INFINITY={.1701412@39},
COMPLEX= {RECORD!POINTER(CMPLX)},
NEW!COMPLEX= {NEW!RECORD(CMLX)};
RECORD!CLASS CMPLX (REAL Re.lm);
Internal Real Procedure Factorial (Integer x);
BEGIN "Factorial Recursive (self-calling) procedure"
return(if x=O then I else x*factorial (x-1));
END "Factorial Recursive (self-calling) procedure";
Internal Real procedure Lagrange!Interpolation (Real Array x,y;
Integer n;
Real xx);
begin"Lagrange!Interpolation"
real array 1 [1:n]
integer i,j;
real yy;
for i:=1 step 1 until n do
begin "Lagrange looping i=/=j"
for j:=1 step 1 until n do if (j neq i) then
l[i] =l[i]*(xx-x(j])/(x[i]-x[j])
yy:=y[i]*l[i]+yy;
end "Lagrange looping i=/=j";
return (yy);
end"Lagrange!Interpolation";
Internal Real Procedure Integral!of!Polynomial!Term(real a,b,n,u1,u2);
begin "Integral Of Polynomial Term"
comment Integral of (A+BU)+N from U1 to U2;
return(if n neq -1 then((a+b*u2)f(n+1))/(n+1)/b-((a+b*u1)+(n+1))/(n+1)/b
else 1/b*(log(a+b*u2)-log(a+b*u1));
end "lntegral Of Polynomial Term";
Internal Real Procedure Average(real array a:integer n);
begin "avg"
real sum; integer m;m:=n;
while (n:=n-1) do sum:=sum+a[n];
return(sum/m);
end "avg";
Internal Real Procedure Norm (Real Array Vector;Integer Rank);
BEGIN "Norm of vector "
REAL Scalar!Norm;
while (rank:=rank-1) do Scalar!Norm:=Scalar!Norm+Vector[rank]+2;
return(Scalar!Norm^.5);
END "Norm of vector ";
Internal Procedure Quadratic(real a,b,c;
reference real first!root, second!root;
reference string root!type);
begin "Quadratic"
real g,d;
g:=b^2-4*a*c;
if (a=O) then
begin "real1"
first!root:=second!root:=-c/b;
root!type:="single real root";
end "real1";
else
begin "two-root"
if (g>0) then
begin "real2"
d:=sqrt(g);
first!root:=(-b-d)/(2.*a);
second!root:=( b-d)/(2.*a);
root!type:="double real roots";
end "real2"
else
begin "cmp2"
first!root:=-b/2./a;
second!root:=sqrt(-g)/2/a;
root!type:="imaginary first +/-second*i";
end "cmp2";
end "two-root";
return;
end "Quadratic";
Internal Real procedure matsum (real array x; integer dim);
begin"matsum" real ssum;
for dim:=dim step -l until 1 do ssum:=ssum+x [dim] ;return(ssum)
end"matsum";
Internal Real procedure
Summation!Of!Function(INTEGER Start!Index,End!Index,Interval;
REAL X;
REAL PROCEDURE Funct);
begin"Summation!Of!Function"
integer icount;real g;
for icount:=start!index step Interval until End!lndex do
g:=g+funct(x,icount);
return (g)
end"Summation!Of!Function";
Internal Real Procedure secant!root(REAL X0,X1,Accuracy;REAL PROCEDURE F);
begin"sec" real x2;
while (abs(f(xO))geq accuracy ) do
x0:=(f(x1:=x0)*(x2:=x1)*f(x2)*x1)
/(f(x1)-f(x2));
return (xO)
end"sec";
Internal Real procedure Newton!Raphson(REAL X0,X1,Accuracy;REAL PROCEDURE F);
Begin "Newton Raphson Root Determination Reverse diff'
Real x2;
While (abs(f(x2)-f(x1)) geq accuracy) do
x0:=((x2:=x1)-(x1:=x0))*f(x1)/(f(x1)-f(x2))+x1;
Return (xO)
End"Newton Raphson Root Determination Reverse diff";
Internal Real procedure Regula!Falsi(REAL XO,Xl,Accuracy;REAL PROCEDURE F):
Begin"Regula Falsi According to Firouztale, Spencer & Wright"
Real x2;
While (abs(f(x2)-f(x1)) geq accuracy) do
x0:=((x2:=x1)-(x1:=x0))*f(x1)/(f(x1)-f(x2))+x1;
Return (x0);
End"Regula Falsi According to Firouztale, Spencer & Wright";
lnternai Procedure Power!Method( Real array a;
Reference Real Array Eigen!Vector;
Integer Rank;
Reference Real Eigen!Value;
Real Accuracy);
BEGIN "POWERM: Power Method of Mises for Eigen Systems"
Comment Panagiotopoulos Strikes Again
Real Array Temp!Eigen!Pseudo[1:Rank];
Integer I,J;
Real Old!Eigen!Value;
While ( Eigen!Value-Old!Eigen!Value>Accuracy) do
BEGIN "Power Method Inner Block"
Old!Eigen!Value:=Eigen!Value;
Eigen!Value:=0;
For I:=1 step 1 until Rank do
BEGIN "Mises"
For J:=1 step 1 until Rank do
Temp!Eigen!Pseudo[i]:=Temp!Eigen!Pseudo[i]
+a[j,3]*Eigen!Vector[3];
Eigen!value:=Eigen!Value+Temp!Eigen!Pseudo[i]^2;
END "Mises";
Eigen!Value:=Eigen!Value^.5;
For I:=1 step 1 until Rank do
Eigen!Vector[i]:=Temp!Eigen!Pseudo[i]/Eigen!Value:
END "Power Method Inner Block";
If (Eigen!Vector[1) neq 0) then
While(Rank:=Rank-1) do
Eigen!Vector[Rank]:=Eigen!Vector[Rank]/Eigen!Vector[1];
Return;
END "POWERM: Power Method of Mises for Eigen Systems";
Internal Procedure RUNGEK(REFERENCE REAL ARRAY X;
INTEGER N!Equations;
REAL PROCEDURE DRV;
REAL Delta!X, End!x);
BEGIN "RUNGEK: Runge-Kutta-Gill fourth-order Integration Program"
Integer i,j,k;
Real XX;
Comment : 1)Derived from Taylor expansion of ODE
2)DRV must return derivative of Jth function
3)X has independant variable&'solutions
4)IC's must be entered as x[O, _];
for xx:=x[O,C] step DeltaIX until EndIX do
BEGIN "Integration Stepping"
Real Array C[1:4,1:N!Equations];
i:=i+1;x[i,0]:=xx;
for j:=1 step 1 until N!Equations do
BEGIN "Gill constants"
c[1,j]:=DRV(x[i,0],x[i,j],j);
c[2,j]:=DRV(x[i,0]+.5*Delta!X,x[i,j]+.5*Delta!X*C[1,j],j);
c[3,j]:=DRV(x[i,0]+.5*Delta!X,x[i,j]+.2071068*Delta!X*C[1,j]
+.2928932*Delta!X*C[2,j] ,j);
c[4,j]:=DRV(x[i,0]+Delta!X,x[i,j]-.7071068*Delta!X*C[1,j]
+1.7071068*Delta!X*C[3,j] ,j);
c[i+1,j]:=x[i,j]
+Delta!X/6*
( C[1,j]
+.2928932*C[2,j]
+3.4142136*c[3,j]
+C[4,j]
);
END "Gill constants":
for j:=1 step 1 until N!Equations do
If(ABS( (C[3,j]-C[2,j])
/(C[2,j]-C[1,j])
)
>.01
)
Then Delta!X:=Delta!x/10;
END "Integration Stepping";
END "RUNGEK: Runge-Kutta-Gill fourth-order Integration Program";
Internal Real procedure realrd(string s);
begin "realrd"
! input real constants and assign to real variables;
real f; string a;
print(13&10&"Please input ",s,"> ");
a:=inchwl;f:=realscan(a,32);
return (f)
end "realrd";
Internal Procedure read!arai (reference real array a; integer n;string name);
begin"araird" string ary; integer i;
!similarly input arrays;
print (13&10&"Type in array,"&name&", elements sep by space,",
"<cr> when finished;",n,"elements "&13&10&" > ");
ary:=inchwl;
for i:=1 step 1 until n do a[i]:=realscan(ary,32);
end"araird";
Internal Procedure File!Read!Arai (REFERENCE REAL Array A;
INTEGER N,Channel,BreakTable);
begin"array-file" STRING Ary; INTEGER I;Ary:=NULL;
for i:=0 step 1 until n do
begin "ary-in"
if ary=null then ary:=input(Channel,Breaktable);
a[i]:=realscan(ary,32)
end "ary-in";
end"array-file";
Internal Simple Procedure Matrix!Print(REAL ARRAY A;
INTEGER Down,Accross;
STRING Name);
BEGIN "Print a two-dimensional array" Integer I,J;
Print (10&13&9&"Matrix:"&Name&13&10);
For I:=1 step 1 until Down do
BEGIN "Matrix looping"
for j:=1 step 1 until accross do
Print (a[i,j] ,32&null)
Print (13&10&Null);
END "Matrix looping";
END "Print a two-dimensional array";
Internal Real Procedure FIBONACCI(REAL Aa,Bb,Accuracy;
REAL PROCEDURE F);
begin"fibonacci"
! MAXIMISATION -give neg fct for min;
real new!fibonacci,fibonacci!old,temporary,tau,l,r,a,b;
fibonacci!old:=0;new!fibonacci:=1;l:=a:=aa;r:=b:=bb;
while abs((l-r)/r)>accuracy do
begin"fibonacci section"
new!fibonacci:=fibonacci!old+new!fibonacci
+O*(fibonacci!old:=new!fibonacci);
tau:=fibonacci!old/new!fibonacci;
l:=b-tau*abs(b-a);r:=a+tau*abs (b-a);
if f(l)>f(r) then b:=r else a:=l;
if a<aa or b>bb then done;
end"fibonacci section";
return((a+b)/2);
end"fibonacci";
Internal Procedure gauss!jordan (reference real array a,bb;integer n);
begin "gauss" real array c[1:n,1:n+1],b,d,e [1:n+1]
integer i,ipl,j,m,k,l,ko,km;real r,rl,s;label unitise,skip;
comment The Panagiotopoulos-Snyder matrix zonker;
for i:=1 step 1 until n do for j:=i step 1 until n do c[i,j]:=a[i,j];
for i:=1 step 1 until n do c[i,n+1]:=bb[i];
for i:=1 step 1 until n do
begin "gj1"
r:=c[i,i]; if (i geq n) then go unitise; m:=i;ip1:=i+1;
for j:=ip1 step 1 until n do
begin "gj2"
r1:=abs(c[m,i]); s:=abs(c[j,i]);if (r1<s) then m:=j;
end "gj2";
if(m=i) then go unitise; r:=c[m,i];
for j:=i step 1 until n+1 do
begin "gj3" b[j]:=c[i,j];c[i,j]:=c[m,j];c[m,j]:=b[j];end "gj3";
unitise: for k:=1 step 1 until n+1 do b[k]:=c[i,k]/r;
for k:=1 step 1 until n do
begin "gj4" s:=c[k,i]; if (s=0) then go skip;
for l:=1 step 1 until n+1 do d[l]:=-s*b[l];
for ko:=1 step 1 until n+1 do e[ko]:=c[k,ko]+d[ko];
for km:=1 step 1 until n+1 do c[k,km]:=e[km]; skip:
end "gj4";
for k:=1 step 1 until n+1 do c[i,k]:=b[k];
end "gj1";
for i:=1 step 1 until n do bb[i]:=c[i,n+1];
end "gauss";
Internal Real Procedure Simpson (real delta!x;real array y; integer n!rows);
begin "simpsn"
define oddevn(i)={i-2*(i%2)};
real s;integer i;
for i:=0 step 1 until n!rows do
if i=0 or i=n!rows then s:=s+y[i] else
if oddevn(i) then s:=s+4*y[i] else
s:=s+2*y[i];
s:=s*delta!x/3;return(s);
end"simpsn";
Internal Real Procedure Max!Array(Real Array Our!Vector;
Integer Dimension)
BEGIN "Maximum of Vector" Real Maximum; Maximum:=-INFINITY;
Dimension:=Dimension+1;
While (Dimension:=Dimension-1) do
if Our!Vector[Dimension]>Maximum then Maximum:=Our!Vector[Dimension];
Return (Maximum)
END "Maximum of Vector";
Internal Real Procedure Min!Array(Real Array Our!Vector;
Integer Dimension);
BEGIN "Minimum of Vector" Real Minimum: Minimum:=INFlNITY;
Dimension:=Dimension+1;
While (Dimension:=Dimension-1) do
if Our!Vector[Dimension]<Minimum then Minimum:=Our!Vector[Dimension];
Return (Minimum);
END "Minimum of Vector";
Internal Procedure Recirculation(Real Array X,R,T;
Reference Real Array Recirculation;
Integer Vector!Length);
begin "Cardiac Recirculation"
Integer I;
begin "loop"
Real Total!Mass,Mass!Integral,Max!Mass,Delta!T;
Integer I!T,Chan;
Real Array!Temp[0:Vector!Length];
delta!t:=t[2]-t[1];
total!mass:=(mass!integral:=simpson(delta!x,vector!length))
/delta!t/vector!length;
max!Mass:=Max!Array(x,vector!length);
for i:=0 step 1 until vector!length do
begin "recirc fro temp"
for i!t:=1 step 1 until vector!length-! do
if (i-i!t)geq 0 then
temp[i!t]:=x[i!t]*r[i-i!t]/max!mass;
recirculation[i]:=simpson(delta!t,temp,i);
end "recirc fro temp";
end"loop";
end"Cardiac Recirculation";
Internal Procedure Instaplot(REAL ARRAY X,Y,Window!Points,Text!XY!Location;
INTEGER Dimension,Picture!Number,Plotter!Type;
STRING PLX!File!Name,Text;
BOOLEAN Insta!Window,Want!Axes,Want!PLX!File,Want!Text,
Want!Ticks!Labeled,New!Picture);
BEGIN "InstaFlot"
integer i;
if not New!Picture
then dappend(Picture!Number)
else
begin "Open New Picture"
dini (Plotter!Type,0,0,0);dget;
If Insta!Window then
BEGIN "Instant Windowing"
Window!Points[1]:=Min!Array(x,dimension);
Window!Points[2]:=Max!Array(x,dimension);
Window!Points[3]:=Min!Array(y,dimension);
Window!Points[4]:=Max!Array(y,dimension);
END "Instant Windowing";
dwind(Window!Points[1],Window!Points[2],
Window!Points[3],Window!Points[4]);
dopen(Picture!Number);
if Want!Axes then
BEGIN "Plot Axes"
dvect(Window!Points[1],Window!Points[3],
Window!Points[2],Window!Points[3]);
dvect(Window!Points[1],Window!Points[3],
Window!Points[1],Window!Points[4]);
dvect(Window!Points[1],Window!Points[4],
Window!Points[2],Window!Points[4]);
dvect(Window!Points[2],Window!Points[3],
Window!Points[2],Window!Points[4]);
END "Plot Axes";
if Want!Ticks!Labeled then
BEGIN "Make Tick Marks"
Real X!Removed,Y!Removed;
X!Removed:=(Window!Points[2]
-Window!Points[1])
*(-.075);
Y!Removed:=(Window!Points[4]
-Window!Points [3])
*(-.075);
Setformat(3,2);
DMove(
Window!Points[1],
Window!Points[3]
-Y!Removed);
DText(CVG(
(Window!Points[1])));
DMove(
(Window!Points[2]
+Window!Points[1])/2,
Window!Points[3]
-Y!Removed);
DText(CVG((Window!Points[2]
+Window!Points[1])/2
));
DMove(
Window!Points[2],
Window!Points[3]
-Y!Removed);
DText(CVG((Window!Points[2])
));
DMove(
Window!Points[1]
-X!Removed,
Window!Points[3]);
DText(CVG((Window!Points[3])
));
DMove(
Window!Points[1]
-X!Removed,
(Window!Points[4]
+Window!Points[3])/2);
DText(CVG(
(Window!Points[4]
+Window!Points[3])/2
));
DMove(
Window!Points[1]
-X!Removed,
Window!Points[4]);
DText(CVG((Window!Points[4])
));
SetFormat(13,6);
END "Make Tick Marks";
end "Open New Picture";
for i:=1 step 1 until Dimension do ddraw (x[i] ,y[i]);
BEGIN "Close Plot"
INTEGER Terminal!Channel,File!Hexnumber;
REAL ARRAY Buffer!Pointer[1:128];
if Want!PLX!file then
BEGIN "PLX quiet term step one"
terminal!channel:=openfile("TTY:","RW");
Fiie!Hexnumber:=cvsix(PLX!File!Name)
!Terminal , Shut Up!!;
sfmod(Terminal!Channel,
'400000000000 LOR rfmod(Terminal!Channel));
END "PLX quiet term step one";
If Want!Text then
BEGIN "Make Text"
dmove(Text!XY!Location[1],Text!XY!Location[2]);
dtext(Text);
END "Make Text";
dget;dpost(Picture!Number); ddone;
If Want!PLX!File then dplot(Buffer!Pointer[1],File!Hexnumber);
dunpost(Picture!Number);
If Want!PLX!File then
BEGIN "PLX quiet term step two't
sfmod (Terminal!Channel,
'377777777777 LAND rfmod(Terminal!Channel))
cfile(Terminal!Channel);
END "PLX quiet term step two";
END "Close Plot";
END "Instaplot";
Internal SIMPLE PROCEDURE cc!trapper(procedure trap!proc);
begin "assembler routine"
start!code comment start assembler code;
movei 1,'400000; comment move octal immed to 1;
hrlzi 3,'400000; comment half word instuction;
epcap; comment enable capabilities;
end; comment end assembler code;
Psimap(34,trap!proc,0,1); comment chan 34 calls trap!proc;
Enable(34); comment enable channel 34;
Ati(34,3); comment use warning if ASCII #3;
end"assembler routine";
Internal simple procedure sink;
BEGIN "SINK:CUBI logout procedure"
start!code comment start assembler code;
dtach; comment detach current job;
movni 1,1; comment move neg immmed of 1 to 1
to specify current job;
lgout; comment now logout detached job
which spec'd by -1/curr;
end; comment end assembler code;
END "SINK:CUBI logout procedure";
Internal Procedure Structural Section;
BEGIN "Structural Section Properties, Interactive"
INTEGER No!Coordinates;
Real Array x,y,inertia[O:4]
REAL height,area,area!x,top!section!mod,bottom!section!mod;
no!coordinates:=realrd("Number of Coordinates");
height:=realrd("Section Height");
x[0]:=realrd("Starting X Coordinates");
y[0]:=realrd("Starting Y Coordinates");
While(no!coordinates:=no!coordinates-1) do
BEGIN "Centroid Averaging"
x[1]:=realrd("Next X");
y[1]:=realrd("Next Y");
x[2] =x[1]+x[0];
x[3]:=x[1]-x[0];
y[2]:=y[1]-y[0];
y[3]:=y[1]+y[0];
area:=area-y[2]*x[2]/2;
x[4]:=x[4]+y[2]/ 8*(x[2]^2+x[3]^2);
inertia[1]:=inertia[1]
+x[3]*y[3]/24*(y[2]^2+y[3]^2);
inertia[2]:=inertia[2]
-x[2]*y[2]/24*(x[2]^2+x[3]^2);
x[0]:=x[1];y[0]:=y[1];area!x:=-area;
END "Centroid Averaging";
x[4]:=x[4]/area!x;y[4]:=y[4]/area;
inertia[3]:=inertia[1]-area*y[4]^2;
inertia[4]:=inertia[2]-area*x[4]^2;
y[0]:=height-y[4];
top!section!mod:=inertia[3]/y[4];
bottom!section!mod:=inertia[3]/y[0];
PRINT(13&10&"STRUCTURAL SECTION PROPERTIES:
**ASSUMING HOMOGENEOUS DIMENSIONS** -
AREA:", area,
13&10&"X-Axis Moment of Inertia:" ,inertia[3],
13&10&"Centroid to Top:'t,y[0],
13&10&"Centroid to Bottom:",y[4],
13&10&"Bottom Section Modulus:",bottom!section!mod,
13&10&"Top Section Modulus:",top!Section!mod,12&0);
END "Structural Section Properties, Interactive";
Internal Procedure Monte!CarIo(REFERENCE REAL ARRAY X;
INTEGER Order, N!Step);
BEGIN "Monte Carlo Random Path"
INTEGER N,Power10;
x[1]:=1;
Power10:=10+Order;
For n:=1 step 1 until (n!step-1) do
x[n+1]:=(x[n]*4782969) mod power1O;
return;
END "Monte Carlo Random Path";
Internal Procedure Binary!Monte!Carlo(REFERENCE REAL ARRAY X;
INTEGER Order, N!Step,Large!Number)
BEGIN "Monte Carlo Random Path ,Binary"
INTEGER N,Power2;
x[1] =1;
Power2:=2+Order;
Large!Number:=Large!Number*8-3;
For n:=1 step 1 until (n!step-1) do
BEGIN "Binary Assembler Operations"
real xx,xxx;
xx:=x[n]
Start!Code comment begin assembler;
move 1,xx; comment put x[n] in 1;
imul 1,large!number; comment xx*larg..;
idiv 1,power2; comment div by power2;
movem 2,xxx; comment remainder in 2;
end; comment end assembler code;
x[n+1]:=xxx;
END "Binary Assembler Operations";
return;
END "Monte Carlo Random Path ,Binary";
Internal Procedure Gradient!Search(REAL PROCEDURE F;
REFERENCE REAL ARRAY X;
REAL Delta,Accuracy,Initial!Slope;
REAL ARRAY X!Upper,X!Lower;
INTEGER Vector!Length;
REFERENCE BOOLEAN Out!Of!Bounds);
BEGIN "Gradient Search acc. to Prof. Amir N. Nahavandi"
COMMENT : 0) Panagiotopoulos baffles them again!
1) F should return dependant variable for vector set of
all independant variables passed.
2) Delta passed is used as norm of Delta!X.
3) lf search goes out of bounds then nothing is returned and
Out!Of!Bounds is set true.
4) X is the vector of independant variables.
5) Vector Slope is the gradient for the current iteration.
It is initially set so that each element is equal to
a guessed slope, lnitial!Slope.
6) Algorithm due to class notes of ChE451OY81 as taught by
Professor Amir N. Nahavandi of CU ChE,AP/ENE&ME.
;
Real array X!Old,Slope,Delta!X[1:Vector!Length];
Boolean Converged,First!Iteration;
Out!Of!Bounds:=false;
Converged:=false;
First!Iteration:=true;
While not converged do
BEGIN "Search by max'd slope"
Real Slope!Norm,Delta!X!Norm,Avg!Delta!X;
Integer I,J;
If first!iteration
then
for i:=1 step 1 until Vector!Length do
Slope [i]:=Initial!Slope
else
for i:=1 step 1 until Vector!Length do
Slope[i]:=(f(x)-f(x!old))/Delta!X[i];
First!Iteration:=false;
Slope!Norm:=Delta!X!Norm:=Avg!Delta!X:=0;
Slope!Norm:=NORM(Slope,Vector!Length);
Delta!X[1] :=slope[1]/slope!norm*delta;
for i:=2 step 1 until vector!length do
delta!X[i]:=delta!x[1]*slope[i]/slope[1]
comment by condition of maximised slope, i.e.,
by setting derivatives of hyperplane
slopes to zero;
Delta!X!Norm:=NORM(Delta!X,Vector!Length);
Avg!Delta!X:=AVERAGE(Delta!X,Vector!Length);
If (avg!delta!x*vector!length/delta!x!norm<accuracy)
then converged:=true;
For i:=1 step 1 until vector!length do
if x[i]>x!upper[i] or x[i]<x!lower[i]then
out!of!bounds:=true;
for i:=1 step 1 until vector!length do
x[i]:=(x!old[i]:=x[i])+delta!x[i];
END "Search by max'd slope";
if not out!of!bounds then return;
END "Gradient Search acc. to Prof. Amir N. Nahavandi";
Internal REAL PROCEDURE GAMMA(REAL X);
BEGIN "gamma function"
return (
if x=1 then 1
else if x=.5 then 2*atan(1)^.5
else if x<0 then gamma(x+1)/x
else if x=0 then 1
else if (x mod 1)=0 then (x-1)*gamma(x-1)
else 1/(x+.57721*x^2-.65588*x^3-.042*x^4+.16654*x^5-.0422*x^6
-.00962*x^7+.0072189*x^8)
);
END "gamma function";
Internal Real Procedure Legendre!F(real x; integer n)
BEGIN" Legendre P polynomial trig expsn 3 terms acc. to Abramowitz&Stegun"
real p,theta;
theta:=arcos(x);
return (2^(2*n+1)*factorial(n)^2/factorial(2*n+1)
*(sin((n+1)*theta)+(n+1)/ (2*n+3)*sin((N+3)*theta)
+1.5*(n+1)*(n+1)/ (2*n+3)/ (2*n+5)*sin((n+5)*theta)
END "Legendre P polynomial trig expsn 3 terms acc. to Abramowitz&Stegun";
Internal REAL PROCEDURE Weighted!Gauss!Legendre!Integration(REAL ARRAY X!Vector;
INTEGER Vector!Length;
REAL PROCEDURE F;
REAL X!Initial,X!Final);
BEGIN "Gauss-Legendre Integration acc. to Abramowitz & Stegun"
REAL Sum;INTEGER I;REAL x;
for I:=1 step 1 until Vector!Length do
BEGIN "Gauss Legendre Weighting Factor Summing"
x:=(-x!Initial-X!Final +2*X!Vector[i])/(X!Final-X!Initial);
sum:=sum+f(x)*2/(1-x^2)/Legendre!p(X!Vector!Length)^2;
END "Gauss Legendre Weighting Factor Summing";
Return (sum/2*(X!Final-X!Initial));
END "Gauss-Legendre Integration acc. to Abramowitz & Stegun";
INTERNAL COMPLEX PROCEDURE C!Conjugate(COMPLEX z);
BEGIN "Complex Conjugate"
CMPLX:Im[Z]:=-CMPLX:Im[Z]
END "Complex Conjugate";
INTERNAL COMPLEX PROCEDURE C!Add(COMPLEX z1,z2);
BEGIN "Complex Addition"
COMPLEX Z; Z:=NEW!COMPLEX;
CMPLX:Re[Z]:=CMPLX:Re[Z1]+CMPLX:Re[Z2];
CMPLX:Im[Z]:=CMPLX:Im[Z1]+CMPLX:Im[Z2];
RETURN (Z);
END "Complex Addition";
INTERNAL COMPLEX PROCEDURE C!Subtract(COMPLEX z1,z2);
BEGIN "Complex Subtraction"
COMPLEX Z; Z:=NEW!COMPLEX;
CMPLX:Re[Z]:=CMPLX:Re[Z1]-CMPLX:Re[Z2];
CMPLX:Im[Z]:=CMPLX:Im[Z1]-CMPLX:Im[Z2];
RETURN(Z);
END "Complex Subtraction";
INTERNAL REAL PROCEDURE C!Abs(COMPLEX Z);
BEGIN "Absolute Value of Complex Number"
RETURN((CMPLX:Re[Z]^2 + CMPLX:Im[Z]^2)^.5);
END "Absolute Value of Complex Number";
INTERNAL COMPLEX PROCEDURE C!Multiply(COMPLEX Z1,Z2);
BEGIN "Complex Multiplication"
COMPLEX Z; Z:=NEW!COMPLEX;
CMPLX:Re[Z]:=CMPLX:Re[Z1]*CMPLX:Re[Z2]-CMPLX:Im[Z1]*CMPLX:Im[Z2];
CMPLX:Im[Z]:=CMPLX:Re[Z1]*CMPLX:Im[Z2]+CMPLX:Re[Z2]*CMPLX:Im[Z1];
RETURN(Z);
END "Complex Multiplication":
INTERNAL COMPLEX PROCEDURE C!Divide(COMPLEX Z1,Z2);
BEGIN "Complex Division"
REAL Denominator; COMPLEX Z; Z:=NEW!COMPLEX;
Denominator:=CMPLX:Re[Z2]^2-CMPLX:Im[Z2]^2
CMPLX:Re[Z]:=(CMPLX:Re[z1]*(CMPLX:Re[z2]+CMPLX:Im[z1]*CMPLX:Im[Z2]))
/Denominator;
CMPLX:Im[Z]:=(CMPLX:Im[z1]*(CMPLX:Re[z2]-CMPLX:Re[z1]*CMPLX:Im[Z2]))
/Denominator;
RETURN (z)
END "Complex Division";
INTERNAL COMPLEX PROCEDURE C!EXP(COMPLEX z):
BEGIN "Complex Eulerian Antilogarithm"
REAL R!Exponentiated,Theta; COMPLEX ZZ; ZZ:=NEW!COMPLEX;
R!Exponentiated:=EXP(C!Abs(Z))
Theta:=ATAN(-CMPLX:Im[Z]/CMPLX:Re[Z]);
CMPLX:Re[ZZ]:=COS(Theta)*R!Exponentiated;
CMPLX:Im[ZZ]:=SIN(Theta)*R!Exponentiated;
RETURN (ZZ);
END "Complex Eulerian Antilogarithm";
INTERNAL STRING PROCEDURE CvCmplx(COMPLEX z);
BEGIN "Complex to String Conversion for output"
RETURN (Cvs(CMPLX:Re[Z])
& (If CMPLX:Im[Z] < 0 then " - " else " + ") &
Cvs(CMPLX:Im[Z])&"i");
END "Complex to String Conversion for output";
END "SBRPC1";