path: root/src/algebra/d02routine.spad.pamphlet

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra d02routine.spad}
\author{Brian Dupee}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{domain D02BBFA d02bbfAnnaType}
<<domain D02BBFA d02bbfAnnaType>>=
)abbrev domain D02BBFA d02bbfAnnaType
++ Author: Brian Dupee
++ Date Created: February 1995
++ Date Last Updated: January 1996
++ Basic Operations: 
++ Description:
++ \axiomType{d02bbfAnnaType} is a domain of 
++ \axiomType{OrdinaryDifferentialEquationsInitialValueProblemSolverCategory}
++ for the NAG routine D02BBF, a ODE routine which uses an  
++ Runge-Kutta method to solve a system of differential 
++ equations.  The function \axiomFun{measure} measures the
++ usefulness of the routine D02BBF for the given problem.  The 
++ function \axiomFun{ODESolve} performs the integration by using 
++ \axiomType{NagOrdinaryDifferentialEquationsPackage}.


d02bbfAnnaType():OrdinaryDifferentialEquationsSolverCategory == Result add  -- Runge Kutta

  EDF	==> Expression DoubleFloat
  LDF	==> List DoubleFloat
  MDF	==> Matrix DoubleFloat
  DF	==> DoubleFloat
  F	==> Float
  FI	==> Fraction Integer
  EFI	==> Expression Fraction Integer
  SOCDF	==> Segment OrderedCompletion DoubleFloat
  VEDF	==> Vector Expression DoubleFloat
  VEF	==> Vector Expression Float
  EF	==> Expression Float
  VDF	==> Vector DoubleFloat
  VMF	==> Vector MachineFloat
  MF	==> MachineFloat
  ODEA	==> Record(xinit:DF,xend:DF,fn:VEDF,yinit:LDF,intvals:LDF,_
                      g:EDF,abserr:DF,relerr:DF)
  RSS	==> Record(stiffnessFactor:F,stabilityFactor:F)
  INT	==> Integer
  EF2	==> ExpressionFunctions2

  import d02AgentsPackage, NagOrdinaryDifferentialEquationsPackage
  import AttributeButtons

  accuracyCF(ode:ODEA):F ==
    b := getButtonValue("d02bbf","accuracy")$AttributeButtons
    accuracyIntensityValue := combineFeatureCompatibility(b,accuracyIF(ode))
    accuracyIntensityValue > 0.999 => 0$F
    0.8*exp(-((10*accuracyIntensityValue)**3)$F/266)$F

  stiffnessCF(stiffnessIntensityValue:F):F ==
    b := getButtonValue("d02bbf","stiffness")$AttributeButtons
    0.5*exp(-(2*combineFeatureCompatibility(b,stiffnessIntensityValue))**2)$F

  stabilityCF(stabilityIntensityValue:F):F ==
    b := getButtonValue("d02bbf","stability")$AttributeButtons
    0.5 * cos(combineFeatureCompatibility(b,stabilityIntensityValue))$F

  expenseOfEvaluationCF(ode:ODEA):F ==
    b := getButtonValue("d02bbf","expense")$AttributeButtons
    expenseOfEvaluationIntensityValue := 
      combineFeatureCompatibility(b,expenseOfEvaluationIF(ode))
    0.35+0.2*exp(-(2.0*expenseOfEvaluationIntensityValue)**3)$F
    
  measure(R:RoutinesTable,args:ODEA) ==
    m := getMeasure(R,d02bbf :: Symbol)$RoutinesTable
    ssf := stiffnessAndStabilityOfODEIF args
    m := combineFeatureCompatibility(m,[accuracyCF(args),
            stiffnessCF(ssf.stiffnessFactor),
              expenseOfEvaluationCF(args),
                stabilityCF(ssf.stabilityFactor)])
    [m,"Runge-Kutta Merson method"]

  ODESolve(ode:ODEA) ==
    i:LDF := ode.intvals
    M := inc(# i)$INT
    irelab := 0$INT
    if positive?(a := ode.abserr) then 
      inc(irelab)$INT
    if positive?(r := ode.relerr) then
      inc(irelab)$INT
    if positive?(a+r) then
      tol:DF := a + r
    else
      tol := float(1,-4,10)$DF
    asp7:Union(fn:FileName,fp:Asp7(FCN)) :=
      [retract(vedf2vef(ode.fn)$ExpertSystemToolsPackage)$Asp7(FCN)]
    asp8:Union(fn:FileName,fp:Asp8(OUTPUT)) := 
      [coerce(ldf2vmf(i)$ExpertSystemToolsPackage)$Asp8(OUTPUT)]
    d02bbf(ode.xend,M,# ode.fn,irelab,ode.xinit,matrix([ode.yinit])$MDF,
             tol,-1,asp7,asp8)

@
\section{domain D02BHFA d02bhfAnnaType}
<<domain D02BHFA d02bhfAnnaType>>=
)abbrev domain D02BHFA d02bhfAnnaType
++ Author: Brian Dupee
++ Date Created: February 1995
++ Date Last Updated: January 1996
++ Basic Operations: 
++ Description:
++ \axiomType{d02bhfAnnaType} is a domain of 
++ \axiomType{OrdinaryDifferentialEquationsInitialValueProblemSolverCategory}
++ for the NAG routine D02BHF, a ODE routine which uses an  
++ Runge-Kutta method to solve a system of differential 
++ equations.  The function \axiomFun{measure} measures the
++ usefulness of the routine D02BHF for the given problem.  The 
++ function \axiomFun{ODESolve} performs the integration by using 
++ \axiomType{NagOrdinaryDifferentialEquationsPackage}.

