/usr/share/axiom-20170501/input/tutchap67.input

--Copyright The Numerical Algorithms Group Limited 1996.
(vecA,vecB) : Vector Integer
vecA := vector [3,0,4]
vecB := vector [2,4,-2]
vecA + vecB
vecA - vecB
5*vecA
vecB*7
vecB/2
1/2 * vecB
vecA * (1/2)
dot(vecA,vecB)
sqrt dot(vecA,vecA)
magnitude vecA
direction x == 1/magnitude x * x
direction vecA
magnitude direction vecB
p := vector [u1*t,u2*t,u3*t-1/2*g*t^2]
D(p,t)
DV(s,t) == map(x+->D(x,t),s)
v := DV(p,t) -- the velocity,
a := DV(v,t) -- acceleration
j := DV(a,t) -- and jerk
)clear properties all
s := operator 's;
sSol := solve(D(s(t),t,2) - aF + r/m*D(s(t),t),s,t=0,[0,u])
u := vector [ux,uy]
aF := vector [0,-g]
function(sSol ,'s,'t)
s t
function(numerator sSol,'n,'t)
function(1/denominator sSol*n(t),'s,'t)
s t
outputGeneral 6
map(x+->eval(x,[m=1,ux=20,uy=10,r=0.1,g=9.8]),s t)
draw(curve(%.1,%.2),t=0..2)
map(x+->eval(x,[m=1,ux=20,uy=10,r=0,g=9.8]),s t)
map(x+->limit(x,r=0),s t)
map(x+->eval(x,[m=1,ux=20,uy=10,g=9.8]),%::Vector Expression Integer)
draw(curve(%.1,%.2),t=0..2)
Integer has Ring
PositiveInteger has Ring 
has(Fraction Integer,DivisionRing)
has(Integer,DivisionRing)
List has Group
List ? has Group
Field has Ring                   
Group has Ring
vecC := vector [0,x +-> x,[a,b,c]]  
matA := matrix [[x,0],[7.3,%i]]     
vecC(1)
vecC.2
vecC 3
vecD : Vector Integer := [1,2,3,4,5,6]
matB : Matrix Integer := [[1,2,3],[4,5,6]]
vector [1,2,3] :: Matrix Integer
matA(1,2)
matrix [[matA]]
matrix [[matA::SquareMatrix(2,Polynomial Complex Float)]]
matrix [[squareMatrix matA]]
matC := matrix [[a,b,c],[d,e,f]]
lvec := vector [2,3]            
rvec := vector [4,5,6]
lvec * matC
matC * rvec
lrvec := vector [1,2]
lrvec * ((matrix [[a,b],[c,d]] * lrvec) :: Matrix Polynomial Integer)
vecD := new(5,0)          
new(3,3,0)$Matrix Integer
Z3 := %;
I3 := Z3; I3(1,1) := 1; I3(2,2) := 1; I3(3,3) := 1;
I3
I3 := Z3; for k in [1,2,3] repeat I3(k,k) := 1
I3
I3 := Z3; for k in 1..3 repeat I3(k,k) := 1
I3
expand(-7/2..7/3)
for i in 11..20 repeat _
           ( print i; _
             if prime? i then messagePrint(" (That was prime.)")$OutputForm )
poly := 0;
for i in [1,2,3,4] for c in ['a,'b,'c,'d] repeat poly := poly + i*c
poly
vecC := vector [n^2 for n in 1..3]
hilbert3 := matrix [[1/(i+j) for i in 1..3] for j in 1..3]
I3 := matrix [[((m,n)+->if m=n then 1 else 0)(i,j) _
                                             for i in 1..3] for j in 1..3]
diagonalMatrix [1,1,1]
inverse hilbert3                            
% * hilbert3                                
matC := matrix [[a,b],[c,d]]
inverse matC
determinant matC
(x-a)^2 + (y-b)^2 - r^2 = 0
row1 := [1,x,y,x^2+y^2]
row2 := [eval(row1.i,[x,y],[0,1]) for i in 1..4]
row3 := [eval(row1.i,[x,y],[1,0]) for i in 1..4];
row4 := [eval(row1.i,[x,y],[-1,-1]) for i in 1..4];
determinant [row1,row2,row3,row4] = 0
solve([[1/(i+j) for i in 1..3] for j in 1..3],[3,5,7])
matD := matrix [[1,2,3,4],[2,3,4,5],[3,4,5,6],[4,5,6,7]]
c := vector [5,6,7,8]
solve(matD,c)
horizConcat(matD,c)
rank %
rank matD
solve([[2,3,4,5],[3,4,5,6],[4,5,6,7]],[5,6,7])  
subMatrix(matD,1,3,1,4)
rank horizConcat(matD,vector [5,6,7,9])  
solve(matD,[5,6,7,9])            
hilbert3 :: Matrix DoubleFloat -- continuing the previous session
% * inverse %                             
matrix [[1/(i+j) for i in 1..11] for j in 1..11]::Matrix DoubleFloat;
badUnit := % * inverse %;
diagEls := set [%(i,i) for i in 1..11];                                
min diagEls
max diagEls
offDiags := empty()$Set DoubleFloat
for i in 1..11 repeat _
            for j in 1..11 | i ~= j repeat _
               offDiags := union(offDiags,badUnit(i,j))
min offDiags
max offDiags
hilbert11 := matrix [[1/(i+j) for i in 1..11] for j in 1..11];
% * inverse %
detHilbert3 := determinant hilbert3
detHilbert11 := determinant hilbert11
% :: DoubleFloat                                              
determinant(hilbert11::Matrix DoubleFloat)                    
test3 := hilbert3 :: Matrix Polynomial Fraction Integer;
test3(1,1) := (1 + eps)/2;
determinant test3        
(% - detHilbert3)/detHilbert3
for i in 1..3 repeat for j in 1..3 repeat _
            test3(i,j) := hilbert3(i,j) + (2*randnum(2) - 1)*eps
test3
(determinant test3 - detHilbert3)/detHilbert3
error3 := matrix [[eps[i,j] for i in 1..3] for j in 1..3]
test3 := hilbert3 + t*error3;                                          
detErr := (determinant test3 - detHilbert3)/detHilbert3;
detErrReduced := coefficient(%,'t,1)
coefficient(detErr,'t,0)
epses := variables detErrReduced
coefs := coefficients detErrReduced
index := first sort((i,j)+->abs coefs.i > abs coefs.j, expand(1..9)) 
epses.5
coefs.5
sort((i,j)+->abs i > abs j, coefs).2
axiom-test 20170501-3 / usr / share / axiom-20170501 / input / tutchap67.input
