Archivi Categorie: maxima

Maxima – 139 – Matrici e algebra lineare – 6

AphexTwin⁩

Continuo da qui, copio dal Reference Manual, PDF scaricabile da qui, sono a p.415.

matrixexp (M)
matrixexp (M, n)
matrixexp (M, V)
Calculates the matrix exponential e^(M*V). Instead of the vector V a number n can be specified as the second argument. If this argument is omitted matrixexp replaces it by 1.

The matrix exponential of a matrix M can be expressed as a power series: e^M=sum(M^k/k!,0,inf).

matrixmap (f, M)
Returns a matrix with element i,j equal to f(M[i,j]).

See also map, fullmap, fullmapl, and apply.

matrixp (expr)
Returns true if expr is a matrix, otherwise false.

matrix_element_add
Default value: +.

matrix_element_add is the operation invoked in place of addition in a matrix multiplication. matrix_element_add can be assigned any n-ary operator (that is, a function which handles any number of arguments). The assigned value may be the name of an operator enclosed in quote marks, the name of a function, or a lambda expression.

See also matrix_element_mult and matrix_element_transpose.

(%i1) matrix_element_add: "*"$

(%i2) matrix_element_mult: "^"$

(%i3) aa: matrix ([a, b, c], [d, e, f]);
                                  [ a  b  c ]
(%o3)                             [         ]
                                  [ d  e  f ]
(%i4) bb: matrix ([u, v, w], [x, y, z]);
                                  [ u  v  w ]
(%o4)                             [         ]
                                  [ x  y  z ]
(%i5) aa . transpose (bb);
                            [  u  v  w   x  y  z ]
                            [ a  b  c   a  b  c  ]
(%o5)                       [                    ]
                            [  u  v  w   x  y  z ]
                            [ d  e  f   d  e  f  ]

matrix_element_mult
Default value: *.

matrix_element_mult is the operation invoked in place of multiplication in a matrix multiplication. matrix_element_mult can be assigned any binary operator. The assigned value may be the name of an operator enclosed in quote marks, the name of a function, or a lambda expression.

The dot operator . is a useful choice in some contexts.

See also matrix_element_add and matrix_element_transpose.

(%i1) matrix_element_add: lambda ([[x]], sqrt (apply ("+", x)))$

(%i2) matrix_element_mult: lambda ([x, y], (x - y)^2)$

(%i3) [a, b, c] . [x, y, z];
                                 2          2          2
(%o3)                sqrt((c - z)  + (b - y)  + (a - x) )
(%i4) aa: matrix ([a, b, c], [d, e, f]);
                                  [ a  b  c ]
(%o4)                             [         ]
                                  [ d  e  f ]
(%i5) bb: matrix ([u, v, w], [x, y, z]);
                                  [ u  v  w ]
(%o5)                             [         ]
                                  [ x  y  z ]
(%i6) aa . transpose (bb);
(%o6)
 [             2          2          2               2          2          2  ]
 [ sqrt((c - w)  + (b - v)  + (a - u) )  sqrt((c - z)  + (b - y)  + (a - x) ) ]
 [                                                                            ]
 [             2          2          2               2          2          2  ]
 [ sqrt((f - w)  + (e - v)  + (d - u) )  sqrt((f - z)  + (e - y)  + (d - x) ) ]

matrix_element_transpose
Default value: false.

matrix_element_transpose is the operation applied to each element of a matrix when it is transposed. matrix_element_mult can be assigned any unary operator. The assigned value may be the name of an operator enclosed in quote marks, the name of a function, or a lambda expression.

When matrix_element_transpose equals transpose, the transpose function is applied to every element. When matrix_element_transpose equals nonscalars, the transpose function is applied to every nonscalar element. If some element is an atom, the nonscalars option applies transpose only if the atom is declared nonscalar, while the transpose option always applies transpose.

The default value, false, means no operation is applied.

See also matrix_element_add and matrix_element_mult.

(%i1) declare (a, nonscalar)$

(%i2) transpose ([a, b]);
                                     [ a ]
(%o2)                                [   ]
                                     [ b ]
(%i3) matrix_element_transpose: nonscalars$

(%i4) transpose ([a, b]);
                               [ transpose(a) ]
(%o4)                          [              ]
                               [      b       ]
(%i5) matrix_element_transpose: transpose$

(%i6) transpose ([a, b]);
                               [ transpose(a) ]
(%o6)                          [              ]
                               [ transpose(b) ]
(%i7) matrix_element_transpose: lambda ([x], realpart(x) - %i*imagpart(x))$

(%i8) m: matrix ([1 + 5*%i, 3 - 2*%i], [7*%i, 11]);
                            [ 5 %i + 1  3 - 2 %i ]
(%o8)                       [                    ]
                            [   7 %i       11    ]
(%i9) transpose (m);
                             [ 1 - 5 %i  - 7 %i ]
(%o9)                        [                  ]
                             [ 2 %i + 3    11   ]

mattrace (M)
Returns the trace (that is, the sum of the elements on the main diagonal) of the square matrix M.

mattrace is called by ncharpoly, an alternative to Maxima’s charpoly.

load ("nchrpl") loads this function.

minor (M, i, j)
Returns the i, j minor of the matrix M. That is, M with row i and column j removed.

ncharpoly (M, x)
Returns the characteristic polynomial of the matrix M with respect to x. This is an alternative to Maxima’s charpoly.

ncharpoly works by computing traces of powers of the given matrix, which are known to be equal to sums of powers of the roots of the characteristic polynomial. From these quantities the symmetric functions of the roots can be calculated, which are nothing more than the coefficients of the characteristic polynomial. charpoly works by forming the determinant of x * ident [n] - a. Thus ncharpoly wins, for example, in the case of large dense matrices filled with integers, since it avoids polynomial arithmetic altogether.

load ("nchrpl") loads this file.

newdet (M)
Computes the determinant of the matrix M by the Johnson-Gentleman tree minor algorithm. newdet returns the result in CRE form.

permanent (M)
Computes the permanent of the matrix M by the Johnson-Gentleman tree minor algorithm. A permanent is like a determinant but with no sign changes. permanent returns the result in CRE form. See also newdet.

rank (M)
Computes the rank of the matrix M. That is, the order of the largest non-singular subdeterminant of M.

rank may return the wrong answer if it cannot determine that a matrix element that is equivalent to zero is indeed so.

ratmx
Default value: false.

When ratmx is false, determinant and matrix addition, subtraction, and multiplication are performed in the representation of the matrix elements and cause the result of matrix inversion to be left in general representation.

When ratmx is true, the 4 operations mentioned above are performed in CRE form and the result of matrix inverse is in CRE form. Note that this may cause the elements to be expanded (depending on the setting of ratfac) which might not always be desired.

Maxima – 138 – Matrici e algebra lineare – 5

sl91

Continuo da qui, copio dal Reference Manual, PDF scaricabile da qui, sono a p.411.

gramschmidt (x)
gramschmidt (x, F)
Carries out the Gram-Schmidt orthogonalization algorithm on x, which is either a matrix or a list of lists. x is not modified by gramschmidt. The inner product employed by gramschmidt is F, if present, otherwise the inner product is the function innerproduct.

If x is a matrix, the algorithm is applied to the rows of x. If x is a list of lists, the algorithm is applied to the sublists, which must have equal numbers of elements. In either case, the return value is a list of lists, the sublists of which are orthogonal and span the same space as x. If the dimension of the span of x is less than the number of rows or sublists, some sublists of the return value are zero.