d02bhfAnnaType():OrdinaryDifferentialEquationsSolverCategory == Result add  -- Runge Kutta
  EDF	==> Expression DoubleFloat
  LDF	==> List DoubleFloat
  MDF	==> Matrix DoubleFloat
  DF	==> DoubleFloat
  F	==> Float
  FI	==> Fraction Integer
  EFI	==> Expression Fraction Integer
  SOCDF	==> Segment OrderedCompletion DoubleFloat
  VEDF	==> Vector Expression DoubleFloat
  VEF	==> Vector Expression Float
  EF	==> Expression Float
  VDF	==> Vector DoubleFloat
  VMF	==> Vector MachineFloat
  MF	==> MachineFloat
  ODEA	==> Record(xinit:DF,xend:DF,fn:VEDF,yinit:LDF,intvals:LDF,_
                      g:EDF,abserr:DF,relerr:DF)
  RSS	==> Record(stiffnessFactor:F,stabilityFactor:F)
  INT	==> Integer
  EF2	==> ExpressionFunctions2

  import d02AgentsPackage, NagOrdinaryDifferentialEquationsPackage
  import AttributeButtons

  accuracyCF(ode:ODEA):F ==
    b := getButtonValue("d02bhf","accuracy")$AttributeButtons
    accuracyIntensityValue := combineFeatureCompatibility(b,accuracyIF(ode))
    accuracyIntensityValue > 0.999 => 0$F
    0.8*exp(-((10*accuracyIntensityValue)**3)$F/266)$F

  stiffnessCF(stiffnessIntensityValue:F):F ==
    b := getButtonValue("d02bhf","stiffness")$AttributeButtons
    0.5*exp(-(2*combineFeatureCompatibility(b,stiffnessIntensityValue))**2)$F

  stabilityCF(stabilityIntensityValue:F):F ==
    b := getButtonValue("d02bhf","stability")$AttributeButtons
    0.5 * cos(combineFeatureCompatibility(b,stabilityIntensityValue))$F

  expenseOfEvaluationCF(ode:ODEA):F ==
    b := getButtonValue("d02bhf","expense")$AttributeButtons
    expenseOfEvaluationIntensityValue := 
      combineFeatureCompatibility(b,expenseOfEvaluationIF(ode))
    0.35+0.2*exp(-(2.0*expenseOfEvaluationIntensityValue)**3)$F
    
  measure(R:RoutinesTable,args:ODEA) ==
    m := getMeasure(R,d02bhf :: Symbol)$RoutinesTable
    ssf := stiffnessAndStabilityOfODEIF args
    m := combineFeatureCompatibility(m,[accuracyCF(args),
            stiffnessCF(ssf.stiffnessFactor),
              expenseOfEvaluationCF(args),
                stabilityCF(ssf.stabilityFactor)])
    [m,"Runge-Kutta Merson method"]

  ODESolve(ode:ODEA) ==
    irelab := 0$INT
    if positive?(a := ode.abserr) then 
      inc(irelab)$INT
    if positive?(r := ode.relerr) then
      inc(irelab)$INT
    if positive?(a+r) then
      tol := max(a,r)$DF
    else
      tol:DF := float(1,-4,10)$DF
    asp7:Union(fn:FileName,fp:Asp7(FCN)) := 
      [retract(e:VEF := vedf2vef(ode.fn)$ExpertSystemToolsPackage)$Asp7(FCN)]
    asp9:Union(fn:FileName,fp:Asp9(G)) := 
      [retract(edf2ef(ode.g)$ExpertSystemToolsPackage)$Asp9(G)]
    d02bhf(ode.xend,# e,irelab,0$DF,ode.xinit,matrix([ode.yinit])$MDF,
             tol,-1,asp9,asp7)

