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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
|
\documentclass{article}
\usepackage{axiom}
\begin{document}
\title{\$SPAD/src/input tutChap67.input}
\author{The Axiom Team}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{License}
<<license>>=
--Copyright The Numerical Algorithms Group Limited 1996.
@
<<*>>=
<<license>>
(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
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}
|