factor is called at each stage of the algorithm to simplify intermediate results. As a consequence, the return value may contain factored integers.

load(eigen) loads this function.

Gram-Schmidt algorithm using default inner product function.

(%i1) load ("eigen")$

(%i2) x: matrix ([1, 2, 3], [9, 18, 30], [12, 48, 60]);
                                [ 1   2   3  ]
                                [            ]
(%o2)                           [ 9   18  30 ]
                                [            ]
                                [ 12  48  60 ]
(%i3) y: gramschmidt (x);
                              2      2            4     3
                             3      3   3 5      2  3  2  3
(%o3)         [[1, 2, 3], [- ---, - --, ---], [- ----, ----, 0]]
                             2 7    7   2 7       5     5
(%i4) map (innerproduct, [y[1], y[2], y[3]], [y[2], y[3], y[1]]);
(%o4)                              [0, 0, 0]

Gram-Schmidt algorithm using a specified inner product function.

(%i5) ip (f, g) := integrate (f * g, u, a, b);
(%o5)                 ip(f, g) := integrate(f g, u, a, b)
(%i6) y : gramschmidt([1, sin(u), cos(u)], ip), a= -%pi/2, b=%pi/2;
                                      %pi cos(u) - 2
(%o6)                     [1, sin(u), --------------]
                                           %pi
(%i7) map (ip, [y[1], y[2], y[3]], [y[2], y[3], y[1]]), a= -%pi/2, b=%pi/2;
(%o7)                              [0, 0, 0]

ident (n)
Returns an n by n identity matrix.

innerproduct (x, y)
inprod (x, y)
Returns the inner product (also called the scalar product or dot product) of x and y, which are lists of equal length, or both 1-column or 1-row matrices of equal length.

The return value is conjugate (x) . y, where . is the noncommutative multiplication operator.

load ("eigen") loads this function.

inprod is a synonym for innerproduct.

invert_by_adjoint (M)
Returns the inverse of the matrix M. The inverse is computed by the adjoint method. invert_by_adjoint honors the ratmx and detout flags, the same as invert.

invert (M)
Returns the inverse of the matrix M. The inverse is computed via the LU decomposition.

When ratmx is true, elements of M are converted to canonical rational expressions (CRE), and the elements of the return value are also CRE.

When ratmx is false, elements of M are not converted to a common representation. In particular, float and bigfloat elements are not converted to rationals.

When detout is true, the determinant is factored out of the inverse. The global flags doallmxops and doscmxops must be false to prevent the determinant from being absorbed into the inverse. xthru can multiply the determinant into the inverse.

invert does not apply any simplifications to the elements of the inverse apart from the default arithmetic simplifications. ratsimp and expand can apply additional simplifications. In particular, when M has polynomial elements, expand(invert(M)) might be preferable.

invert(M) is equivalent to M^^-1.

list_matrix_entries (M)
Returns a list containing the elements of the matrix M.

(%i8) list_matrix_entries(matrix([a,b],[c,d]));
(%o8)                            [a, b, c, d]

lmxchar
Default value: [.
lmxchar is the character displayed as the left delimiter of a matrix. See also rmxchar.

(%i9) lmxchar: "|"$

(%i10) matrix ([a, b, c], [d, e, f], [g, h, i]);
                                  | a  b  c ]
                                  |         ]
(%o10)                            | d  e  f ]
                                  |         ]
                                  | g  h  i ]

matrix (row_1, ..., row_n)
Returns a rectangular matrix which has the rows row_1, ..., row_n. Each row is a list of expressions. All rows must be the same length.

The operations + (addition), - (subtraction), * (multiplication), and / (division), are carried out element by element when the operands are two matrices, a scalar and a matrix, or a matrix and a scalar. The operation ^ (exponentiation, equivalently **) is carried out element by element if the operands are a scalar and a matrix or a matrix and a scalar, but not if the operands are two matrices. All operations are normally carried out in full, including . (noncommutative multiplication).

Matrix multiplication is represented by the noncommutative multiplication operator ... The corresponding noncommutative exponentiation operator is ^^. For a matrix A, A.A = A^^2 and A^^-1 is the inverse of A, if it exists. A^^-1 is equivalent to invert(A).

There are switches for controlling simplification of expressions involving dot and matrix-list operations. These are doallmxops, domxexpt, domxmxops, doscmxops, and doscmxplus.

There are additional options which are related to matrices. These are: lmxchar, rmxchar, ratmx, listarith, detout, scalarmatrix and sparse.

There are a number of functions which take matrices as arguments or yield matrices as return values. See eigenvalues, eigenvectors, determinant, charpoly, genmatrix, addcol, addrow, copymatrix, transpose, echelon, and rank.

Construction of matrices from lists.

(%i1) x: matrix ([17, 3], [-8, 11]);
                                  [ 17   3  ]
(%o1)                             [         ]
                                  [ - 8  11 ]
(%i2) y: matrix ([%pi, %e], [a, b]);
                                  [ %pi  %e ]
(%o2)                             [         ]
                                  [  a   b  ]

Addition, element by element.

(%i3) x + y;
                             [ %pi + 17  %e + 3 ]
(%o3)                        [                  ]
                             [  a - 8    b + 11 ]

Subtraction, element by element.

(%i4) x - y;
                             [ 17 - %pi   3 - %e ]
(%o4)                        [                   ]
                             [ (- a) - 8  11 - b ]

Multiplication, element by element.

(%i5) x * y;
                               [ 17 %pi  3 %e ]
(%o5)                          [              ]
                               [ - 8 a   11 b ]

Division, element by element.

(%i6) x / y;
                               [ 17       - 1 ]
                               [ ---  3 %e    ]
                               [ %pi          ]
(%o6)                          [              ]
                               [   8    11    ]
                               [ - -    --    ]
                               [   a    b     ]

Matrix to a scalar exponent, element by element.

(%i7) x ^ 3;
                                [ 4913    27  ]
(%o7)                           [             ]
                                [ - 512  1331 ]

Scalar base to a matrix exponent, element by element.

(%i8) exp(y);
                                [   %pi    %e ]
                                [ %e     %e   ]
(%o8)                           [             ]
                                [    a     b  ]
                                [  %e    %e   ]

Matrix base to a matrix exponent. This is not carried out element by element. See also matrixexp.

(%i9) x ^ y;
                                       [ %pi  %e ]
                                       [         ]
                                       [  a   b  ]
                            [ 17   3  ]
(%o9)                       [         ]
                            [ - 8  11 ]

Noncommutative matrix multiplication.

(%i10) x . y;
                         [ 3 a + 17 %pi  3 b + 17 %e ]
(%o10)                   [                           ]
                         [ 11 a - 8 %pi  11 b - 8 %e ]
(%i11) y . x;
                       [ 17 %pi - 8 %e  3 %pi + 11 %e ]
(%o11)                 [                              ]
                       [  17 a - 8 b     11 b + 3 a   ]

Noncommutative matrix exponentiation. A scalar base b to a matrix power M is carried out element by element and so b^^m is the same as b^m.

(%i12) x ^^ 3;
                               [  3833   1719 ]
(%o12)                         [              ]
                               [ - 4584  395  ]
(%i13) %e ^^ y;
                                   [ %pi  %e ]
                                  <[         ]>
                                   [  a   b  ]