@
\section{domain D02CJFA d02cjfAnnaType}
<<domain D02CJFA d02cjfAnnaType>>=
)abbrev domain D02CJFA d02cjfAnnaType
++ Author: Brian Dupee
++ Date Created: February 1995
++ Date Last Updated: January 1996
++ Basic Operations: 
++ Description:
++ \axiomType{d02cjfAnnaType} is a domain of 
++ \axiomType{OrdinaryDifferentialEquationsInitialValueProblemSolverCategory}
++ for the NAG routine D02CJF, a ODE routine which uses an  
++ Adams-Moulton-Bashworth method to solve a system of differential 
++ equations.  The function \axiomFun{measure} measures the
++ usefulness of the routine D02CJF for the given problem.  The 
++ function \axiomFun{ODESolve} performs the integration by using 
++ \axiomType{NagOrdinaryDifferentialEquationsPackage}.

d02cjfAnnaType():OrdinaryDifferentialEquationsSolverCategory == Result add  -- Adams
  EDF	==> Expression DoubleFloat
  LDF	==> List DoubleFloat
  MDF	==> Matrix DoubleFloat
  DF	==> DoubleFloat
  F	==> Float
  FI	==> Fraction Integer
  EFI	==> Expression Fraction Integer
  SOCDF	==> Segment OrderedCompletion DoubleFloat
  VEDF	==> Vector Expression DoubleFloat
  VEF	==> Vector Expression Float
  EF	==> Expression Float
  VDF	==> Vector DoubleFloat
  VMF	==> Vector MachineFloat
  MF	==> MachineFloat
  ODEA	==> Record(xinit:DF,xend:DF,fn:VEDF,yinit:LDF,intvals:LDF,_
                      g:EDF,abserr:DF,relerr:DF)
  RSS	==> Record(stiffnessFactor:F,stabilityFactor:F)
  INT	==> Integer
  EF2	==> ExpressionFunctions2

  import d02AgentsPackage, NagOrdinaryDifferentialEquationsPackage

  accuracyCF(ode:ODEA):F ==
    b := getButtonValue("d02cjf","accuracy")$AttributeButtons
    accuracyIntensityValue := combineFeatureCompatibility(b,accuracyIF(ode))
    accuracyIntensityValue > 0.9999 => 0$F
    0.6*(cos(accuracyIntensityValue*(pi()$F)/2)$F)**0.755

  stiffnessCF(ode:ODEA):F ==
    b := getButtonValue("d02cjf","stiffness")$AttributeButtons
    ssf := stiffnessAndStabilityOfODEIF ode
    stiffnessIntensityValue := 
      combineFeatureCompatibility(b,ssf.stiffnessFactor)
    0.5*exp(-(1.1*stiffnessIntensityValue)**3)$F

  measure(R:RoutinesTable,args:ODEA) ==
    m := getMeasure(R,d02cjf :: Symbol)$RoutinesTable
    m := combineFeatureCompatibility(m,[accuracyCF(args), stiffnessCF(args)])
    [m,"Adams method"]

  ODESolve(ode:ODEA) ==
    i:LDF := ode.intvals
    if empty?(i) then
      i := [ode.xend]
    M := inc(# i)$INT
    if positive?((a := ode.abserr)*(r := ode.relerr))$DF then
      ire:String := "D"
    else 
      if positive?(a) then
        ire:String := "A"
      else 
        ire:String := "R"
    tol := max(a,r)$DF
    asp7:Union(fn:FileName,fp:Asp7(FCN)) := 
      [retract(e:VEF := vedf2vef(ode.fn)$ExpertSystemToolsPackage)$Asp7(FCN)]
    asp8:Union(fn:FileName,fp:Asp8(OUTPUT)) := 
      [coerce(ldf2vmf(i)$ExpertSystemToolsPackage)$Asp8(OUTPUT)]
    asp9:Union(fn:FileName,fp:Asp9(G)) := 
      [retract(edf2ef(ode.g)$ExpertSystemToolsPackage)$Asp9(G)]
    d02cjf(ode.xend,M,# e,tol,ire,ode.xinit,matrix([ode.yinit])$MDF,
             -1,asp9,asp7,asp8)

@
\section{domain D02EJFA d02ejfAnnaType}
<<domain D02EJFA d02ejfAnnaType>>=
)abbrev domain D02EJFA d02ejfAnnaType
++ Author: Brian Dupee
++ Date Created: February 1995
++ Date Last Updated: January 1996
++ Basic Operations: 
++ Description:
++ \axiomType{d02ejfAnnaType} is a domain of 
++ \axiomType{OrdinaryDifferentialEquationsInitialValueProblemSolverCategory}
++ for the NAG routine D02EJF, a ODE routine which uses a backward 
++ differentiation formulae method to handle a stiff system 
++ of differential equations.  The function \axiomFun{measure} measures
++ the usefulness of the routine D02EJF for the given problem.  The 
++ function \axiomFun{ODESolve} performs the integration by using 
++ \axiomType{NagOrdinaryDifferentialEquationsPackage}.

