/usr/share/axiom-20170501/src/algebra/D02EJFA.spad is in axiom-source 20170501-3.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 | )abbrev domain D02EJFA d02ejfAnnaType
++ Author: Brian Dupee
++ Date Created: February 1995
++ Date Last Updated: January 1996
++ References:
++ Dupe95 Using Computer Algebra to Choose and Apply Numerical Routines
++ Dewa92 Using Computer Algebra to Select Numerical Algorithms
++ 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() : SIG == CODE where
-- 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
SIG ==> OrdinaryDifferentialEquationsSolverCategory
CODE ==> Result add
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)
|