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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
|
\documentclass{article}
\usepackage{axiom}
\begin{document}
\title{\$SPAD/src/input lodo.input}
\author{The Axiom Team}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{License}
<<license>>=
--Copyright The Numerical Algorithms Group Limited 1991.
@
<<*>>=
<<license>>
---------------------------------- 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)
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}
|