d02ejfAnnaType():OrdinaryDifferentialEquationsSolverCategory == Result add -- BDF "Stiff"
  EDF	==> Expression DoubleFloat
  LDF	==> List DoubleFloat
  MDF	==> Matrix DoubleFloat
  DF	==> DoubleFloat
  F	==> Float
  FI	==> Fraction Integer
  EFI	==> Expression Fraction Integer
  SOCDF	==> Segment OrderedCompletion DoubleFloat
  VEDF	==> Vector Expression DoubleFloat
  VEF	==> Vector Expression Float
  EF	==> Expression Float
  VDF	==> Vector DoubleFloat
  VMF	==> Vector MachineFloat
  MF	==> MachineFloat
  ODEA	==> Record(xinit:DF,xend:DF,fn:VEDF,yinit:LDF,intvals:LDF,_
                      g:EDF,abserr:DF,relerr:DF)
  RSS	==> Record(stiffnessFactor:F,stabilityFactor:F)
  INT	==> Integer
  EF2	==> ExpressionFunctions2

  import d02AgentsPackage, NagOrdinaryDifferentialEquationsPackage

  accuracyCF(ode:ODEA):F ==
    b := getButtonValue("d02ejf","accuracy")$AttributeButtons
    accuracyIntensityValue :=  combineFeatureCompatibility(b,accuracyIF(ode))
    accuracyIntensityValue > 0.999 => 0$F
    0.5*exp(-((10*accuracyIntensityValue)**3)$F/250)$F

  intermediateResultsCF(ode:ODEA):F ==
    intermediateResultsIntensityValue := intermediateResultsIF(ode)
    i := 0.5 * exp(-(intermediateResultsIntensityValue/1.649)**3)$F
    a := accuracyCF(ode)
    i+(0.5-i)*(0.5-a)

  stabilityCF(ode:ODEA):F ==
    b := getButtonValue("d02ejf","stability")$AttributeButtons
    ssf := stiffnessAndStabilityOfODEIF ode
    stabilityIntensityValue := 
      combineFeatureCompatibility(b,ssf.stabilityFactor)
    0.68 - 0.5 * exp(-(stabilityIntensityValue)**3)$F

  expenseOfEvaluationCF(ode:ODEA):F ==
    b := getButtonValue("d02ejf","expense")$AttributeButtons
    expenseOfEvaluationIntensityValue := 
      combineFeatureCompatibility(b,expenseOfEvaluationIF(ode))
    0.5 * exp(-(1.7*expenseOfEvaluationIntensityValue)**3)$F
    
  systemSizeCF(args:ODEA):F ==
    (1$F - systemSizeIF(args))/2.0

  measure(R:RoutinesTable,args:ODEA) ==
    arg := copy args
    m := getMeasure(R,d02ejf :: Symbol)$RoutinesTable
    m := combineFeatureCompatibility(m,[intermediateResultsCF(arg),
           accuracyCF(arg),
             systemSizeCF(arg),
               expenseOfEvaluationCF(arg),
                 stabilityCF(arg)])
    [m,"BDF method for Stiff Systems"]

  ODESolve(ode:ODEA) ==
    i:LDF := ode.intvals
    m := inc(# i)$INT
    if positive?((a := ode.abserr)*(r := ode.relerr))$DF then
      ire:String := "D"
    else 
      if positive?(a) then
        ire:String := "A"
      else 
        ire:String := "R"
    if positive?(a+r)$DF then
      tol := max(a,r)$DF
    else 
      tol := float(1,-4,10)$DF
    asp7:Union(fn:FileName,fp:Asp7(FCN)) := 
      [retract(e:VEF := vedf2vef(ode.fn)$ExpertSystemToolsPackage)$Asp7(FCN)]
    asp31:Union(fn:FileName,fp:Asp31(PEDERV)) := 
      [retract(e)$Asp31(PEDERV)]
    asp8:Union(fn:FileName,fp:Asp8(OUTPUT)) := 
      [coerce(ldf2vmf(i)$ExpertSystemToolsPackage)$Asp8(OUTPUT)]
    asp9:Union(fn:FileName,fp:Asp9(G)) :=
      [retract(edf2ef(ode.g)$ExpertSystemToolsPackage)$Asp9(G)]
    n:INT := # ode.yinit
    iw:INT := (12+n)*n+50
    ans := d02ejf(ode.xend,m,n,ire,iw,ode.xinit,matrix([ode.yinit])$MDF,
             tol,-1,asp9,asp7,asp31,asp8)

@
\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 D02BBFA d02bbfAnnaType>>
<<domain D02BHFA d02bhfAnnaType>>
<<domain D02CJFA d02cjfAnnaType>>
<<domain D02EJFA d02ejfAnnaType>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}