(%o13)                          %e

Nota: sul manuale riporta:

(%i13) %e ^^ y;
                         [   %pi    %e ]
                         [ %e     %e   ]
(%o13)                   [             ]
                         [    a     b  ]
                         [  %e    %e   ]

che è ovviamente equivalente; forse qualche flag…

A matrix raised to a -1 exponent with noncommutative exponentiation is the matrix inverse, if it exists.

(%i14) x ^^ -1;
                                [ 11      3  ]
                                [ ---  - --- ]
                                [ 211    211 ]
(%o14)                          [            ]
                                [  8    17   ]
                                [ ---   ---  ]
                                [ 211   211  ]
(%i15) x . (x ^^ -1);
                                   [ 1  0 ]
(%o15)                             [      ]
                                   [ 0  1 ]

Maxima – 137 – Matrici e algebra lineare – 4

Dqpv

Continuo da qui, copio dal Reference Manual, PDF scaricabile da qui, sono a p.406.

eigenvalues (M)
eivals (M)
Returns a list of two lists containing the eigenvalues of the matrix M. The first sublist of the return value is the list of eigenvalues of the matrix, and the second sublist is the list of the multiplicities of the eigenvalues in the corresponding order.

eivals is a synonym for eigenvalues.

eigenvalues calls the function solve to find the roots of the characteristic polynomial of the matrix. Sometimes solve may not be able to find the roots of the polynomial; in that case some other functions in this package (except innerproduct, unitvector, columnvector and gramschmidt) will not work. Sometimes solve may find only a subset of the roots of the polynomial. This may happen when the factoring of the polynomial contains polynomials of degree 5 or more. In such cases a warning message is displayed and the only the roots found and their corresponding multiplicities are returned.

In some cases the eigenvalues found by solve may be complicated expressions. (This may happen when solve returns a not-so-obviously real expression for an eigenvalue which is known to be real.) It may be possible to simplify the eigenvalues using some other functions.

The package eigen.mac is loaded automatically when eigenvalues or eigenvectors is referenced. If eigen.mac is not already loaded, load ("eigen") loads it. After loading, all functions and variables in the package are available.

eigenvectors (M)
eivects (M)
Computes eigenvectors of the matrix M. The return value is a list of two elements. The first is a list of the eigenvalues of M and a list of the multiplicities of the eigenvalues. The second is a list of lists of eigenvectors. There is one list of eigenvectors for each eigenvalue. There may be one or more eigenvectors in each list.

eivects is a synonym for eigenvectors.

The package eigen.mac is loaded automatically when eigenvalues or eigenvectors is referenced. If eigen.mac is not already loaded, load ("eigen") loads it. After loading, all functions and variables in the package are available.

Note that eigenvectors internally calls eigenvalues to obtain eigenvalues. So, when eigenvalues returns a subset of all the eigenvalues, the eigenvectors returns the corresponding subset of the all the eigenvectors, with the same warning displayed as eigenvalues.

The flags that affect this function are:

  • nondiagonalizable is set to true or false depending on whether the matrix is nondiagonalizable or diagonalizable after eigenvectors returns.
  • hermitianmatrix when true, causes the degenerate eigenvectors of the Hermitian matrix to be orthogonalized using the Gram-Schmidt algorithm.
  • knowneigvals when true causes the eigen package to assume the eigenvalues of the matrix are known to the user and stored under the global name listeigvals.
  • listeigvals should be set to a list similar to the output eigenvalues.

The function algsys is used here to solve for the eigenvectors. Sometimes if the eigenvalues are messy, algsys may not be able to find a solution. In some cases, it may be possible to simplify the eigenvalues by first finding them using eigenvalues command and then using other functions to reduce them to something simpler. Following simplification, eigenvectors can be called again with the knowneigvals flag set to true.

See also eigenvalues.

A matrix which has just one eigenvector per eigenvalue.

(%i1) M1 : matrix ([11, -1], [1, 7]);
                                  [ 11  - 1 ]
(%o1)                             [         ]
                                  [ 1    7  ]
(%i2) [vals, vecs] : eigenvectors (M1);
(%o2) [[[9 - sqrt(3), sqrt(3) + 9], [1, 1]],
                                      [[[1, sqrt(3) + 2]], [[1, 2 - sqrt(3)]]]]
(%i3) for i thru length (vals[1]) do disp (val[i] = vals[1][i],
          mult[i] = vals[2][i], vec[i] = vecs[i]);
                              val  = 9 - sqrt(3)
                                 1

                                   mult  = 1
                                       1

                           vec  = [[1, sqrt(3) + 2]]
                              1

                              val  = sqrt(3) + 9
                                 2

                                   mult  = 1
                                       2

                           vec  = [[1, 2 - sqrt(3)]]
                              2

(%o3)                                done

A matrix which has two eigenvectors for one eigenvalue (namely 2).

(%i4) M1 : matrix ([0, 1, 0, 0], [0, 0, 0, 0], [0, 0, 2, 0], [0, 0, 0, 2]);
                                [ 0  1  0  0 ]
                                [            ]
                                [ 0  0  0  0 ]
(%o4)                           [            ]
                                [ 0  0  2  0 ]
                                [            ]
                                [ 0  0  0  2 ]
(%i5) [vals, vecs] : eigenvectors (M1);
(%o5) [[[0, 2], [2, 2]], [[[1, 0, 0, 0]], [[0, 0, 1, 0], [0, 0, 0, 1]]]]
(%i6) for i thru length (vals[1]) do disp (val[i] = vals[1][i],
          mult[i] = vals[2][i], vec[i] = vecs[i]);
                                   val  = 0
                                      1

                                   mult  = 2
                                       1

                             vec  = [[1, 0, 0, 0]]
                                1

                                   val  = 2
                                      2

                                   mult  = 2
                                       2

                      vec  = [[0, 0, 1, 0], [0, 0, 0, 1]]
                         2

(%o6)                                done

ematrix (m, n, x, i, j)
Returns an m by n matrix, all elements of which are zero except for the [i, j] element which is x.

entermatrix (m, n)
Returns an m by n matrix, reading the elements interactively.

If n is equal to m, Maxima prompts for the type of the matrix (diagonal, symmetric, antisymmetric, or general) and for each element. Each response is terminated by a semicolon ; or dollar sign $.

If n is not equal to m, Maxima prompts for each element.

The elements may be any expressions, which are evaluated. entermatrix evaluates its arguments.

(%i7) n: 3$

(%i8) m: entermatrix (n, n)$

Is the matrix  1. Diagonal  2. Symmetric  3. Antisymmetric  4. General
Answer 1, 2, 3 or 4 :
1$
Row 1 Column 1:
(a+b)^n$
Row 2 Column 2:
(a+b)^(n+1)$
Row 3 Column 3:
(a+b)^(n+2)$

Matrix entered.
(%i9) m;
                       [        3                     ]
                       [ (b + a)      0         0     ]
                       [                              ]
(%o9)                  [                  4           ]
                       [    0      (b + a)      0     ]
                       [                              ]
                       [                            5 ]
                       [    0         0      (b + a)  ]

genmatrix (a, i_2, j_2, i_1, j_1)
genmatrix (a, i_2, j_2, i_1)
genmatrix (a, i_2, j_2)
Returns a matrix generated from a, taking element a[i_1, j_1] as the upper-left element and a[i_2, j_2] as the lower-right element of the matrix. Here a is a declared array (created by array but not by make_array) or an undeclared array, or an array function, or a lambda expression of two arguments. (An array function is created like other functions with := or define, but arguments are enclosed in square brackets instead of parentheses.)

