diff options
Diffstat (limited to 'src/algebra/d02agents.spad.pamphlet')
| -rw-r--r-- | src/algebra/d02agents.spad.pamphlet | 424 |
1 files changed, 424 insertions, 0 deletions
diff --git a/src/algebra/d02agents.spad.pamphlet b/src/algebra/d02agents.spad.pamphlet new file mode 100644 index 00000000..d3e18ec4 --- /dev/null +++ b/src/algebra/d02agents.spad.pamphlet @@ -0,0 +1,424 @@ +\documentclass{article} +\usepackage{axiom} +\begin{document} +\title{\$SPAD/src/algebra d02agents.spad} +\author{Brian Dupee} +\maketitle +\begin{abstract} +\end{abstract} +\eject +\tableofcontents +\eject +\section{domain ODEIFTBL ODEIntensityFunctionsTable} +<<domain ODEIFTBL ODEIntensityFunctionsTable>>= +)abbrev domain ODEIFTBL ODEIntensityFunctionsTable +++ Author: Brian Dupee +++ Date Created: May 1994 +++ Date Last Updated: January 1996 +++ Basic Operations: showTheIFTable, insert! +++ Description: +++ \axiom{ODEIntensityFunctionsTable()} provides a dynamic table and a set of +++ functions to store details found out about sets of ODE's. + +ODEIntensityFunctionsTable(): E == I where + LEDF ==> List Expression DoubleFloat + LEEDF ==> List Equation Expression DoubleFloat + EEDF ==> Equation Expression DoubleFloat + VEDF ==> Vector Expression DoubleFloat + MEDF ==> Matrix Expression DoubleFloat + MDF ==> Matrix DoubleFloat + EDF ==> Expression DoubleFloat + DF ==> DoubleFloat + F ==> Float + INT ==> Integer + CDF ==> Complex DoubleFloat + LDF ==> List DoubleFloat + LF ==> List Float + S ==> Symbol + LS ==> List Symbol + MFI ==> Matrix Fraction Integer + LFI ==> List Fraction Integer + FI ==> Fraction Integer + ODEA ==> Record(xinit:DF,xend:DF,fn:VEDF,yinit:LDF,intvals:LDF,g:EDF,abserr:DF,relerr:DF) + ON ==> Record(additions:INT,multiplications:INT,exponentiations:INT,functionCalls:INT) + RVE ==> Record(val:EDF,exponent:INT) + RSS ==> Record(stiffnessFactor:F,stabilityFactor:F) + ATT ==> Record(stiffness:F,stability:F,expense:F,accuracy:F,intermediateResults:F) + ROA ==> Record(key:ODEA,entry:ATT) + + E ==> with + showTheIFTable:() -> $ + ++ showTheIFTable() returns the current table of intensity functions. + clearTheIFTable : () -> Void + ++ clearTheIFTable() clears the current table of intensity functions. + keys : $ -> List(ODEA) + ++ keys(tab) returns the list of keys of f + iFTable: List Record(key:ODEA,entry:ATT) -> $ + ++ iFTable(l) creates an intensity-functions table from the elements + ++ of l. + insert!:Record(key:ODEA,entry:ATT) -> $ + ++ insert!(r) inserts an entry r into theIFTable + showIntensityFunctions:ODEA -> Union(ATT,"failed") + ++ showIntensityFunctions(k) returns the entries in the + ++ table of intensity functions k. + + I ==> add + Rep := Table(ODEA,ATT) + import Rep + + theIFTable:$ := empty()$Rep + + showTheIFTable():$ == + theIFTable + + clearTheIFTable():Void == + theIFTable := empty()$Rep + void()$Void + + iFTable(l:List Record(key:ODEA,entry:ATT)):$ == + theIFTable := table(l)$Rep + + insert!(r:Record(key:ODEA,entry:ATT)):$ == + insert!(r,theIFTable)$Rep + + keys(t:$):List ODEA == + keys(t)$Rep + + showIntensityFunctions(k:ODEA):Union(ATT,"failed") == + search(k,theIFTable)$Rep + +@ +\section{package D02AGNT d02AgentsPackage} +<<package D02AGNT d02AgentsPackage>>= +)abbrev package D02AGNT d02AgentsPackage +++ Author: Brian Dupee +++ Date Created: May 1994 +++ Date Last Updated: January 1997 +++ Basic Operations: stiffnessFactor, jacobian +++ Description: +++ \axiom{d02AgentsPackage} contains a set of computational agents +++ for use with Ordinary Differential Equation solvers. +d02AgentsPackage(): E == I where + LEDF ==> List Expression DoubleFloat + LEEDF ==> List Equation Expression DoubleFloat + EEDF ==> Equation Expression DoubleFloat + VEDF ==> Vector Expression DoubleFloat + MEDF ==> Matrix Expression DoubleFloat + MDF ==> Matrix DoubleFloat + EDF ==> Expression DoubleFloat + DF ==> DoubleFloat + F ==> Float + INT ==> Integer + CDF ==> Complex DoubleFloat + LDF ==> List DoubleFloat + LF ==> List Float + S ==> Symbol + LS ==> List Symbol + MFI ==> Matrix Fraction Integer + LFI ==> List Fraction Integer + FI ==> Fraction Integer + ODEA ==> Record(xinit:DF,xend:DF,fn:VEDF,yinit:LDF,intvals:LDF,g:EDF,abserr:DF,relerr:DF) + ON ==> Record(additions:INT,multiplications:INT,exponentiations:INT,functionCalls:INT) + RVE ==> Record(val:EDF,exponent:INT) + RSS ==> Record(stiffnessFactor:F,stabilityFactor:F) + ATT ==> Record(stiffness:F,stability:F,expense:F,accuracy:F,intermediateResults:F) + ROA ==> Record(key:ODEA,entry:ATT) + + E ==> with + combineFeatureCompatibility: (F,F) -> F + ++ combineFeatureCompatibility(C1,C2) is for interacting attributes + combineFeatureCompatibility: (F,LF) -> F + ++ combineFeatureCompatibility(C1,L) is for interacting attributes + sparsityIF: MEDF -> F + ++ sparsityIF(m) calculates the sparsity of a jacobian matrix + jacobian: (VEDF,LS) -> MEDF + ++ jacobian(v,w) is a local function to make a jacobian matrix + eval: (MEDF,LS,VEDF) -> MEDF + ++ eval(mat,symbols,values) evaluates a multivariable matrix at given values + ++ for each of a list of variables + stiffnessAndStabilityFactor: MEDF -> RSS + ++ stiffnessAndStabilityFactor(me) calculates the stability and + ++ stiffness factor of a system of first-order differential equations + ++ (by evaluating the maximum difference in the real parts of the + ++ negative eigenvalues of the jacobian of the system for which O(10) + ++ equates to mildly stiff wheras stiffness ratios of O(10^6) are not + ++ uncommon) and whether the system is likely to show any oscillations + ++ (identified by the closeness to the imaginary axis of the complex + ++ eigenvalues of the jacobian). + stiffnessAndStabilityOfODEIF:ODEA -> RSS + ++ stiffnessAndStabilityOfODEIF(ode) calculates the intensity values + ++ of stiffness of a system of first-order differential equations + ++ (by evaluating the maximum difference in the real parts of the + ++ negative eigenvalues of the jacobian of the system for which O(10) + ++ equates to mildly stiff wheras stiffness ratios of O(10^6) are not + ++ uncommon) and whether the system is likely to show any oscillations + ++ (identified by the closeness to the imaginary axis of the complex + ++ eigenvalues of the jacobian). + ++ + ++ It returns two values in the range [0,1]. + systemSizeIF:ODEA -> F + ++ systemSizeIF(ode) returns the intensity value of the size of + ++ the system of ODEs. 20 equations corresponds to the neutral + ++ value. It returns a value in the range [0,1]. + expenseOfEvaluationIF:ODEA -> F + ++ expenseOfEvaluationIF(o) returns the intensity value of the + ++ cost of evaluating the input ODE. This is in terms of the number + ++ of ``operational units''. It returns a value in the range + ++ [0,1].\newline\indent{20} + ++ 400 ``operation units'' -> 0.75 \newline + ++ 200 ``operation units'' -> 0.5 \newline + ++ 83 ``operation units'' -> 0.25 \newline\indent{15} + ++ exponentiation = 4 units , function calls = 10 units. + accuracyIF:ODEA -> F + ++ accuracyIF(o) returns the intensity value of the accuracy + ++ requirements of the input ODE. A request of accuracy of 10^-6 + ++ corresponds to the neutral intensity. It returns a value + ++ in the range [0,1]. + intermediateResultsIF:ODEA -> F + ++ intermediateResultsIF(o) returns a value corresponding to the + ++ required number of intermediate results required and, therefore, + ++ an indication of how much this would affect the step-length of the + ++ calculation. It returns a value in the range [0,1]. + + I ==> add + + import ExpertSystemToolsPackage + + accuracyFactor:ODEA -> F + expenseOfEvaluation:ODEA -> F + eval1:(LEDF,LEEDF) -> LEDF + stiffnessAndStabilityOfODE:ODEA -> RSS + intermediateResultsFactor:ODEA -> F + leastStabilityAngle:(LDF,LDF) -> F + + intermediateResultsFactor(ode:ODEA):F == + resultsRequirement := #(ode.intvals) + (1.0-exp(-(resultsRequirement::F)/50.0)$F) + + intermediateResultsIF(o:ODEA):F == + ode := copy o + (t := showIntensityFunctions(ode)$ODEIntensityFunctionsTable) case ATT => + s := coerce(t)@ATT + negative?(s.intermediateResults)$F => + s.intermediateResults := intermediateResultsFactor(ode) + r:ROA := [ode,s] + insert!(r)$ODEIntensityFunctionsTable + s.intermediateResults + s.intermediateResults + a:ATT := [-1.0,-1.0,-1.0,-1.0,e:=intermediateResultsFactor(ode)] + r:ROA := [ode,a] + insert!(r)$ODEIntensityFunctionsTable + e + + accuracyFactor(ode:ODEA):F == + accuracyRequirements := convert(ode.abserr)@F + if zero?(accuracyRequirements) then + accuracyRequirements := convert(ode.relerr)@F + val := inv(accuracyRequirements)$F + n := log10(val)$F + (1.0-exp(-(n/(2.0))**2/(15.0))$F) + + accuracyIF(o:ODEA):F == + ode := copy o + (t := showIntensityFunctions(ode)$ODEIntensityFunctionsTable) case ATT => + s := coerce(t)@ATT + negative?(s.accuracy)$F => + s.accuracy := accuracyFactor(ode) + r:ROA := [ode,s] + insert!(r)$ODEIntensityFunctionsTable + s.accuracy + s.accuracy + a:ATT := [-1.0,-1.0,-1.0,e:=accuracyFactor(ode),-1.0] + r:ROA := [ode,a] + insert!(r)$ODEIntensityFunctionsTable + e + + systemSizeIF(ode:ODEA):F == + n := #(ode.fn) + (1.0-exp((-n::F/75.0))$F) + + expenseOfEvaluation(o:ODEA):F == + -- expense of evaluation of an ODE -- <0.3 inexpensive - 0.5 neutral - >0.7 very expensive + -- 400 `operation units' -> 0.75 + -- 200 `operation units' -> 0.5 + -- 83 `operation units' -> 0.25 + -- ** = 4 units , function calls = 10 units. + ode := copy o.fn + expenseOfEvaluation(ode) + + expenseOfEvaluationIF(o:ODEA):F == + ode := copy o + (t := showIntensityFunctions(ode)$ODEIntensityFunctionsTable) case ATT => + s := coerce(t)@ATT + negative?(s.expense)$F => + s.expense := expenseOfEvaluation(ode) + r:ROA := [ode,s] + insert!(r)$ODEIntensityFunctionsTable + s.expense + s.expense + a:ATT := [-1.0,-1.0,e:=expenseOfEvaluation(ode),-1.0,-1.0] + r:ROA := [ode,a] + insert!(r)$ODEIntensityFunctionsTable + e + + leastStabilityAngle(realPartsList:LDF,imagPartsList:LDF):F == + complexList := [complex(u,v)$CDF for u in realPartsList for v in imagPartsList] + argumentList := [abs((abs(argument(u)$CDF)$DF)-(pi()$DF)/2)$DF for u in complexList] + sortedArgumentList := sort(argumentList)$LDF + list := [u for u in sortedArgumentList | not zero?(u) ] + empty?(list)$LDF => 0$F + convert(first(list)$LDF)@F + + stiffnessAndStabilityFactor(me:MEDF):RSS == + + -- search first for real eigenvalues of the jacobian (symbolically) + -- if the system isn't too big + r:INT := ncols(me)$MEDF + b:Boolean := ((# me) < 150) + if b then + mc:MFI := map(edf2fi,me)$ExpertSystemToolsPackage2(EDF,FI) + e:LFI := realEigenvalues(mc,1/100)$NumericRealEigenPackage(FI) + b := ((# e) >= r-1)@Boolean + b => + -- if all the eigenvalues are real, find negative ones + e := sort(neglist(e)$ExpertSystemToolsPackage1(FI)) + -- if there are two or more, calculate stiffness ratio + ((n:=#e)>1)@Boolean => [coerce(e.1/e.n)@F,0$F] + -- otherwise stiffness not present + [0$F,0$F] + + md:MDF := map(edf2df,me)$ExpertSystemToolsPackage2(EDF,DF) + + -- otherwise calculate numerically the complex eigenvalues + -- using NAG routine f02aff. + + res:Result := f02aff(r,r,md,-1)$NagEigenPackage + realParts:Union(Any,"failed") := search(rr::Symbol,res)$Result + realParts case "failed" => [0$F,0$F] + realPartsMatrix:MDF := retract(realParts)$AnyFunctions1(MDF) -- array === matrix + imagParts:Union(Any,"failed") := search(ri::Symbol,res)$Result + imagParts case "failed" => [0$F,0$F] + imagPartsMatrix:MDF := retract(imagParts)$AnyFunctions1(MDF) -- array === matrix + imagPartsList:LDF := members(imagPartsMatrix)$MDF + realPartsList:LDF := members(realPartsMatrix)$MDF + stabilityAngle := leastStabilityAngle(realPartsList,imagPartsList) + negRealPartsList := sort(neglist(realPartsList)$ExpertSystemToolsPackage1(DF)) + empty?(negRealPartsList)$LDF => [0$F,stabilityAngle] + ((n:=#negRealPartsList)>1)@Boolean => + out := convert(negRealPartsList.1/negRealPartsList.n)@F + [out,stabilityAngle] -- calculate stiffness ratio + [-convert(negRealPartsList.1)@F,stabilityAngle] + + eval1(l:LEDF,e:LEEDF):LEDF == + [eval(u,e)$EDF for u in l] + + eval(mat:MEDF,symbols:LS,values:VEDF):MEDF == + l := listOfLists(mat) + ledf := entries(values)$VEDF + e := [equation(u::EDF,v)$EEDF for u in symbols for v in ledf] + l := [eval1(w,e) for w in l] + matrix l + + combineFeatureCompatibility(C1:F,C2:F):F == + + -- C1 C2 + -- s(C1,C2) = ----------------------- + -- C1 C2 + (1 - C1)(1 - C2) + + C1*C2/((C1*C2)+(1$F-C1)*(1$F-C2)) + + combineFeatureCompatibility(C1:F,L:LF):F == + + empty?(L)$LF => C1 + C2 := combineFeatureCompatibility(C1,first(L)$LF) + combineFeatureCompatibility(C2,rest(L)$LF) + + jacobian(v:VEDF,w:LS):Matrix EDF == + jacobian(v,w)$MultiVariableCalculusFunctions(S,EDF,VEDF,LS) + + sparsityIF(m:Matrix EDF):F == + l:LEDF :=parts m + z:LEDF := [u for u in l | zero?(u)$EDF] + ((#z)::F/(#l)::F) + + sum(a:EDF,b:EDF):EDF == a+b + + stiffnessAndStabilityOfODE(ode:ODEA):RSS == + odefns := copy ode.fn + ls:LS := [subscript(Y,[coerce(n)])$Symbol for n in 1..# odefns] + yvals := copy ode.yinit + for i in 1..#yvals repeat + zero?(yvals.i) => yvals.i := 0.1::DF + yexpr := [coerce(v)@EDF for v in yvals] + yv:VEDF := vector(yexpr) + j1:MEDF := jacobian(odefns,ls) + ej1:MEDF := eval(j1,ls,yv) + ej1 := eval(ej1,variables(reduce(sum,members(ej1)$MEDF)),vector([(ode.xinit)::EDF])) + ssf := stiffnessAndStabilityFactor(ej1) + stability := 1.0-sqrt((ssf.stabilityFactor)*(2.0)/(pi()$F)) + stiffness := (1.0)-exp(-(ssf.stiffnessFactor)/(500.0)) + [stiffness,stability] + + stiffnessAndStabilityOfODEIF(ode:ODEA):RSS == + odefn := copy ode + (t := showIntensityFunctions(odefn)$ODEIntensityFunctionsTable) case ATT => + s:ATT := coerce(t)@ATT + negative?(s.stiffness)$F => + ssf:RSS := stiffnessAndStabilityOfODE(odefn) + s := [ssf.stiffnessFactor,ssf.stabilityFactor,s.expense, + s.accuracy,s.intermediateResults] + r:ROA := [odefn,s] + insert!(r)$ODEIntensityFunctionsTable + ssf + [s.stiffness,s.stability] + ssf:RSS := stiffnessAndStabilityOfODE(odefn) + s:ATT := [ssf.stiffnessFactor,ssf.stabilityFactor,-1.0,-1.0,-1.0] + r:ROA := [odefn,s] + insert!(r)$ODEIntensityFunctionsTable + ssf + +@ +\section{License} +<<license>>= +--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. +--All rights reserved. +-- +--Redistribution and use in source and binary forms, with or without +--modification, are permitted provided that the following conditions are +--met: +-- +-- - Redistributions of source code must retain the above copyright +-- notice, this list of conditions and the following disclaimer. +-- +-- - Redistributions in binary form must reproduce the above copyright +-- notice, this list of conditions and the following disclaimer in +-- the documentation and/or other materials provided with the +-- distribution. +-- +-- - Neither the name of The Numerical ALgorithms Group Ltd. nor the +-- names of its contributors may be used to endorse or promote products +-- derived from this software without specific prior written permission. +-- +--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER +--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +@ +<<*>>= +<<license>> + +<<domain ODEIFTBL ODEIntensityFunctionsTable>> +<<package D02AGNT d02AgentsPackage>> +@ +\eject +\begin{thebibliography}{99} +\bibitem{1} nothing +\end{thebibliography} +\end{document} |