\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra d02Package.spad} \author{Brian Dupee} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{package ODEPACK AnnaOrdinaryDifferentialEquationPackage} <>= )abbrev package ODEPACK AnnaOrdinaryDifferentialEquationPackage ++ Author: Brian Dupee ++ Date Created: February 1995 ++ Date Last Updated: December 1997 ++ Basic Operations: solve, measure ++ Description: ++ \axiomType{AnnaOrdinaryDifferentialEquationPackage} is a \axiom{package} ++ of functions for the \axiom{category} \axiomType{OrdinaryDifferentialEquationsSolverCategory} ++ with \axiom{measure}, and \axiom{solve}. ++ EDF ==> Expression DoubleFloat LDF ==> List DoubleFloat MDF ==> Matrix DoubleFloat DF ==> DoubleFloat FI ==> Fraction Integer EFI ==> Expression Fraction Integer SOCDF ==> Segment OrderedCompletion DoubleFloat VEDF ==> Vector Expression DoubleFloat VEF ==> Vector Expression Float EF ==> Expression Float LF ==> List Float F ==> Float VDF ==> Vector DoubleFloat VMF ==> Vector MachineFloat MF ==> MachineFloat LS ==> List Symbol ST ==> String LST ==> List String INT ==> Integer RT ==> RoutinesTable ODEA ==> Record(xinit:DF,xend:DF,fn:VEDF,yinit:LDF,intvals:LDF,_ g:EDF,abserr:DF,relerr:DF) IFL ==> List(Record(ifail:Integer,instruction:String)) Entry ==> Record(chapter:String, type:String, domainName: String, defaultMin:F, measure:F, failList:IFL, explList:LST) Measure ==> Record(measure:F,name:String, explanations:List String) AnnaOrdinaryDifferentialEquationPackage(): with solve:(NumericalODEProblem) -> Result ++ solve(odeProblem) is a top level ANNA function to solve numerically a ++ system of ordinary differential equations i.e. equations for the ++ derivatives Y[1]'..Y[n]' defined in terms of X,Y[1]..Y[n], together ++ with starting values for X and Y[1]..Y[n] (called the initial ++ conditions), a final value of X, an accuracy requirement and any ++ intermediate points at which the result is required. ++ ++ It iterates over the \axiom{domains} of ++ \axiomType{OrdinaryDifferentialEquationsSolverCategory} ++ to get the name and other ++ relevant information of the the (domain of the) numerical ++ routine likely to be the most appropriate, ++ i.e. have the best \axiom{measure}. ++ ++ The method used to perform the numerical ++ process will be one of the routines contained in the NAG numerical ++ Library. The function predicts the likely most effective routine ++ by checking various attributes of the system of ODE's and calculating ++ a measure of compatibility of each routine to these attributes. ++ ++ It then calls the resulting `best' routine. solve:(NumericalODEProblem,RT) -> Result ++ solve(odeProblem,R) is a top level ANNA function to solve numerically a ++ system of ordinary differential equations i.e. equations for the ++ derivatives Y[1]'..Y[n]' defined in terms of X,Y[1]..Y[n], together ++ with starting values for X and Y[1]..Y[n] (called the initial ++ conditions), a final value of X, an accuracy requirement and any ++ intermediate points at which the result is required. ++ ++ It iterates over the \axiom{domains} of ++ \axiomType{OrdinaryDifferentialEquationsSolverCategory} contained in ++ the table of routines \axiom{R} to get the name and other ++ relevant information of the the (domain of the) numerical ++ routine likely to be the most appropriate, ++ i.e. have the best \axiom{measure}. ++ ++ The method used to perform the numerical ++ process will be one of the routines contained in the NAG numerical ++ Library. The function predicts the likely most effective routine ++ by checking various attributes of the system of ODE's and calculating ++ a measure of compatibility of each routine to these attributes. ++ ++ It then calls the resulting `best' routine. solve:(VEF,F,F,LF) -> Result ++ solve(f,xStart,xEnd,yInitial) is a top level ANNA function to solve numerically a ++ system of ordinary differential equations i.e. equations for the ++ derivatives Y[1]'..Y[n]' defined in terms of X,Y[1]..Y[n], together ++ with a starting value for X and Y[1]..Y[n] (called the initial ++ conditions) and a final value of X. A default value ++ is used for the accuracy requirement. ++ ++ It iterates over the \axiom{domains} of ++ \axiomType{OrdinaryDifferentialEquationsSolverCategory} contained in ++ the table of routines \axiom{R} to get the name and other ++ relevant information of the the (domain of the) numerical ++ routine likely to be the most appropriate, ++ i.e. have the best \axiom{measure}. ++ ++ The method used to perform the numerical ++ process will be one of the routines contained in the NAG numerical ++ Library. The function predicts the likely most effective routine ++ by checking various attributes of the system of ODE's and calculating ++ a measure of compatibility of each routine to these attributes. ++ ++ It then calls the resulting `best' routine. solve:(VEF,F,F,LF,F) -> Result ++ solve(f,xStart,xEnd,yInitial,tol) is a top level ANNA function to solve ++ numerically a system of ordinary differential equations, \axiom{f}, i.e. ++ equations for the derivatives Y[1]'..Y[n]' defined in terms ++ of X,Y[1]..Y[n] from \axiom{xStart} to \axiom{xEnd} with the initial ++ values for Y[1]..Y[n] (\axiom{yInitial}) to a tolerance \axiom{tol}. ++ ++ It iterates over the \axiom{domains} of ++ \axiomType{OrdinaryDifferentialEquationsSolverCategory} contained in ++ the table of routines \axiom{R} to get the name and other ++ relevant information of the the (domain of the) numerical ++ routine likely to be the most appropriate, ++ i.e. have the best \axiom{measure}. ++ ++ The method used to perform the numerical ++ process will be one of the routines contained in the NAG numerical ++ Library. The function predicts the likely most effective routine ++ by checking various attributes of the system of ODE's and calculating ++ a measure of compatibility of each routine to these attributes. ++ ++ It then calls the resulting `best' routine. solve:(VEF,F,F,LF,EF,F) -> Result ++ solve(f,xStart,xEnd,yInitial,G,tol) is a top level ANNA function to solve ++ numerically a system of ordinary differential equations, \axiom{f}, i.e. ++ equations for the derivatives Y[1]'..Y[n]' defined in terms ++ of X,Y[1]..Y[n] from \axiom{xStart} to \axiom{xEnd} with the initial ++ values for Y[1]..Y[n] (\axiom{yInitial}) to a tolerance \axiom{tol}. ++ The calculation will stop if the function G(X,Y[1],..,Y[n]) evaluates to zero before ++ X = xEnd. ++ ++ It iterates over the \axiom{domains} of ++ \axiomType{OrdinaryDifferentialEquationsSolverCategory} contained in ++ the table of routines \axiom{R} to get the name and other ++ relevant information of the the (domain of the) numerical ++ routine likely to be the most appropriate, ++ i.e. have the best \axiom{measure}. ++ ++ The method used to perform the numerical ++ process will be one of the routines contained in the NAG numerical ++ Library. The function predicts the likely most effective routine ++ by checking various attributes of the system of ODE's and calculating ++ a measure of compatibility of each routine to these attributes. ++ ++ It then calls the resulting `best' routine. solve:(VEF,F,F,LF,LF,F) -> Result ++ solve(f,xStart,xEnd,yInitial,intVals,tol) is a top level ANNA function to solve ++ numerically a system of ordinary differential equations, \axiom{f}, i.e. ++ equations for the derivatives Y[1]'..Y[n]' defined in terms ++ of X,Y[1]..Y[n] from \axiom{xStart} to \axiom{xEnd} with the initial ++ values for Y[1]..Y[n] (\axiom{yInitial}) to a tolerance \axiom{tol}. ++ The values of Y[1]..Y[n] will be output for the values of X in ++ \axiom{intVals}. ++ ++ It iterates over the \axiom{domains} of ++ \axiomType{OrdinaryDifferentialEquationsSolverCategory} contained in ++ the table of routines \axiom{R} to get the name and other ++ relevant information of the the (domain of the) numerical ++ routine likely to be the most appropriate, ++ i.e. have the best \axiom{measure}. ++ ++ The method used to perform the numerical ++ process will be one of the routines contained in the NAG numerical ++ Library. The function predicts the likely most effective routine ++ by checking various attributes of the system of ODE's and calculating ++ a measure of compatibility of each routine to these attributes. ++ ++ It then calls the resulting `best' routine. solve:(VEF,F,F,LF,EF,LF,F) -> Result ++ solve(f,xStart,xEnd,yInitial,G,intVals,tol) is a top level ANNA function to solve ++ numerically a system of ordinary differential equations, \axiom{f}, i.e. ++ equations for the derivatives Y[1]'..Y[n]' defined in terms ++ of X,Y[1]..Y[n] from \axiom{xStart} to \axiom{xEnd} with the initial ++ values for Y[1]..Y[n] (\axiom{yInitial}) to a tolerance \axiom{tol}. ++ The values of Y[1]..Y[n] will be output for the values of X in ++ \axiom{intVals}. The calculation will stop if the function ++ G(X,Y[1],..,Y[n]) evaluates to zero before X = xEnd. ++ ++ It iterates over the \axiom{domains} of ++ \axiomType{OrdinaryDifferentialEquationsSolverCategory} contained in ++ the table of routines \axiom{R} to get the name and other ++ relevant information of the the (domain of the) numerical ++ routine likely to be the most appropriate, ++ i.e. have the best \axiom{measure}. ++ ++ The method used to perform the numerical ++ process will be one of the routines contained in the NAG numerical ++ Library. The function predicts the likely most effective routine ++ by checking various attributes of the system of ODE's and calculating ++ a measure of compatibility of each routine to these attributes. ++ ++ It then calls the resulting `best' routine. solve:(VEF,F,F,LF,EF,LF,F,F) -> Result ++ solve(f,xStart,xEnd,yInitial,G,intVals,epsabs,epsrel) is a top level ANNA function to solve ++ numerically a system of ordinary differential equations, \axiom{f}, i.e. ++ equations for the derivatives Y[1]'..Y[n]' defined in terms ++ of X,Y[1]..Y[n] from \axiom{xStart} to \axiom{xEnd} with the initial ++ values for Y[1]..Y[n] (\axiom{yInitial}) to an absolute error ++ requirement \axiom{epsabs} and relative error \axiom{epsrel}. ++ The values of Y[1]..Y[n] will be output for the values of X in ++ \axiom{intVals}. The calculation will stop if the function ++ G(X,Y[1],..,Y[n]) evaluates to zero before X = xEnd. ++ ++ It iterates over the \axiom{domains} of ++ \axiomType{OrdinaryDifferentialEquationsSolverCategory} contained in ++ the table of routines \axiom{R} to get the name and other ++ relevant information of the the (domain of the) numerical ++ routine likely to be the most appropriate, ++ i.e. have the best \axiom{measure}. ++ ++ The method used to perform the numerical ++ process will be one of the routines contained in the NAG numerical ++ Library. The function predicts the likely most effective routine ++ by checking various attributes of the system of ODE's and calculating ++ a measure of compatibility of each routine to these attributes. ++ ++ It then calls the resulting `best' routine. measure:(NumericalODEProblem) -> Measure ++ measure(prob) is a top level ANNA function for identifying the most ++ appropriate numerical routine from those in the routines table ++ provided for solving the numerical ODE ++ problem defined by \axiom{prob}. ++ ++ It calls each \axiom{domain} of \axiom{category} ++ \axiomType{OrdinaryDifferentialEquationsSolverCategory} in turn to ++ calculate all measures and returns the best i.e. the name of ++ the most appropriate domain and any other relevant information. ++ It predicts the likely most effective NAG numerical ++ Library routine to solve the input set of ODEs ++ by checking various attributes of the system of ODEs and calculating ++ a measure of compatibility of each routine to these attributes. measure:(NumericalODEProblem,RT) -> Measure ++ measure(prob,R) is a top level ANNA function for identifying the most ++ appropriate numerical routine from those in the routines table ++ provided for solving the numerical ODE ++ problem defined by \axiom{prob}. ++ ++ It calls each \axiom{domain} listed in \axiom{R} of \axiom{category} ++ \axiomType{OrdinaryDifferentialEquationsSolverCategory} in turn to ++ calculate all measures and returns the best i.e. the name of ++ the most appropriate domain and any other relevant information. ++ It predicts the likely most effective NAG numerical ++ Library routine to solve the input set of ODEs ++ by checking various attributes of the system of ODEs and calculating ++ a measure of compatibility of each routine to these attributes. == add import ODEA,NumericalODEProblem f2df:F -> DF ef2edf:EF -> EDF preAnalysis:(ODEA,RT) -> RT zeroMeasure:Measure -> Result measureSpecific:(ST,RT,ODEA) -> Record(measure:F,explanations:ST) solveSpecific:(ODEA,ST) -> Result changeName:(Result,ST) -> Result recoverAfterFail:(ODEA,RT,Measure,Integer,Result) -> Record(a:Result,b:Measure) f2df(f:F):DF == (convert(f)@DF)$F ef2edf(f:EF):EDF == map(f2df,f)$ExpressionFunctions2(F,DF) preAnalysis(args:ODEA,t:RT):RT == rt := selectODEIVPRoutines(t)$RT if positive?(# variables(args.g)) then changeMeasure(rt,d02bbf@Symbol,getMeasure(rt,d02bbf@Symbol)*0.8) if positive?(# args.intvals) then changeMeasure(rt,d02bhf@Symbol,getMeasure(rt,d02bhf@Symbol)*0.8) rt zeroMeasure(m:Measure):Result == a := coerce(0$F)$AnyFunctions1(F) text := coerce("Zero Measure")$AnyFunctions1(ST) r := construct([[result@Symbol,a],[method@Symbol,text]])$Result concat(measure2Result m,r)$ExpertSystemToolsPackage measureSpecific(name:ST,R:RT,ode:ODEA):Record(measure:F,explanations:ST) == name = "d02bbfAnnaType" => measure(R,ode)$d02bbfAnnaType name = "d02bhfAnnaType" => measure(R,ode)$d02bhfAnnaType name = "d02cjfAnnaType" => measure(R,ode)$d02cjfAnnaType name = "d02ejfAnnaType" => measure(R,ode)$d02ejfAnnaType error("measureSpecific","invalid type name: " name)$ErrorFunctions measure(Ode:NumericalODEProblem,R:RT):Measure == ode:ODEA := retract(Ode)$NumericalODEProblem sofar := 0$F best := "none" :: ST routs := copy R routs := preAnalysis(ode,routs) empty?(routs)$RT => error("measure", "no routines found")$ErrorFunctions rout := inspect(routs)$RT e := retract(rout.entry)$AnyFunctions1(Entry) meth := empty()$LST for i in 1..# routs repeat rout := extract!(routs)$RT e := retract(rout.entry)$AnyFunctions1(Entry) n := e.domainName if e.defaultMin > sofar then m := measureSpecific(n,R,ode) if m.measure > sofar then sofar := m.measure best := n str:LST := [string(rout.key)$Symbol "measure: " outputMeasure(m.measure)$ExpertSystemToolsPackage " - " m.explanations] else str := [string(rout.key)$Symbol " is no better than other routines"] meth := append(meth,str)$LST [sofar,best,meth] measure(ode:NumericalODEProblem):Measure == measure(ode,routines()$RT) solveSpecific(ode:ODEA,n:ST):Result == n = "d02bbfAnnaType" => ODESolve(ode)$d02bbfAnnaType n = "d02bhfAnnaType" => ODESolve(ode)$d02bhfAnnaType n = "d02cjfAnnaType" => ODESolve(ode)$d02cjfAnnaType n = "d02ejfAnnaType" => ODESolve(ode)$d02ejfAnnaType error("solveSpecific","invalid type name: " n)$ErrorFunctions changeName(ans:Result,name:ST):Result == sy:Symbol := coerce(name "Answer")$Symbol anyAns:Any := coerce(ans)$AnyFunctions1(Result) construct([[sy,anyAns]])$Result recoverAfterFail(ode:ODEA,routs:RT,m:Measure,iint:Integer,r:Result): Record(a:Result,b:Measure) == while positive?(iint) repeat routineName := m.name s := recoverAfterFail(routs,routineName(1..6),iint)$RT s case "failed" => iint := 0 if s = "increase tolerance" then ode.relerr := ode.relerr*(10.0::DF) ode.abserr := ode.abserr*(10.0::DF) if s = "decrease tolerance" then ode.relerr := ode.relerr/(10.0::DF) ode.abserr := ode.abserr/(10.0::DF) (s = "no action")@Boolean => iint := 0 fl := coerce(s)$AnyFunctions1(ST) flrec:Record(key:Symbol,entry:Any):=[failure@Symbol,fl] m2 := measure(ode::NumericalODEProblem,routs) zero?(m2.measure) => iint := 0 r2:Result := solveSpecific(ode,m2.name) m := m2 insert!(flrec,r2)$Result r := concat(r2,changeName(r,routineName))$ExpertSystemToolsPackage iany := search(ifail@Symbol,r2)$Result iany case "failed" => iint := 0 iint := retract(iany)$AnyFunctions1(Integer) [r,m] solve(Ode:NumericalODEProblem,t:RT):Result == ode:ODEA := retract(Ode)$NumericalODEProblem routs := copy(t)$RT m := measure(Ode,routs) zero?(m.measure) => zeroMeasure m r := solveSpecific(ode,n := m.name) iany := search(ifail@Symbol,r)$Result iint := 0$Integer if (iany case Any) then iint := retract(iany)$AnyFunctions1(Integer) if positive?(iint) then tu:Record(a:Result,b:Measure) := recoverAfterFail(ode,routs,m,iint,r) r := tu.a m := tu.b r := concat(measure2Result m,r)$ExpertSystemToolsPackage expl := getExplanations(routs,n(1..6))$RoutinesTable expla := coerce(expl)$AnyFunctions1(LST) explaa:Record(key:Symbol,entry:Any) := ["explanations"::Symbol,expla] r := concat(construct([explaa]),r) iflist := showIntensityFunctions(ode)$ODEIntensityFunctionsTable iflist case "failed" => r concat(iflist2Result iflist, r)$ExpertSystemToolsPackage solve(ode:NumericalODEProblem):Result == solve(ode,routines()$RT) solve(f:VEF,xStart:F,xEnd:F,yInitial:LF,G:EF,intVals:LF,epsabs:F,epsrel:F):Result == d:ODEA := [f2df xStart,f2df xEnd,vector([ef2edf e for e in members f])$VEDF, [f2df i for i in yInitial], [f2df j for j in intVals], ef2edf G,f2df epsabs,f2df epsrel] solve(d::NumericalODEProblem,routines()$RT) solve(f:VEF,xStart:F,xEnd:F,yInitial:LF,G:EF,intVals:LF,tol:F):Result == solve(f,xStart,xEnd,yInitial,G,intVals,tol,tol) solve(f:VEF,xStart:F,xEnd:F,yInitial:LF,intVals:LF,tol:F):Result == solve(f,xStart,xEnd,yInitial,1$EF,intVals,tol) solve(f:VEF,xStart:F,xEnd:F,y:LF,G:EF,tol:F):Result == solve(f,xStart,xEnd,y,G,empty()$LF,tol) solve(f:VEF,xStart:F,xEnd:F,yInitial:LF,tol:F):Result == solve(f,xStart,xEnd,yInitial,1$EF,empty()$LF,tol) solve(f:VEF,xStart:F,xEnd:F,yInitial:LF):Result == solve(f,xStart,xEnd,yInitial,1.0e-4) @ \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}