/usr/lib/open-axiom/input/lodo.input is in open-axiom-test 1.5.0~svn3056+ds-1.
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 111 112 113 | --Copyright The Numerical Algorithms Group Limited 1991.
---------------------------------- lodo.input -----------------------------
-- LODO2(M,A) is the domain of linear ordinary differential operators over
-- an A-module M, where A is a differential ring. This includes the
-- cases of operators which are polynomials in D acting upon scalars or
-- vectors depending on a single variable. The coefficients of the
-- operator polynomials can be integers, rational functions, matrices
-- or elements of other domains.
------------------------------------------------------------------------
-- Differential operators with constant coefficients
------------------------------------------------------------------------
)clear all
RN:=FRAC INT
Dx: LODO2(RN, UP(x,RN))
Dx := D() -- definition of an operator
a := Dx + 1
b := a + 1/2*Dx**2 - 1/2
p: UP(x,RN) := 4*x**2 + 2/3 -- something to work on
a p -- application of an operator to a polynomial
(a*b) p = a b p -- multiplication is defined by this identity
c := (1/9)*b*(a + b)**2 -- exponentiation follows from multiplication
(a**2 - 3/4*b + c) (p + 1) -- general application of operator expressions
------------------------------------------------------------------------
-- Differential operators with rational function coefficients
------------------------------------------------------------------------
)clear all
RFZ := FRAC UP(x,INT)
(Dx, a, b): LODO1 RFZ
Dx := D()
b := 3*x**2*Dx**2 + 2*Dx + 1/x
a := b*(5*x*Dx + 7)
p: RFZ := x**2 + 1/x**2
(a*b - b*a) p -- operator multiplication is not commutative
-- When the coefficients of the operator polynomials come from a field
-- it is possible to define left and right division of the operators.
-- This allows the computation of left and right gcd's via remainder
-- sequences, and also the computation of left and right lcm's.
leftDivide(a,b) -- result is the quotient/remainder pair
a - (b * %.quotient + %.remainder)
rightDivide(a,b)
a - (%.quotient * b + %.remainder)
-- A GCD doesn't necessarily divide a and b on both sides.
e := leftGcd(a,b)
leftRemainder(a, e) -- remainder from left division
rightRemainder(a, e) -- remainder from right division
-- An LCM is not necessarily divisible from both sides.
f := rightLcm(a,b)
leftRemainder(f, b)
rightRemainder(f, b) -- the remainder is non-zero
------------------------------------------------------------------------
--
-- Problem: find the first few coefficients of exp(x)/x**i in
-- Dop phi
-- where
-- Dop := D**3 + G/x**2 * D + H/x**3 - 1
-- phi := sum(s[i]*exp(x)/x**i, i = 0..)
------------------------------------------------------------------------
)clear all
Dx: LODO(EXPR INT, f +-> D(f, x))
Dx := D()
Dop:= Dx**3 + G/x**2*Dx + H/x**3 - 1
n == 3
phi == reduce(+,[subscript(s,[i])*exp(x)/x**i for i in 0..n])
phi1 == Dop(phi) / exp x
phi2 == phi1 *x**(n+3)
phi3 == retract(phi2)@(POLY INT)
pans == phi3 ::UP(x,POLY INT)
pans1 == [coefficient(pans, (n+3-i) :: NNI) for i in 2..n+1]
leq == solve(pans1,[subscript(s,[i]) for i in 1..n])
leq
n==4
leq
n==7
leq
------------------------------------------------------------------------
-- Differential operators with matrix coefficients acting on vectors.
------------------------------------------------------------------------
)clear all
PZ := UP(x,INT); Vect := DPMM(3, PZ, SQMATRIX(3,PZ), PZ);
Modo := LODO2(SQMATRIX(3,PZ), Vect)
p := directProduct([3*x**2 + 1, 2*x, 7*x**3 + 2*x]::(VECTOR(PZ)))@Vect
m := [[x**2, 1, 0], [1, x**4, 0], [0, 0, 4*x**2]]::(SQMATRIX(3,PZ))
-- Vect is a left SM(3,PZ)-module
q: Vect := m * p
-- Operator combination and application
Dx: Modo := D()
a: Modo := 1*Dx + m
b: Modo := m*Dx + 1
a*b
a p
b p
(a+b) (p + q)
|