If j_1 is omitted, it is assumed equal to i 1. If both j_1 and i_1 are omitted, both are assumed equal to 1.

If a selected element i,j of the array is undefined, the matrix will contain a symbolic element a[i,j].

(%i10) h [i, j] := 1 / (i + j - 1);
                                           1
(%o10)                        h     := ---------
                               i, j    i + j - 1
(%i11) genmatrix (h, 3, 3);
                                  [    1  1 ]
                                  [ 1  -  - ]
                                  [    2  3 ]
                                  [         ]
                                  [ 1  1  1 ]
(%o11)                            [ -  -  - ]
                                  [ 2  3  4 ]
                                  [         ]
                                  [ 1  1  1 ]
                                  [ -  -  - ]
                                  [ 3  4  5 ]
(%i12) array (a, fixnum, 2, 2);
(%o12)                                 a
(%i13) a [1, 1] : %e;
(%o13)                                %e
(%i14) a [2, 2] : %pi;
(%o14)                                %pi
(%i15) genmatrix (a, 2, 2);
                                  [ %e   0  ]
(%o15)                            [         ]
                                  [ 0   %pi ]
(%i16) genmatrix (lambda ([i, j], j - i), 3, 3);
                                [  0    1   2 ]
                                [             ]
(%o16)                          [ - 1   0   1 ]
                                [             ]
                                [ - 2  - 1  0 ]
(%i17) genmatrix (B, 2, 2);
                               [ B      B     ]
                               [  1, 1   1, 2 ]
(%o17)                         [              ]
                               [ B      B     ]
                               [  2, 1   2, 2 ]

Maxima – 136 – Matrici e algebra lineare – 3

sl92

Continuo da qui, copio dal Reference Manual, PDF scaricabile da qui, sono a p.400.

determinant (M)
Computes the determinant of M by a method similar to Gaussian elimination.

The form of the result depends upon the setting of the switch ratmx.

There is a special routine for computing sparse determinants which is called when the switches ratmx and sparse are both true.

detout
Default value: false.

When detout is true, the determinant of a matrix whose inverse is computed is factored out of the inverse.

For this switch to have an effect doallmxops and doscmxops should be false (see their descriptions). Alternatively this switch can be given to ev which causes the other two to be set correctly.

(%i1) m: matrix ([a, b], [c, d]);
                                   [ a  b ]
(%o1)                              [      ]
                                   [ c  d ]
(%i2) detout: true$

(%i3) doallmxops: false$

(%i4) doscmxops: false$

(%i5) invert (m);
                                 [  d   - b ]
                                 [          ]
                                 [ - c   a  ]
(%o5)                            ------------
                                  a d - b c

diagmatrix (n, x)
Returns a diagonal matrix of size n by n with the diagonal elements all equal to x.

diagmatrix (n, 1) returns an identity matrix (same as ident (n)).

n must evaluate to an integer, otherwise diagmatrix complains with an error message.

x can be any kind of expression, including another matrix. If x is a matrix, it is not copied; all diagonal elements refer to the same instance, x.

doallmxops
Default value: true.

When doallmxops is true, all operations relating to matrices are carried out. When it is false then the setting of the individual dot switches govern which operations are performed.

domxexpt
Default value: true.

When domxexpt is true, a matrix exponential, exp (M) where M is a matrix, is interpreted as a matrix with element [i,j] equal to exp (m[i,j]). Otherwise exp (M) evaluates to exp (ev(M)).

domxexpt affects all expressions of the form base^power where base is an expression assumed scalar or constant, and power is a list or matrix.

(%i1) m: matrix ([1, %i], [a+b, %pi]);
                                [   1    %i  ]
(%o1)                           [            ]
                                [ b + a  %pi ]
(%i2) domxexpt: false$

(%i3) (1 - c)^m;
                                    [   1    %i  ]
                                    [            ]
                                    [ b + a  %pi ]
(%o3)                        (1 - c)
(%i4) domxexpt: true$

(%i5) (1 - c)^m;
                         [                      %i  ]
                         [    1 - c      (1 - c)    ]
(%o5)                    [                          ]
                         [        b + a         %pi ]
                         [ (1 - c)       (1 - c)    ]

domxmxops
Default value: true.

When domxmxops is true, all matrix-matrix or matrix-list operations are carried out (but not scalar-matrix operations); if this switch is false such operations are not carried out.

domxnctimes
Default value: false.

When domxnctimes is true, non-commutative products of matrices are carried out.

dontfactor
Default value: [].

dontfactor may be set to a list of variables with respect to which factoring is not to occur. (The list is initially empty.) Factoring also will not take place with respect to any variables which are less important, according the variable ordering assumed for canonical rational expression (CRE) form, than those on the dontfactor list.

doscmxops
Default value: false.

When doscmxops is true, scalar-matrix operations are carried out.

doscmxplus
Default value: false.

When doscmxplus is true, scalar-matrix operations yield a matrix result. This switch is not subsumed under doallmxops.

dot0nscsimp
Default value: true.

When dot0nscsimp is true, a non-commutative product of zero and a nonscalar term is simplified to a commutative product.

dot0simp
Default value: true.

When dot0simp is true, a non-commutative product of zero and a scalar term is simplified to a commutative product.

dot1simp
Default value: true.

When dot1simp is true, a non-commutative product of one and another term is simplified to a commutative product.

dotassoc
Default value: true.

When dotassoc is true, an expression (A.B).C simplifies to A.(B.C).

dotconstrules
Default value: true.

When dotconstrules is true, a non-commutative product of a constant and another term is simplified to a commutative product. Turning on this flag effectively turns on dot0simp, dot0nscsimp, and dot1simp as well.

dotdistrib
Default value: false.

When dotdistrib is true, an expression A.(B + C) simplifies to A.B + A.C.

dotexptsimp
Default value: true.

When dotexptsimp is true, an expression A.A simplifies to A^^2.

dotident
Default value: 1.

dotident is the value returned by X^^0.

dotscrules
Default value: false.

When dotscrules is true, an expression A.SC or SC.A simplifies to SC*A and A.(SC*B) simplifies to SC*(A.B).

echelon (M)
Returns the echelon form of the matrix M, as produced by Gaussian elimination. The echelon form is computed from M by elementary row operations such that the first non-zero element in each row in the resulting matrix is one and the column elements under the first one in each row are all zero.

triangularize also carries out Gaussian elimination, but it does not normalize the leading non-zero element in each row.

lu_factor and cholesky are other functions which yield triangularized matrices.

(%i1) M: matrix ([3, 7, aa, bb], [-1, 8, 5, 2], [9, 2, 11, 4]);
                              [  3   7  aa  bb ]
                              [                ]
(%o1)                         [ - 1  8  5   2  ]
                              [                ]
                              [  9   2  11  4  ]
(%i2) echelon (M);
                         [ 1  - 8  - 5      - 2     ]
                         [                          ]
                         [         28       11      ]
                         [ 0   1   --       --      ]
(%o2)                    [         37       37      ]
                         [                          ]
                         [              37 bb - 119 ]
                         [ 0   0    1   ----------- ]
                         [              37 aa - 313 ]

Maxima – 135 – Matrici e algebra lineare – 2

NGC1300

