diff options
Diffstat (limited to 'src/algebra/d02agents.spad.pamphlet')
| -rw-r--r-- | src/algebra/d02agents.spad.pamphlet | 423 |
1 files changed, 0 insertions, 423 deletions
diff --git a/src/algebra/d02agents.spad.pamphlet b/src/algebra/d02agents.spad.pamphlet deleted file mode 100644 index e315970a..00000000 --- a/src/algebra/d02agents.spad.pamphlet +++ /dev/null @@ -1,423 +0,0 @@ -\documentclass{article} -\usepackage{open-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 - - 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} |