\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} <>= )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} <>= )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} <>= --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. @ <<*>>= <> <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}