Continuo da qui, copio dal Reference Manual, PDF scaricabile da qui, sono a p.400.

addcol (M, list_1, ..., list_n)
Appends the column(s) given by the one or more lists (or matrices) onto the matrix M. See also addrow and append.

addrow (M, list_1, ..., list_n)
Appends the row(s) given by the one or more lists (or matrices) onto the matrix M. See also addcol and append.

adjoint (M)
Returns the adjoint of the matrix M. The adjoint matrix is the transpose of the
matrix of cofactors of M.

augcoefmatrix ([eqn_1, ..., eqn_m], [x_1, ..., x_n])
Returns the augmented coefficient matrix for the variables x_1, ..., x_n of the system of linear equations eqn_1, ..., eqn_m. This is the coefficient matrix with a column adjoined for the constant terms in each equation (i.e., those terms not dependent upon x_1, ..., x_n).

(%i1) m: [2*x - (a - 1)*y = 5*b, c + b*y + a*x = 0]$

(%i2) augcoefmatrix (m, [x, y]);
                              [ 2  1 - a  - 5 b ]
(%o2)                         [                 ]
                              [ a    b      c   ]

cauchy_matrix ([x_1, x_2, ..., x_m], [y_1, y_2, ..., y_n])
cauchy_matrix ([x_1, x_2, ..., x_n])
Returns a n by m Cauchy matrix with the elements a[i,j] = 1/(x_i+y_i). The second argument of cauchy_matrix is optional. For this case the elements of the Cauchy
matrix are a[i,j] = 1/(x_i+x_j).

Remark: In the literature the Cauchy matrix can be found defined in two forms. A second definition is a[i,j] = 1/(x_i-y_i).

(%i3) cauchy_matrix([x1,x2],[y1,y2]);
                             [    1        1    ]
                             [ -------  ------- ]
                             [ y1 + x1  y2 + x1 ]
(%o3)                        [                  ]
                             [    1        1    ]
                             [ -------  ------- ]
                             [ y1 + x2  y2 + x2 ]
(%i4) cauchy_matrix([x1,x2]);
                             [   1         1    ]
                             [  ----    ------- ]
                             [  2 x1    x2 + x1 ]
(%o4)                        [                  ]
                             [    1       1     ]
                             [ -------   ----   ]
                             [ x2 + x1   2 x2   ]

charpoly (M, x)
Returns the characteristic polynomial for the matrix M with respect to variable x. That is, determinant (M - diagmatrix (length (M), x)).

(%i1) a: matrix ([3, 1], [2, 4]);
                                   [ 3  1 ]
(%o1)                              [      ]
                                   [ 2  4 ]
(%i2) expand (charpoly (a, lambda));
                                  2
(%o2)                       lambda  - 7 lambda + 10
(%i3) (programmode: true, solve (%));
(%o3)                      [lambda = 5, lambda = 2]
(%i4) matrix ([x1], [x2]);
                                    [ x1 ]
(%o4)                               [    ]
                                    [ x2 ]
(%i5) ev (a . % - lambda*%, %th(2)[1]);
                                 [ x2 - 2 x1 ]
(%o5)                            [           ]
                                 [ 2 x1 - x2 ]
(%i6) %[1, 1] = 0;
(%o6)                            x2 - 2 x1 = 0
(%i7) x2^2 + x1^2 = 1;
                                   2     2
(%o7)                            x2  + x1  = 1
(%i8) solve ([%th(2), %], [x1, x2]);
                   1               2               1             2
(%o8)  [[x1 = - -------, x2 = - -------], [x1 = -------, x2 = -------]]
                sqrt(5)         sqrt(5)         sqrt(5)       sqrt(5)

Nota: il codice riportato nel PDF che sto seguendo da errore; la versione corretta è qui.

coefmatrix ([eqn_1, ..., eqn_m], [x_1, ..., x_n])
Returns the coefficient matrix for the variables x_1, ..., x_n of the system of linear equations eqn_1, ..., eqn_m.

(%i1) coefmatrix([2*x-(a-1)*y+5*b = 0, b*y+a*x = 3], [x,y]);
                                 [ 2  1 - a ]
(%o1)                            [          ]
                                 [ a    b   ]

col (M, i)
Returns the i’th column of the matrix M. The return value is a matrix.

columnvector (L)
covect (L)
Returns a matrix of one column and length (L) rows, containing the elements of the list L.

covect is a synonym for columnvector.

load ("eigen") loads this function.

This is useful if you want to use parts of the outputs of the functions in this package in matrix calculations.

(%i1) load ("eigen")$

(%i2) columnvector ([aa, bb, cc, dd]);
                                    [ aa ]
                                    [    ]
                                    [ bb ]
(%o2)                               [    ]
                                    [ cc ]
                                    [    ]
                                    [ dd ]

copymatrix (M)
Returns a copy of the matrix M. This is the only way to make a copy aside from copying M element by element.

Note that an assignment of one matrix to another, as in m2: m1, does not copy m1. An assignment m2 [i,j]: x or setelmx(x, i, j, m2) also modifies m1 [i,j]. Creating a copy with copymatrix and then using assignment creates a separate, modified copy.

Maxima – 134 – Matrici e algebra lineare – 1

v8HE

Continuo da qui, copio dal Reference Manual, PDF scaricabile da qui, sono a p.399.

Dot
The operator . represents noncommutative multiplication and scalar product. When the operands are 1-column or 1-row matrices a and b, the expression a.b is equivalent to sum (a[i]*b[i], i, 1, length(a)). If a and b are not complex, this is the scalar product, also called the inner product or dot product, of a and b. The scalar product is defined as conjugate(a).b when a and b are complex; innerproduct in the eigen package provides the complex scalar product.

When the operands are more general matrices, the product is the matrix product a and b. The number of rows of b must equal the number of columns of a, and the result has number of rows equal to the number of rows of a and number of columns equal to the number of columns of b.

To distinguish . as an arithmetic operator from the decimal point in a floating point number, it may be necessary to leave spaces on either side. For example, 5.e3 is 5000.0 but 5 . e3 is 5 times e3.

There are several flags which govern the simplification of expressions involving ., namely dot0nscsimp, dot0simp, dot1simp, dotassoc, dotconstrules, dotdistrib, dotexptsimp, dotident, and dotscrules.

Vettori
vect is a package of functions for vector analysis. load ("vect") loads this package, and demo ("vect") displays a demonstration.

The vector analysis package can combine and simplify symbolic expressions including dot products and cross products, together with the gradient, divergence, curl, and Laplacian operators. The distribution of these operators over sums or products is governed by several flags, as are various other expansions, including expansion into components in any specific orthogonal coordinate systems. There are also functions for deriving the scalar or vector potential of a field.

The vect package contains these functions: vectorsimp, scalefactors, express, potential, and vectorpotential.

By default the vect package does not declare the dot operator to be a commutative
operator. To get a commutative dot operator ., the command declare(".", commutative) must be executed.

eigen
The package eigen contains several functions devoted to the symbolic computation of eigenvalues and eigenvectors. Maxima loads the package automatically if one of the functions eigenvalues or eigenvectors is invoked. The package may be loaded explicitly as load ("eigen").

demo ("eigen") displays a demonstration of the capabilities of this package. batch ("eigen") executes the same demonstration, but without the user prompt between successive computations.

The functions in the eigen package are: innerproduct, unitvector, columnvector, gramschmidt, eigenvalues, eigenvectors, uniteigenvectors, and similaritytransform.

