$SPAD/src/input lodo.input

Timothy Daly

December 3, 2016

Abstract

LODO2(M,A) is the domain of linear ordinary differential operators over an Amodule 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.

Contents

2

1 License

— license —

--Copyright The Numerical Algorithms Group Limited 1991.

———-

— * —

)set break resume
)spool lodo.output
)set message test on
)set message auto off
)clear all

———-

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.

2 Differential operators with constant coefficients

— * —

--S 1 of 55
RN:=FRAC INT
--R
--R
--R   (1)  Fraction(Integer)
--R                                                                 Type: Domain
--E 1
--S 2 of 55
Dx: LODO2(RN, UP(x,RN))
--R
--R                                                                   Type: Void
--E 2

———-

Definition of an operator — * —

--S 3 of 55
Dx := D()

--R
--R
--R   (3)  D
--IType: LinearOrdinaryDifferentialOperator2(Fraction(Integer),...
--E 3
--S 4 of 55
a  := Dx  + 1
--R
--R
--R   (4)  D + 1
--IType: LinearOrdinaryDifferentialOperator2(Fraction(Integer),...
--E 4
--S 5 of 55
b  := a + 1/2*Dx**2 - 1/2
--R
--R
--R        1  2       1
--R   (5)  - D  + D + -
--R        2          2
--IType: LinearOrdinaryDifferentialOperator2(Fraction(Integer),...
--E 5

———-

Something to work on — * —

--S 6 of 55
p: UP(x,RN) := 4*x**2 + 2/3
--R
--R
--R          2   2
--R   (6)  4x  + -
--R              3
--R                              Type: UnivariatePolynomial(x,Fraction(Integer))
--E 6

———-

Application of an operator to a polynomial — * —

--S 7 of 55
a p
--R
--R
--R          2        2
--R   (7)  4x  + 8x + -
--R                   3

--R                              Type: UnivariatePolynomial(x,Fraction(Integer))
--E 7

———-

Multiplication is defined by this identity — * —

--S 8 of 55
(a*b) p = a b p
--R
--R
--R          2         37    2         37
--R   (8)  2x  + 12x + --= 2x  + 12x + --
--R                     3               3
--R                    Type: Equation(UnivariatePolynomial(x,Fraction(Integer)))
--E 8

———-

Exponentiation follows from multiplication — * —

--S 9 of 55
c := (1/9)*b*(a + b)**2
--R
--R
--R         1  6    5  5   13  4   19  3   79  2    7     1
--R   (9)  -- D  + -- D  + -- D  + -- D  + -- D  + -- D + -
--R        72      36      24      18      72      12     8
--IType: LinearOrdinaryDifferentialOperator2(Fraction(Integer),...
--E 9

———-

General application of operator expressions — * —

--S 10 of 55
(a**2 - 3/4*b + c) (p + 1)
--R
--R
--R           2   44     541
--R   (10)  3x  + -- x + ---
--R                3      36
--R                              Type: UnivariatePolynomial(x,Fraction(Integer))
--E 10

———-

3 Differential operators with rational function coefficients

— * —

)clear all
--S 11 of 55
RFZ := FRAC UP(x,INT)
--R
--R
--R   (1)  Fraction(UnivariatePolynomial(x,Integer))
--R                                                                 Type: Domain
--E 11
--S 12 of 55
(Dx, a, b): LODO1 RFZ
--R
--R                                                                   Type: Void
--E 12
--S 13 of 55
Dx := D()
--R
--R
--R   (3)  D
--IType: LinearOrdinaryDifferentialOperator1(Fraction(...
--E 13
--S 14 of 55
b := 3*x**2*Dx**2 + 2*Dx + 1/x
--R
--R
--R          2 2        1
--R   (4)  3x D  + 2D + -
--R                     x
--IType: LinearOrdinaryDifferentialOperator1(Fraction(...
--E 14
--S 15 of 55
a := b*(5*x*Dx + 7)
--R
--R
--R           3 3       2        2         7
--R   (5)  15x D  + (51x  + 10x)D  + 29D + -
--R                                        x
--IType: LinearOrdinaryDifferentialOperator1(Fraction(...
--E 15

--S 16 of 55
p: RFZ := x**2 + 1/x**2
--R
--R
--R         4
--R        x  + 1
--R   (6)  ------
--R           2
--R          x
--R                              Type: Fraction(UnivariatePolynomial(x,Integer))
--E 16

———-

Operator multiplication is not commutative — * —

--S 17 of 55
(a*b - b*a) p
--R
--R
--R             4
--R        - 75x  + 540x - 75
--R   (7)  ------------------
--R                 4
--R                x
--R                              Type: Fraction(UnivariatePolynomial(x,Integer))
--E 17

———-

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.

The result is the quotient/remainder pair — * —

--S 18 of 55
leftDivide(a,b)
--R
--R
--R   (8)  [quotient= 5x D + 7,remainder= 0]
--IType: Record(quotient: LinearOrdinaryDifferentialOperator1(...
--E 18
--S 19 of 55
a - (b * %.quotient + %.remainder)
--R
--R
--R   (9)  0

--IType: LinearOrdinaryDifferentialOperator1(Fraction(...
--E 19
--S 20 of 55
rightDivide(a,b)
--R
--R
--R                                              5
--R   (10)  [quotient= 5x D + 7,remainder= 10D + -]
--R                                              x
--IType: Record(quotient: LinearOrdinaryDifferentialOperator1(...
--E 20
--S 21 of 55
a - (%.quotient * b + %.remainder)
--R
--R
--R   (11)  0
--IType: LinearOrdinaryDifferentialOperator1(Fraction(...
--E 21

———-

A GCD doesn’t necessarily divide a and b on both sides.

— * —

--S 22 of 55
e := leftGcd(a,b)
--R
--R
--R           2 2        1
--R   (12)  3x D  + 2D + -
--R                      x
--IType: LinearOrdinaryDifferentialOperator1(Fraction(UnivariatePolynomial(x,Integer)))
--E 22

———-

Remainder from left division — * —

--S 23 of 55
leftRemainder(a, e)
--R
--R
--R   (13)  0
--IType: LinearOrdinaryDifferentialOperator1(Fraction(...
--E 23

———- Remainder from right division

— * —

--S 24 of 55
rightRemainder(a, e)
--R
--R
--R               5
--R   (14)  10D + -
--R               x
--IType: LinearOrdinaryDifferentialOperator1(Fraction(...
--E 24

———-

An LCM is not necessarily divisible from both sides. — * —

--S 25 of 55
f := rightLcm(a,b)
--R
--R
--R            3 3       2        2         7
--R   (15)  15x D  + (51x  + 10x)D  + 29D + -
--R                                         x
--IType: LinearOrdinaryDifferentialOperator1(Fraction(...
--E 25
--S 26 of 55
leftRemainder(f, b)
--R
--R
--R   (16)  0
--IType: LinearOrdinaryDifferentialOperator1(Fraction(...
--E 26

———-

The remainder is non-zero — * —

--S 27 of 55
rightRemainder(f, b)
--R
--R
--R               5
--R   (17)  10D + -
--R               x
--IType: LinearOrdinaryDifferentialOperator1(Fraction(...
--E 27

———-

4 Problem 1

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
--S 28 of 55
Dx: LODO(EXPR INT, f +-> D(f, x))
--R
--R                                                                   Type: Void
--E 28
--S 29 of 55
Dx := D()
--R
--R
--R   (2)  D
--IType: LinearOrdinaryDifferentialOperator(Expression(Integer),...
--E 29
--S 30 of 55
Dop:= Dx**3 + G/x**2*Dx + H/x**3 - 1
--R
--R
--R                       3
--R         3    G     - x  + H
--R   (3)  D  + -- D + --------
--R              2         3
--R             x         x
--IType: LinearOrdinaryDifferentialOperator(Expression(Integer),...
--E 30
--S 31 of 55
n == 3
--R
--R                                                                   Type: Void
--E 31
--S 32 of 55
phi == reduce(+,[subscript(s,[i])*exp(x)/x**i for i in 0..n])
--R
--R                                                                   Type: Void
--E 32
--S 33 of 55

phi1 ==  Dop(phi) / exp x
--R
--R                                                                   Type: Void
--E 33
--S 34 of 55
phi2 == phi1 *x**(n+3)
--R
--R                                                                   Type: Void
--E 34
--S 35 of 55
phi3 == retract(phi2)@(POLY INT)
--R
--R                                                                   Type: Void
--E 35
--S 36 of 55
pans == phi3 ::UP(x,POLY INT)
--R
--R                                                                   Type: Void
--E 36
--S 37 of 55
pans1 == [coefficient(pans, (n+3-i) :: NNI) for i in 2..n+1]
--R
--R                                                                   Type: Void
--E 37
--S 38 of 55
leq == solve(pans1,[subscript(s,[i]) for i in 1..n])
--R
--R                                                                   Type: Void
--E 38
--S 39 of 55
leq
--R
--R   Compiling body of rule n to compute value of type PositiveInteger
--R   Compiling body of rule phi to compute value of type Expression(
--R      Integer)
--R   Compiling body of rule phi1 to compute value of type Expression(
--R      Integer)
--R   Compiling body of rule phi2 to compute value of type Expression(
--R      Integer)
--R   Compiling body of rule phi3 to compute value of type Polynomial(
--R      Integer)
--R   Compiling body of rule pans to compute value of type
--R      UnivariatePolynomial(x,Polynomial(Integer))
--R   Compiling body of rule pans1 to compute value of type List(

--R      Polynomial(Integer))
--R   Compiling body of rule leq to compute value of type List(List(
--R      Equation(Fraction(Polynomial(Integer)))))
--I   Compiling function G3586 with type Integer -> Boolean
--R
--R   (12)
--R                           2                                3        2
--R         s G     3s H + s G  + 6s G     (9s G + 54s )H + s G  + 18s G  + 72s G
--R          0        0     0       0         0       0      0        0        0
--R   [[s = ---,s = ------------------,s = --------------------------------------]]
--R      1   3   2          18          3                    162
--R                    Type: List(List(Equation(Fraction(Polynomial(Integer)))))
--E 39
--S 40 of 55
n==4
--R
--R   Compiled code for n has been cleared.
--R   Compiled code for leq has been cleared.
--R   Compiled code for pans1 has been cleared.
--R   Compiled code for phi2 has been cleared.
--R   Compiled code for phi has been cleared.
--R   Compiled code for phi3 has been cleared.
--R   Compiled code for phi1 has been cleared.
--R   Compiled code for pans has been cleared.
--R   1 old definition(s) deleted for function or rule n
--R                                                                   Type: Void
--E 40
--S 41 of 55
leq
--R
--R   Compiling body of rule n to compute value of type PositiveInteger
--R   Compiling body of rule phi to compute value of type Expression(
--R      Integer)
--R   Compiling body of rule phi1 to compute value of type Expression(
--R      Integer)
--R   Compiling body of rule phi2 to compute value of type Expression(
--R      Integer)
--R   Compiling body of rule phi3 to compute value of type Polynomial(
--R      Integer)
--R   Compiling body of rule pans to compute value of type
--R      UnivariatePolynomial(x,Polynomial(Integer))
--R   Compiling body of rule pans1 to compute value of type List(
--R      Polynomial(Integer))
--R   Compiling body of rule leq to compute value of type List(List(
--R      Equation(Fraction(Polynomial(Integer)))))
--R
--R   (14)
--R   [

--R                             2
--R          s G      3s H + s G  + 6s G
--R           0         0     0       0
--R     [s = ---, s = ------------------,
--R       1   3    2          18
--R                              3        2
--R          (9s G + 54s )H + s G  + 18s G  + 72s G
--R             0       0      0        0        0
--R      s = --------------------------------------,
--R       3                    162
--R
--R       s  =
--R        4
--R                  2         2                          4        3         2
--R             27s H  + (18s G  + 378s G + 1296s )H + s G  + 36s G  + 396s G
--R                0         0         0         0      0        0         0
--R           +
--R             1296s G
--R                  0
--R        /
--R           1944
--R       ]
--R     ]
--R                    Type: List(List(Equation(Fraction(Polynomial(Integer)))))
--E 41
--S 42 of 55
n==7
--R
--R   Compiled code for n has been cleared.
--R   Compiled code for leq has been cleared.
--R   Compiled code for pans1 has been cleared.
--R   Compiled code for phi2 has been cleared.
--R   Compiled code for phi has been cleared.
--R   Compiled code for phi3 has been cleared.
--R   Compiled code for phi1 has been cleared.
--R   Compiled code for pans has been cleared.
--R   1 old definition(s) deleted for function or rule n
--R                                                                   Type: Void
--E 42
--S 43 of 55
leq
--R
--R   Compiling body of rule n to compute value of type PositiveInteger
--R   Compiling body of rule phi to compute value of type Expression(
--R      Integer)
--R   Compiling body of rule phi1 to compute value of type Expression(
--R      Integer)
--R   Compiling body of rule phi2 to compute value of type Expression(

--R      Integer)
--R   Compiling body of rule phi3 to compute value of type Polynomial(
--R      Integer)
--R   Compiling body of rule pans to compute value of type
--R      UnivariatePolynomial(x,Polynomial(Integer))
--R   Compiling body of rule pans1 to compute value of type List(
--R      Polynomial(Integer))
--R   Compiling body of rule leq to compute value of type List(List(
--R      Equation(Fraction(Polynomial(Integer)))))
--R
--R   (16)
--R   [
--R                             2
--R          s G      3s H + s G  + 6s G
--R           0         0     0       0
--R     [s = ---, s = ------------------,
--R       1   3    2          18
--R                              3        2
--R          (9s G + 54s )H + s G  + 18s G  + 72s G
--R             0       0      0        0        0
--R      s = --------------------------------------,
--R       3                    162
--R
--R       s  =
--R        4
--R                  2         2                          4        3         2
--R             27s H  + (18s G  + 378s G + 1296s )H + s G  + 36s G  + 396s G
--R                0         0         0         0      0        0         0
--R           +
--R             1296s G
--R                  0
--R        /
--R           1944
--R       ,
--R
--R       s  =
--R        5
--R                               2         3          2
--R             (135s G + 2268s )H  + (30s G  + 1350s G  + 16416s G + 38880s )H
--R                  0         0          0          0           0          0
--R           +
--R                5        4          3          2
--R             s G  + 60s G  + 1188s G  + 9504s G  + 25920s G
--R              0        0          0          0           0
--R        /
--R           29160
--R       ,
--R
--R       s  =
--R        6

--R                   3          2                        2
--R             405s H  + (405s G  + 18468s G + 174960s )H
--R                 0          0           0           0
--R           +
--R                   4          3           2                                6
--R             (45s G  + 3510s G  + 88776s G  + 777600s G + 1166400s )H + s G
--R                 0          0           0            0            0      0
--R           +
--R                  5          4           3           2
--R             90s G  + 2628s G  + 27864s G  + 90720s G
--R                0          0           0           0
--R        /
--R           524880
--R       ,
--R
--R       s  =
--R        7
--R                                 3
--R             (2835s G + 91854s )H
--R                   0          0
--R           +
--R                    3           2                            2
--R             (945s G  + 81648s G  + 2082996s G + 14171760s )H
--R                  0           0             0             0
--R           +
--R                   5          4            3             2
--R             (63s G  + 7560s G  + 317520s G  + 5554008s G  + 34058880s G)H
--R                 0          0            0             0              0
--R           +
--R                7         6          5           4             3              2
--R             s G  + 126s G  + 4788s G  + 25272s G  - 1744416s G  - 26827200s G
--R              0         0          0           0             0              0
--R           +
--R             - 97977600s G
--R                        0
--R        /
--R           11022480
--R       ]
--R     ]
--R                    Type: List(List(Equation(Fraction(Polynomial(Integer)))))
--E 43

———-

5 Differential operators with matrix coefficients acting on vectors

— * —

)clear all
--S 44 of 55
PZ := UP(x,INT); Vect := DPMM(3, PZ, SQMATRIX(3,PZ), PZ);
--R
--R
--R                                                                 Type: Domain
--E 44
--S 45 of 55
Modo := LODO2(SQMATRIX(3,PZ), Vect);
--R
--R
--R                                                                 Type: Domain
--E 45
--S 46 of 55
p := directProduct([3*x**2 + 1, 2*x, 7*x**3 + 2*x]::(VECTOR(PZ)))@Vect
--R
--R
--R           2          3
--R   (3)  [3x  + 1,2x,7x  + 2x]
--IType: DirectProductMatrixModule(3,UnivariatePolynomial(x,Integer),...
--E 46
--S 47 of 55
m := [[x**2, 1, 0], [1, x**4, 0], [0, 0, 4*x**2]]::(SQMATRIX(3,PZ))
--R
--R
--R        + 2         +
--R        |x   1    0 |
--R        |           |
--R   (4)  |     4     |
--R        |1   x    0 |
--R        |           |
--R        |          2|
--R        +0   0   4x +
--R                        Type: SquareMatrix(3,UnivariatePolynomial(x,Integer))
--E 47

———-

Vect is a left SM(3,PZ)-module — * —

--S 48 of 55
q: Vect := m * p
--R
--R

--R           4    2        5     2        5     3
--R   (5)  [3x  + x  + 2x,2x  + 3x  + 1,28x  + 8x ]
--IType: DirectProductMatrixModule(3,UnivariatePolynomial(x,Integer),...
--E 48

———-

Operator combination and application — * —

--S 49 of 55
Dx:  Modo := D()
--R
--R
--R   (6)  D
--IType: LinearOrdinaryDifferentialOperator2(SquareMatrix(3,...
--E 49
--S 50 of 55
a:   Modo := 1*Dx  + m
--R
--R
--R            + 2         +
--R            |x   1    0 |
--R            |           |
--R   (7)  D + |     4     |
--R            |1   x    0 |
--R            |           |
--R            |          2|
--R            +0   0   4x +
--IType: LinearOrdinaryDifferentialOperator2(SquareMatrix(3,...
--E 50
--S 51 of 55
b:   Modo := m*Dx  + 1
--R
--R
--R        + 2         +
--R        |x   1    0 |    +1  0  0+
--R        |           |    |       |
--R   (8)  |     4     |D + |0  1  0|
--R        |1   x    0 |    |       |
--R        |           |    +0  0  1+
--R        |          2|
--R        +0   0   4x +
--IType: LinearOrdinaryDifferentialOperator2(SquareMatrix(3,...
--E 51
--S 52 of 55
a*b

--R
--R
--R   (9)
--R   + 2         +     + 4              4    2                  +    + 2         +
--R   |x   1    0 |     |x  + 2x + 2    x  + x            0      |    |x   1    0 |
--R   |           | 2   |                                        |    |           |
--R   |     4     |D  + |   4    2     8     3                   |D + |     4     |
--R   |1   x    0 |     |  x  + x     x  + 4x  + 2        0      |    |1   x    0 |
--R   |           |     |                                        |    |           |
--R   |          2|     |                              4         |    |          2|
--R   +0   0   4x +     +     0            0        16x  + 8x + 1+    +0   0   4x +
--IType: LinearOrdinaryDifferentialOperator2(SquareMatrix(3,...
--E 52
--S 53 of 55
a p
--R
--R
--R            4    2        5     2        5     3      2
--R   (10)  [3x  + x  + 8x,2x  + 3x  + 3,28x  + 8x  + 21x  + 2]
--IType: DirectProductMatrixModule(3,UnivariatePolynomial(x,Integer),...
--E 53
--S 54 of 55
b p
--R
--R
--R            3     2       4         4     3     2
--R   (11)  [6x  + 3x  + 3,2x  + 8x,84x  + 7x  + 8x  + 2x]
--IType: DirectProductMatrixModule(3,UnivariatePolynomial(x,Integer),...
--E 54
--S 55 of 55
(a+b) (p + q)
--R
--R
--R   (12)
--R      6      5      4      3      2
--R   [3x  + 14x  + 17x  + 22x  + 10x  + 18x + 6,
--R      9      8     6      5      4      3     2
--R    2x  + 10x  + 3x  + 10x  + 16x  + 12x  + 7x  + 18x + 6,
--R        7       6      5       4      3      2
--R    112x  + 560x  + 88x  + 320x  + 23x  + 53x  + 2x + 2]
--IType: DirectProductMatrixModule(3,UnivariatePolynomial(x,Integer),...
--E 55
)spool
)lisp (bye)

———-

References^{[1] nothing}

19