$SPAD/src/input lodesys.input

Timothy Daly December 3, 2016

Abstract

Contents

2

1 License

— license —

--Copyright The Numerical Algorithms Group Limited 1991.

———-

— * —

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

———-

There are 2 different ways to input a homogeneous 1st order system of linear ordinary differential equations of the form dy/dt = M y where y is a vector of unknown functions of t. the first is simply solve(M, t) which will be understood to be a differential system:

— * —

--S 1 of 13
M := matrix [[ 1+4*t,  -5*t,   7*t,  -8*t,   8*t,  -6*t],_
             [ -10*t, 1+9*t, -14*t,  16*t, -16*t,  12*t],_
             [  -5*t,   5*t, 1-8*t,   8*t,  -8*t,   6*t],_
             [  10*t, -10*t,  14*t,1-17*t,  16*t, -12*t],_
             [   5*t,  -5*t,   7*t,  -8*t, 1+7*t,  -6*t],_
             [  -5*t,   5*t,  -7*t,   8*t,  -8*t, 1+5*t]]
--R
--R
--R        +4t + 1   - 5t      7t       - 8t       8t     - 6t +
--R        |                                                   |
--R        |- 10t   9t + 1   - 14t       16t     - 16t    12t  |
--R        |                                                   |
--R        | - 5t     5t    - 8t + 1     8t       - 8t     6t  |
--R   (1)  |                                                   |
--R        | 10t    - 10t     14t     - 17t + 1   16t    - 12t |
--R        |                                                   |
--R        |  5t     - 5t      7t       - 8t     7t + 1   - 6t |
--R        |                                                   |
--R        + - 5t     5t      - 7t       8t       - 8t   5t + 1+
--R                                            Type: Matrix(Polynomial(Integer))
--E 1

———-

The original system in Barkatou’s AAECC paper is t^{2}dy/dt = M ∗ y — * —

--S 2 of 13
sol := solve(inv(t**2) * M, t)
--R
--R
--R   (2)
--R           1          1         1        1       1         1
--R         - -        - -       - -      - -     - -       - -
--R      5    t     5    t    5    t   5    t  5    t    5    t
--R   [[t %e   ,- 2t %e   ,- t %e   ,2t %e   ,t %e   ,- t %e   ],
--R         1        1       1      1     1       1           1     1
--R       - -      - -     - -    - -   - -     - -         - -   - -
--R         t        t       t      t     t       t           t     t
--R     %e      4%e      %e    2%e    %e      %e         7%e    %e
--R    [-----,- ------,- -----,------,-----,- -----], [0,------,-----,0,0,0],
--R       t       5t       t      t     t       t          5t     t
--R              1       1               1         1               1           1
--R            - -     - -             - -       - -             - -         - -
--R              t       t               t         t               t           t
--R         8%e      %e             8%e        %e             6%e          %e
--R    [0,- ------,0,-----,0,0], [0,------,0,0,-----,0], [0,- ------,0,0,0,-----]]
--R           5t       t              5t         t              5t           t
--R                           Type: Union(List(Vector(Expression(Integer))),...)
--E 2

———-

Verify the solutions — * —

--S 3 of 13
[t**2 * map(h +-> D(h, t), v) - M * v for v in sol]
--R
--R
--R   (3)
--R   [[0,0,0,0,0,0], [0,0,0,0,0,0], [0,0,0,0,0,0], [0,0,0,0,0,0], [0,0,0,0,0,0],
--R    [0,0,0,0,0,0]]
--R                                      Type: List(Vector(Expression(Integer)))
--E 3

———-

The second way is to type each equation using a separate operator for each unknown: — * —

--S 4 of 13
x := operator x
--R
--R
--R   (4)  x
--R                                                          Type: BasicOperator
--E 4

--S 5 of 13
y := operator y
--R
--R
--R   (5)  y
--R                                                          Type: BasicOperator
--E 5
--S 6 of 13
sys := [D(x t, t) = x t + sqrt 3 * y t, D(y t, t) = sqrt 3 * x t - y t]
--R
--R
--R          ,      +-+             ,               +-+
--R   (6)  [x (t)= \|3 y(t) + x(t),y (t)= - y(t) + \|3 x(t)]
--R
--R                                    Type: List(Equation(Expression(Integer)))
--E 6
--S 7 of 13
solve(sys, [x, y], t).basis
--R
--R
--R                 2t               - 2t
--R            2t %e       - 2t   3%e
--R   (7)  [[%e  ,----],[%e    ,- -------]]
--R                +-+               +-+
--R               \|3               \|3
--R                                      Type: List(Vector(Expression(Integer)))
--E 7

———-

Similarly there are 2 different ways to input the inhomogeneous system dy/dt = M y + v where v is a given vector of functions.

The first is solve(M, v, t): — * —

--S 8 of 13
v := vector [1, (-29*t + 19)/5, -1, t + 1, - 2*t + 3, -1]
--R
--R
--R             29     19
--R   (8)  [1,- -- t + --,- 1,t + 1,- 2t + 3,- 1]
--R              5      5
--R                                  Type: Vector(Polynomial(Fraction(Integer)))
--E 8

———- Get a particular solution to t^{2}dy/dt = My + v

— * —

--S 9 of 13
solp := solve(inv(t**2) * M, inv(t**2) * v, t).particular
--R
--R
--R               19
--R   (9)  [- 1,- --,1,- 1,- 3,1]
--R                5
--R                                            Type: Vector(Expression(Integer))
--E 9

———-

Verify the particular solution — * —

--S 10 of 13
t**2 * map(h +-> D(h, t), solp) - M * solp - v
--R
--R
--R   (10)  [0,0,0,0,0,0]
--R                                            Type: Vector(Expression(Integer))
--E 10

———-

The second way is by listing the equations:

— * —

--S 11 of 13
z := operator z
--R
--R
--R   (11)  z
--R                                                          Type: BasicOperator
--E 11
--S 12 of 13
sys := [D(x t, t) = y t + z t + t, D(y t, t) = x t + z t, D(z t, t) = x t + y t]
--R
--R
--R           ,                      ,                  ,
--R   (12)  [x (t)= z(t) + y(t) + t,y (t)= z(t) + x(t),z (t)= y(t) + x(t)]
--R
--R                                    Type: List(Equation(Expression(Integer)))
--E 12
--S 13 of 13
solve(sys, [x, y, z], t).particular

--R
--R
--R          2t - 3 - 2t + 1 - 2t + 1
--R   (13)  [------,--------,--------]
--R             4       4        4
--R                                            Type: Vector(Expression(Integer))
--E 13
)spool
)lisp (bye)

———-

References^{[1] nothing}

8