Maxima – 133 – Funzioni numeriche – 5

DrRr

Continuo da qui, copio dal Reference Manual, PDF scaricabile da qui, sono a p.395.

ploteq (exp, options...)
Plots equipotential curves for exp, which should be an expression depending on two variables. The curves are obtained by integrating the differential equation that define the orthogonal trajectories to the solutions of the autonomous system obtained from the gradient of the expression given. The plot can also show the integral curves for that gradient system (option fieldlines).

This program also requires Xmaxima, even if its run from a Maxima session in a console, since the plot will be created by the Tk scripts in Xmaxima. By default, the plot region will be empty until the user clicks in a point (or gives its coordinate with in the set-up menu or via the trajectory_at option).

Most options accepted by plotdf can also be used for ploteq and the plot interface is the same that was described in plotdf.

(%i1) V: 900/((x+1)^2+y^2)^(1/2)-900/((x-1)^2+y^2)^(1/2)$

(%i2) ploteq(V,[x,-2,2],[y,-2,2],[fieldlines,"blue"])$

133-0

Clicking on a point will plot the equipotential curve that passes by that point (in red) and the orthogonal trajectory (in blue).

rk (ODE, var, initial, domain)
rk ([ODE1, . . . , ODEm], [v1, ..., vm], [init1, ..., initm], domain)
The first form solves numerically one first-order ordinary differential equation, and the second form solves a system of m of those equations, using the 4th order Runge-Kutta method. var represents the dependent variable. ODE must be an expression that depends only on the independent and dependent variables and defines the derivative of the dependent variable with respect to the independent variable.

The independent variable is specified with domain, which must be a list of four elements as, for instance [t, 0, 10, 0.1]the first element of the list identifies the independent variable, the second and third elements are the initial and final values for that variable, and the last element sets the increments that should be used within that interval.

If m equations are going to be solved, there should be m dependent variables v1, v2, ..., vm. The initial values for those variables will be init1, init2, ..., initm. There will still be just one independent variable defined by domain, as in the previous case.

ODE1, ..., ODEm are the expressions that define the derivatives of each dependent variable in terms of the independent variable. The only variables that may appear in those expressions are the independent variable and any of the dependent variables. It is important to give the derivatives ODE1, ..., ODEm in the list in exactly the same order used for the dependent variables; for instance, the third element in the list will be interpreted as the derivative of the third dependent variable.

The program will try to integrate the equations from the initial value of the independent variable until its last value, using constant increments. If at some step one of the dependent variables takes an absolute value too large, the integration will be interrupted at that point. The result will be a list with as many elements as the number of iterations made. Each element in the results list is itself another list with m+1 elements: the value of the independent variable, followed by the values of the dependent variables corresponding to that point.

To solve numerically the differential equation

133-1

With initial value x(t=0) = 1, in the interval of t from 0 to 8 and with increments of 0.1 for t, use:

(%i3) results: rk(t-x^2,x,1,[t,0,8,0.1])$

(%i4) plot2d ([discrete, results])$

133-2

the results will be saved in the list results and the plot will show the solution
obtained, with t on the horizontal axis and x on the vertical axis.

To solve numerically the system:

133-3

for t between 0 and 4, and with values of -1.25 and 0.75 for x and y at t=0:

(%i5) sol: rk([4-x^2-4*y^2,y^2-x^2+1],[x,y],[-1.25,0.75],[t,0,4,0.02])$

(%i6) plot2d ([discrete,makelist([p[1],p[3]],p,sol)], [xlabel,"t"],[ylabel,"y"])$

133-4

Maxima – 132 – Funzioni numeriche – 4

ftc2

Continuo da qui, copio dal Reference Manual, PDF scaricabile da qui, sono a p.390.

Funzioni per soluzioni numeriche di equazioni differenziali

The Ordinary Differential Equations (ODE) solved by the functions in this section should have the form,

132-0

which is a first-order ODE. Higher order differential equations of order n must be written as a system of n first-order equations of that kind. For instance, a second-order ODE should be written as a system of two equations

132-1

The first argument in the functions will be a list with the expressions on the right-side of the ODE’s. The variables whose derivatives are represented by those expressions should be given in a second list. In the case above those variables are x and y. The independent variable, t in the examples above, might be given in a separated option. If the expressions given do not depend on that independent variable, the system is called autonomous.

plotdf (dydx, options...)
plotdf (dvdu, [u, v], options...)
plotdf ([dxdt, cdydt], options...)
plotdf ([dudt,cdvdt], [u,cv], options...)
The function plotdf creates a two-dimensional plot of the direction field (also called slope field) for a first-order Ordinary Differential Equation (ODE) or a system of two autonomous first-order ODE’s.

plotdf requires Xmaxima, even if its run from a Maxima session in a console, since the plot will be created by the Tk scripts in Xmaxima. If Xmaxima is not installed plotdf will not work.

dydx, dxdt and dydt are expressions that depend on x and y. dvdu, dudt and dvdt are expressions that depend on u and v. In addition to those two variables, the expressions can also depend on a set of parameters, with numerical values given with the parameters option (the option syntax is given below), or with a range of allowed values specified by a sliders option.

Several other options can be given within the command, or selected in the menu. Integral curves can be obtained by clicking on the plot, or with the option trajectory_at.

The direction of the integration can be controlled with the direction option, which can have values of forward, backward or both. The number of integration steps is given by nsteps; at each integration step the time increment will be adjusted automatically to produce displacements much smaller than the size of the plot window.

The numerical method used is 4th order Runge-Kutta with variable time steps.

Plot window menu:
The menu bar of the plot window has the following seven icons:

  • An X. Can be used to close the plot window.
  • A wrench and a screwdriver. Opens the configuration menu with several fields that show the ODE(s) in use and various other settings. If a pair of coordinates are entered in the field Trajectory at and the enter key is pressed, a new integral curve will be shown, in addition to the ones already shown.
  • Two arrows following a circle. Replots the direction field with the new settings defined in the configuration menu and replots only the last integral curve that was previously plotted.
  • Hard disk drive with an arrow. Used to save a copy of the plot, in Postscript format, in the file specified in a field of the box that appears when that icon is clicked.
  • Magnifying glass with a plus sign. Zooms in the plot.
  • Magnifying glass with a minus sign. Zooms out the plot. The plot can be displaced by holding down the right mouse button while the mouse is moved.
  • Icon of a plot. Opens another window with a plot of the two variables in terms of time, for the last integral curve that was plotted.

Plot options:
Options can also be given within the plotdf itself, each one being a list of two or more elements. The first element in each option is the name of the option, and the remainder is the value or values assigned to the option.

The options which are recognized by plotdf are the following:

  • nsteps defines the number of steps that will be used for the independent variable, to compute an integral curve. The default value is 100.
  • direction defines the direction of the independent variable that will be followed to compute an integral curve. Possible values are forward, to make the independent variable increase nsteps times, with increments tstep, backward, to make the independent variable decrease, or both that will lead to an integral curve that extends nsteps forward, and nsteps backward. The keywords right and left can be used as synonyms for forward and backward. The default value is both.
  • tinitial defines the initial value of variable t used to compute integral curves.
    Since the differential equations are autonomous, that setting will only appear in
    the plot of the curves as functions of t. The default value is 0.
  • versus_t is used to create a second plot window, with a plot of an integral curve, as two functions x, y, of the independent variable t. If versus_t is given any value different from 0, the second plot window will be displayed. The second plot window includes another menu, similar to the menu of the main plot window. The default value is 0.
  • trajectory at defines the coordinates xinitial and yinitial for the starting point of an integral curve. The option is empty by default.
  • parameters defines a list of parameters, and their numerical values, used in the definition of the differential equations. The name and values of the parameters must be given in a string with a comma-separated sequence of pairs name=value.
  • sliders defines a list of parameters that will be changed interactively using slider buttons, and the range of variation of those parameters. The names and ranges of the parameters must be given in a string with a comma-separated sequence of elements name=min:max.
  • xfun defines a string with semi-colon-separated sequence of functions of x to be displayed, on top of the direction field. Those functions will be parsed by Tcl and not by Maxima.
  • x should be followed by two numbers, which will set up the minimum and maximum values shown on the horizontal axis. If the variable on the horizontal axis is not x, then this option should have the name of the variable on the horizontal axis. The default horizontal range is from -10 to 10.
  • y should be followed by two numbers, which will set up the minimum and maximum values shown on the vertical axis. If the variable on the vertical axis is not y, then this option should have the name of the variable on the vertical axis. The default vertical range is from -10 to 10.
  • xaxislabel will be used to identify the horizontal axis. Its default value is the name of the first state variable.
  • yaxislabel will be used to identify the vertical axis. Its default value is the name of the second state variable.
  • number of arrows should be set to a square number and defines the approximate density of the arrows being drawn. The default value is 225.

Esempi
To show the direction field of the differential equation y' = exp(−x) + y and the solution that goes through (2, −0.1):

(%i1) plotdf(exp(-x)+y,[trajectory_at,2,-0.1])$
/bin/sh: 1: /usr/bin/xmaxima: not found

OOPS! devo installare Xmaxima, con sudo apt install xmaxima.

(%i1) plotdf(exp(-x)+y,[trajectory_at,2,-0.1])$

132-2

To obtain the direction field for the equation dif f (y, x) = x − y^2 and the solution with initial condition y(−1) = 3, we can use the command:

(%i2) plotdf(x-y^2,[xfun,"sqrt(x);-sqrt(x)"],
             [trajectory_at,-1,3], [direction,forward],
             [y,-5,5], [x,-4,16])$

132-3

The following example shows the direction field of a harmonic oscillator, defined by the two equations dz/dt = v and dv/dt = −k ∗ z/m, and the integral curve through (z, v) = (6, 0), with a slider that will allow you to change the value of m interactively (k is fixed at 2):

(%i3) plotdf([v,-k*z/m], [z,v], [parameters,"m=2,k=2"],
             [sliders,"m=1:5"], [trajectory_at,6,0])$

132-4

To plot the direction field of the Duffing equation, m∗x''+c∗x'+k∗x+b∗x^3 = 0, we introduce the variable y = x' and use:

(%i4) plotdf([y,-(k*x + c*y + b*x^3)/m],
             [parameters,"k=-1,m=1.0,c=0,b=1"],
             [sliders,"k=-2:2,m=-1:1"],[tstep,0.1])$

132-5

The direction field for a damped pendulum, including the solution for the given initial conditions, with a slider that can be used to change the value of the mass m, and with a plot of the two state variables as a function of time:

(%i5) plotdf([w,-g*sin(a)/l - b*w/m/l], [a,w],
             [parameters,"g=9.8,l=0.5,m=0.3,b=0.05"],
             [trajectory_at,1.05,-9],[tstep,0.01],
             [a,-10,2], [w,-14,14], [direction,forward],
             [nsteps,300], [sliders,"m=0.1:1"], [versus_t,1])$

132-6
132-7

Maxima – 131 – Funzioni numeriche – 3

Dqpf

Continuo da qui, copio dal Reference Manual, PDF scaricabile da qui, sono a p.387.

Funzioni per soluzioni numeriche di equazioni

horner (expr, x)
horner (expr)
Returns a rearranged representation of expr as in Horner’s rule, using x as the main variable if it is specified. x may be omitted in which case the main variable of the canonical rational expression form of expr is used.

horner sometimes improves stability if expr is to be numerically evaluated. It is also useful if Maxima is used to generate programs to be run in Fortran. See also stringout.

(%i1) expr: 1e-155*x^2 - 5.5*x + 5.2e155;
                                  2
(%o1)                   1.0E-155 x  - 5.5 x + 5.2E+155
(%i2) expr2: horner (%, x), keepfloat: true;
(%o2)                1.0 ((1.0E-155 x - 5.5) x + 5.2E+155)
(%i3) ev (expr, x=1e155);
(%o3)               1.0E-155 false - 2.999999999999997E+154
(%i4) ev (expr2, x=1e155);
(%o4)                       7.000000000000006E+154

find_root (expr, x, a, b, [abserr, relerr])
find_root (f, a, b, [abserr, relerr])
bf_find_root (expr, x, a, b, [abserr, relerr])
bf_find_root (f, a, b, [abserr, relerr])
find_root_error
find_root_abs
find_root_rel
Finds a root of the expression expr or the function f over the closed interval [a, b].

The expression expr may be an equation, in which case find_root seeks a root of lhs(expr) - rhs(expr).

Given that Maxima can evaluate expr or f over [a, b] and that expr or f is continuous, find_root is guaranteed to find the root, or one of the roots if there is more than one.

find_root initially applies binary search. If the function in question appears to be smooth enough, find_root applies linear interpolation instead.

bf_find_root is a bigfloat version of find_root. The function is computed using bigfloat arithmetic and a bigfloat result is returned. Otherwise, bf_find_root is identical to find_root, and the following description is equally applicable to bf_find_root.

The accuracy of find_root is governed by abserr and relerr, which are optional keyword arguments to find_root. These keyword arguments take the form key=val.

The keyword arguments are:
abserr Desired absolute error of function value at root. Default is find_root_abs.
relerr Desired relative error of root. Default is find_root_rel.

find_root stops when the function in question evaluates to something less than or equal to abserr, or if successive approximants x_0, x_1 differ by no more than relerr * max(abs(x_0), abs(x_1)). The default values of find_root_abs and find_root_rel are both zero.

find_root expects the function in question to have a different sign at the endpoints of the search interval. When the function evaluates to a number at both endpoints and these numbers have the same sign, the behavior of find_root is governed by find_root_error. When find_root_error is true, find_root prints an error message.

Otherwise find_root returns the value of find_root_error. The default value of find_root_error is true.

If f evaluates to something other than a number at any step in the search algorithm, find_root returns a partially-evaluated find_root expression.

The order of a and b is ignored; the region in which a root is sought is [min(a, b), max(a, b)].

(%i5) f(x) := sin(x) - x/2;
                                               x
(%o5)                         f(x) := sin(x) - -
                                               2
(%i6) find_root (sin(x) - x/2, x, 0.1, %pi);
(%o6)                          1.895494267033981
(%i7) find_root (sin(x) = x/2, x, 0.1, %pi);
(%o7)                          1.895494267033981
(%i8) find_root (f(x), x, 0.1, %pi);
(%o8)                          1.895494267033981
(%i9) find_root (f, 0.1, %pi);
(%o9)                          1.895494267033981
(%i10) find_root (exp(x) = y, x, 0, 100);
                                   x
(%o10)                 find_root(%e  = y, x, 0.0, 100.0)
(%i11) log (10.0);
(%o11)                         2.302585092994046
(%i12) fpprec:32;
(%o12)                                32
(%i13) bf_find_root (exp(x) = y, x, 0, 100), y = 10;
(%o13)                2.3025850929940456840179914546844b0
(%i14) log(10b0);
(%o14)                2.3025850929940456840179914546844b0

newton (expr, x, x_0, eps)
Returns an approximate solution of expr = 0 by Newton’s method, considering expr to be a function of one variable, x. The search begins with x = x_0 and proceeds until abs(expr) < eps (with expr evaluated at the current value of x).

newton allows undefined variables to appear in expr, so long as the termination test abs(expr) < eps evaluates to true or false. Thus it is not necessary that expr evaluate to a number.

load(newton1) loads this function.

See also realroots, allroots, find_root and mnewton-pkg.

(%i15) load ("newton1");
(%o15)        /usr/share/maxima/5.41.0/share/numeric/newton1.mac
(%i16) newton (cos (u), u, 1, 1/100);
(%o16)                         1.57067527716125
(%i17) ev (cos (u), u = %);
(%o17)                       1.210496333503352E-4
(%i18) ev (x^2 - a^2, x = %);
                                                  2
(%o18)                     1.46530137342506E-8 - a

Maxima – 130 – Funzioni numeriche – 2

feynmann

Continuo da qui, copio dal Reference Manual, PDF scaricabile da qui, sono a p.384.

fft (x)
Computes the complex fast Fourier transform. x is a list or array (named or unnamed) which contains the data to transform. The number of elements must be a power of 2. The elements must be literal numbers (integers, rationals, floats, or bigfloats) or symbolic constants, or expressions a + b*%i where a and b are literal numbers or symbolic constants.

fft returns a new object of the same type as x, which is not modified. Results are always computed as floats or expressions a + b*%i where a and b are floats. If bigfloat precision is needed the function bf_fft can be used instead as a drop-in replacement of fft that is slower, but supports bfloats. In addition if it is known that the input consists of only real values (no imaginary parts), real_fft can be used which is potentially faster.

The discrete Fourier transform is defined as follows. Let y be the output of the transform. Then for k from 0 through n - 1,
y[k] = (1/n) sum(x[j] exp(+2 %i %pi j k / n), j, 0, n - 1)

As there are various sign and normalization conventions possible, this definition of the transform may differ from that used by other mathematical software.

When the data x are real, real coefficients a and b can be computed such that
x[j] = sum(a[k]*cos(2*%pi*j*k/n)+b[k]*sin(2*%pi*j*k/n), k, 0, n/2)
with
a[0] = realpart (y[0])
b[0] = 0
and, for k from 1 through n/2 - 1,
a[k] = realpart (y[k] + y[n - k])
b[k] = imagpart (y[n - k] - y[k])
and
a[n/2] = realpart (y[n/2])
b[n/2] = 0

load(fft) loads this function.

See also inverse_fft (inverse transform), recttopolar, and polartorect. See real_fft for FFTs of a real-valued input, and bf_fft and bf_real_fft for operations on bigfloat values.

(%i1) load ("fft") $

(%i2) fpprintprec : 4 $

(%i3) L : [1, 2, 3, 4, -1, -2, -3, -4] $

(%i4) L1 : fft (L);
(%o4) [0.0, 1.81 %i - 0.1035, 0.0, 0.3106 %i + 0.6035, 0.0,
                                 0.6035 - 0.3106 %i, 0.0, (- 1.81 %i) - 0.1035]
(%i5) L2 : inverse_fft (L1);
(%o5) [1.0, 4.44E-16 %i + 2.0, 3.0, 4.0 - 4.44E-16 %i, - 1.0,
                               (- 4.44E-16 %i) - 2.0, - 3.0, 4.44E-16 %i - 4.0]
(%i6) lmax (abs (L2 - L));
(%o6)                              4.44E-16

Complex data.

(%i1) load ("fft") $

(%i2) fpprintprec : 4 $

(%i3) L : [1, 1 + %i, 1 - %i, -1, -1, 1 - %i, 1 + %i, 1] $

(%i4) L1 : fft (L);
(%o4) [0.5, 0.5, 0.25 %i - 0.25, (- 0.3535 %i) - 0.3535, 0.0, 0.5,
                                        (- 0.25 %i) - 0.25, 0.3535 %i + 0.3535]
(%i5) L2 : inverse_fft (L1);
(%o5) [1.0, 1.0 %i + 1.0, 1.0 - 1.0 %i, - 1.0, - 1.0, 1.0 - 1.0 %i,
                                                             1.0 %i + 1.0, 1.0]
(%i6) lmax (abs (L2 - L));
(%o6)                                 0.0

Computation of sine and cosine coefficients.

(%i1) load ("fft") $

(%i2) fpprintprec : 4 $

(%i3) L : [1, 2, 3, 4, 5, 6, 7, 8] $

(%i4) n : length (L) $

(%i5) x : make_array (any, n) $

(%i6) fillarray (x, L) $

(%i7) y : fft (x) $

(%i8) a : make_array (any, n/2 + 1) $

(%i9) b : make_array (any, n/2 + 1) $

(%i10) a[0] : realpart (y[0]) $

(%i11) b[0] : 0 $

(%i12) for k : 1 thru n/2 - 1 do
         (a[k] : realpart (y[k] + y[n - k]),
         b[k] : imagpart (y[n - k] - y[k]));
(%o12)                               done
(%i13) a[n/2] : y[n/2] $

(%i14) b[n/2] : 0 $

(%i15) listarray (a);
(%o15)                 [4.5, - 1.0, - 1.0, - 1.0, - 0.5]
(%i16) listarray (b);
(%o16)                    [0, 2.414, 1.0, 0.4142, 0]
(%i17) f(j) := sum (a[k]*cos(2*%pi*j*k/n) + b[k]*sin(2*%pi*j*k/n), k, 0, n/2) $

(%i18) makelist (float (f (j)), j, 0, n - 1);
(%o18)             [1.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0]

real_fft (x)
Computes the fast Fourier transform of a real-valued sequence x. This is equivalent to performing fft(x), except that only the first N/2+1 results are returned, where N is the length of x. N must be power of two.

No check is made that x contains only real values.

The symmetry properites of the Fourier transform of real sequences to reduce he complexity. In particular the first and last output values of real_fft are purely real.

For larger sequencyes, real_fft may be computed more quickly than fft.

Since the output length is short, the normal inverse_fft cannot be directly used. Use inverse_real_fft to compute the inverse.

inverse_real_fft (y)
Computes the inverse Fourier transfrom of y, which must have a length of N/2+1 where N is a power of two. That is, the input x is expected to be the output of real_fft.

No check is made to ensure that the input has the correct format. (The first and last
elements must be purely real.)

bf_inverse_fft (y)
Computes the inverse complex fast Fourier transform. This is the bigfloat version of inverse_fft that converts the input to bigfloats and returns a bigfloat result.

bf_fft (y)
Computes the forward complex fast Fourier transform. This is the bigfloat version of fft that converts the input to bigfloats and returns a bigfloat result.

bf_real_fft (x)
Computes the forward fast Fourier transform of a real-valued input returning a bigfloat result. This is the bigfloat version of real_fft.

bf_inverse_real_fft (y)
Computes the inverse fast Fourier transfrom with a real-valued bigfloat output. This is the bigfloat version of inverse_real_fft.