SymPy – 26 – conclusioni

Continuo da qui.

L’esame di SymPy finisce qui. La documentazione continua, anzi entra nel vivo dopo il tutorial che ho esaminato. Ma (il bello di essere ormai fuori dal mondo del lavoro) ci sono altre cose –über-sexy– da seguire, prossimamente.

SymPy è OK. Essere integrato con Python, anzi essere Python, è un punto a suo favore (anzi di più!). Essere gestibile all’interno della REPL è un’altro punto. Per quel che ho visto conviene usare la REPL di Python o quella di IPython (come faccio io) e non quella integrata nel Web disponibile sul sito di SymPy: non è aggiornata e soggetta a rallentamenti non dovuti alla connessione.

Uno dei motivi che mi hanno portato alla scoperta di SymPy è la possibilità di visualizzare i risultati in formato LaTeX; quasi, c’è bisogno di qualche edizione dell’output.

Occorre usare l’output di print().

Si può inserire in WordPress in questa forma:

f{\left (x \right )} - 2 \frac{d}{d x} f{\left (x \right )} + 
\frac{d^{2}}{d x^{2}} f{\left (x \right )} = \sin{\left (x \right )}

tra i tag $latex#formula#$ e appare così:

In alternativa si può utilizzare LaTeX Equation Editor.

Si può quindi scaricare l’immagine, cosa che faccio abitualmente (non sono tanto bravo con LaTeX), dopo aver modificato le opzioni come indicato in figura.

Un’alternativa (über-sexy, anzi di più) è Maxima (e wxMaxima, o GMaxima, o …); ne ho accennato in passato, sarebbe un’altra cosa da approfondire. Avendo tempo.

🤢

SymPy – 25 – manipolazione avanzata di espressioni – 3

Continuo da qui, copio qui.

Percorrere l’albero (tree)
With this knowledge [post precedente], let’s look at how we can recurse through an expression tree. The nested nature of args is a perfect fit for recursive functions. The base case will be empty args. Let’s write a simple function that goes through an expression and prints all the args at each level.

See how nice it is that () signals leaves in the expression tree. We don’t even have to write a base case for our recursion; it is handled automatically by the for loop.

Let’s test our function.

Can you guess why we called our function pre? We just wrote a pre-order traversal function for our expression tree. See if you can write a post-order traversal function.

Such traversals are so common in SymPy that the generator functions preorder_traversal and postorder_traversal are provided to make such traversals easy. We could have also written our algorithm as

Prevenire la valutazione dell’espressione
There are generally two ways to prevent the evaluation, either pass an evaluate=False parameter while constructing the expression, or create an evaluation stopper by wrapping the expression with UnevaluatedExpr.

If you don’t remember the class corresponding to the expression you want to build (operator overloading usually assumes evaluate=True), just use sympify and pass a string:

Note that evaluate=False won’t prevent future evaluation in later usages of the expression:

That’s why the class UnevaluatedExpr comes handy. UnevaluatedExpr is a method provided by SymPy which lets the user keep an expression unevaluated. By unevaluated it is meant that the value inside of it will not interact with the expressions outside of it to give simplified outputs. For example:

The x remaining alone is the x wrapped by UnevaluatedExpr. To release it:

Other examples:

A point to be noted is that UnevaluatedExpr cannot prevent the evaluation of an expression which is given as argument. For example:

Remember that expr2 will be evaluated if included into another expression. Combine both of the methods to prevent both inside and outside evaluations:

UnevalutedExpr is supported by SymPy printers and can be used to print the result in different output forms. For example

In order to release the expression and get the evaluated LaTeX form, just use doit():

Uhmmmm… un po’ pasticciato ma OK 😊

SymPy – 24 – manipolazione avanzata di espressioni – 2

Continuo da qui, copio qui.

Ricorsione attraverso l’albero delle espressioni
Now that you know how expression trees work in SymPy, let’s look at how to dig our way through an expression tree. Every object in SymPy has two very important attributes, func, and args.

func
func is the head of the object. For example, (x*y).func is Mul. Usually it is the same as the class of the object (though there are exceptions to this rule).

Two notes about func. First, the class of an object need not be the same as the one used to create it. For example

We created Add(x, x), so we might expect expr.func to be Add, but instead we got Mul. Why is that? Let’s take a closer look at expr.

Add(x, x), i.e., x + x, was automatically converted into Mul(2, x), i.e., 2*x, which is a Mul. SymPy classes make heavy use of the __new__ class constructor, which, unlike __init__, allows a different class to be returned from the constructor.

Second, some classes are special-cased, usually for efficiency reasons [note 3].

Note 3: Classes like One and Zero are singletonized, meaning that only one object is ever created, no matter how many times the class is called. This is done for space efficiency, as these classes are very common. For example, Zero might occur very often in a sparse matrix represented densely. As we have seen, NegativeOne occurs any time we have -x or 1/x. It is also done for speed efficiency because singletonized objects can be compared by is. The unique objects for each singletonized class can be accessed from the S object.

For the most part, these issues will not bother us. The special classes Zero, One, NegativeOne, and so on are subclasses of Integer, so as long as you use isinstance, it will not be an issue.

args
args are the top-level arguments of the object. (x*y).args would be (x, y). Let’s look at some examples

From this, we can see that expr == Mul(3, y**2, x). In fact, we can see that we can completely reconstruct expr from its func and its args.

Note that although we entered 3*y**2*x, the args are (3, x, y**2). In a Mul, the Rational coefficient will come first in the args, but other than that, the order of everything else follows no special pattern. To be sure, though, there is an order.

Mul’s args are sorted, so that the same Mul will have the same args. But the sorting is based on some criteria designed to make the sorting unique and efficient that has no mathematical significance.

The srepr form of our expr is Mul(3, x, Pow(y, 2)). What if we want to get at the args of Pow(y, 2). Notice that the y**2 is in the third slot of expr.args, i.e., expr.args[2].

So to get the args of this, we call expr.args[2].args.

Now what if we try to go deeper. What are the args of y. Or 2. Let’s see.

They both have empty args. In SymPy, empty args signal that we have hit a leaf of the expression tree.

So there are two possibilities for a SymPy expression. Either it has empty args, in which case it is a leaf node in any expression tree, or it has args, in which case, it is a branch node of any expression tree. When it has args, it can be completely rebuilt from its func and its args. This is expressed in the key invariant.

Key invariant
Every well-formed SymPy expression must either have empty args or satisfy expr == expr.func(*expr.args).

(Recall that in Python if a is a tuple, then f(*a) means to call f with arguments from the elements of a, e.g., f(*(1, 2, 3)) is the same as f(1, 2, 3).)

This key invariant allows us to write simple algorithms that walk expression trees, change them, and rebuild them into new expressions. Tutto questo nel prossimo post.

SymPy – 23 – manipolazione avanzata di espressioni – 1

Continuo da qui, copio qui.

In this section, we discuss some ways that we can perform advanced manipulation of expressions.

Capire l’albero delle espressioni
Before we can do this, we need to understand how expressions are represented in SymPy. A mathematical expression is represented as a tree. Let us take the expression 2x+xy, i.e., 2**x + x*y. We can see what this expression looks like internally by using srepr.

Quick Tip: To play with the srepr form of expressions in the SymPy Live shell, change the output format to Repr in the settings.

The easiest way to tear this apart is to look at a diagram of the expression tree:

Note: The above diagram was made using Graphviz and the dotprint function.

First, let’s look at the leaves of this tree. Symbols are instances of the class Symbol. While we have been doing

x = symbols('x')

we could have also done

x = Symbol('x')

Either way, we get a Symbol with the name “x” [We have been using symbols instead of Symbol because it automatically splits apart strings into multiple Symbols. symbols('x y z') returns a tuple of three Symbols. Symbol('x y z') returns a single Symbol called x y z]. For the number in the expression, 2, we got Integer(2). Integer is the SymPy class for integers. It is similar to the Python built-in type int, except that Integer plays nicely with other SymPy types.

When we write 2**x, this creates a Pow object. Pow is short for “power”.

We could have created the same object by calling Pow(2, x)

Note that in the srepr output, we see Integer(2), the SymPy version of integers, even though technically, we input 2, a Python int. In general, whenever you combine a SymPy object with a non-SymPy object via some function or operation, the non-SymPy object will be converted into a SymPy object. The function that does this is sympify [Technically, it is an internal function called sympify, which differs from sympify in that it does not convert strings. x + '2' is not allowed].

We have seen that 2**x is represented as Pow(2, x). What about x*y? As we might expect, this is the multiplication of x and y. The SymPy class for multiplication is Mul.

Thus, we could have created the same object by writing Mul(x, y).

Now we get to our final expression, 2**x + x*y. This is the addition of our last two objects, Pow(2, x), and Mul(x, y). The SymPy class for addition is Add, so, as you might expect, to create this object, we use Add(Pow(2, x), Mul(x, y)).

SymPy expression trees can have many branches, and can be quite deep or quite broad. Here is a more complicated example

Here is a diagram

This expression reveals some interesting things about SymPy expression trees. Let’s go through them one by one.

Let’s first look at the term x**2. As we expected, we see Pow(x, 2). One level up, we see we have Mul(-1, Pow(x, 2)). There is no subtraction class in SymPy. x - y is represented as x + -y, or, more completely, x + -1*y, i.e., Add(x, Mul(-1, y)).

Next, look at 1/y. We might expect to see something like Div(1, y), but similar to subtraction, there is no class in SymPy for division. Rather, division is represented by a power of -1. Hence, we have Pow(y, -1). What if we had divided something other than 1 by y, like x/y? Let’s see.

We see that x/y is represented as x*y**-1, i.e., Mul(x, Pow(y, -1)).

Finally, let’s look at the sin(x*y)/2 term. Following the pattern of the previous example, we might expect to see Mul(sin(x*y), Pow(Integer(2), -1)). But instead, we have Mul(Rational(1, 2), sin(x*y)). Rational numbers are always combined into a single term in a multiplication, so that when we divide by 2, it is represented as multiplying by 1/2.

Finally, one last note. You may have noticed that the order we entered our expression and the order that it came out from srepr or in the graph were different. You may have also noticed this phenomenon earlier in the tutorial. For example

This because in SymPy, the arguments of the commutative operations Add and Mul are stored in an arbitrary (but consistent!) order, which is independent of the order inputted (if you’re worried about noncommutative multiplication, don’t be. In SymPy, you can create noncommutative Symbols using Symbol('A', commutative=False), and the order of multiplication for noncommutative Symbols is kept the same as the input). Furthermore, as we shall see in the next section, the printing order and the order in which things are stored internally need not be the same either.

In general, an important thing to keep in mind when working with SymPy expression trees is this: the internal representation of an expression and the way it is printed need not be the same. The same is true for the input form. If some expression manipulation algorithm is not working in the way you expected it to, chances are, the internal representation of the object is different from what you thought it was.

Quick Tip: The way an expression is represented internally and the way it is printed are often not the same.

OK (quasi) panico 😯

SymPy – 22 – matrici – 2

Continuo da qui, copio qui.

Metodi avanzati
determinante
To compute the determinant of a matrix, use det.

reduced row echelon form (RREF)
To put a matrix into reduced row echelon form, use rref. rref returns a tuple of two elements. The first is the reduced row echelon form, and the second is a tuple of indices of the pivot columns.

Note: The first element of the tuple returned by rref is of type Matrix. The second is of type list.

nullspace
To find the nullspace of a matrix, use nullspace. nullspace returns a list of column vectors that span the nullspace of the matrix.

columnspace
To find the columnspace of a matrix, use columnspace. columnspace returns a list of column vectors that span the columnspace of the matrix.

aurovalori, autovettori e diagonalizzazione
To find the eigenvalues of a matrix, use eigenvals. eigenvals returns a dictionary of eigenvalue:algebraic multiplicity pairs (similar to the output of roots).

This means that M has eigenvalues -2, 3, and 5, and that the eigenvalues -2 and 3 have algebraic multiplicity 1 and that the eigenvalue 5 has algebraic multiplicity 2.

To find the eigenvectors of a matrix, use eigenvects. eigenvects returns a list of tuples of the form (eigenvalue:algebraic multiplicity, [eigenvectors]).

This shows us that, for example, the eigenvalue 5 also has geometric multiplicity 2, because it has two eigenvectors. Because the algebraic and geometric multiplicities are the same for all the eigenvalues, M is diagonalizable.

To diagonalize a matrix, use diagonalize. diagonalize returns a tuple (P,D), where D is diagonal and M=PDP−1.

Note that since eigenvects also includes the eigenvalues, you should use it instead of eigenvals if you also want the eigenvectors. However, as computing the eigenvectors may often be costly, eigenvals should be preferred if you only wish to find the eigenvalues.

If all you want is the characteristic polynomial, use charpoly. This is more efficient than eigenvals, because sometimes symbolic roots can be expensive to calculate.

Quick Tip: lambda is a reserved keyword in Python, so to create a Symbol called λ, while using the same names for SymPy Symbols and Python variables, use lamda (without the b). It will still pretty print as λ. (Unicode per λ: 955 0x3bb).

SymPy – 21 – matrici – 1

Continuo da qui, copio qui.

To make a matrix in SymPy, use the Matrix object. A matrix is constructed by providing a list of row vectors that make up the matrix. For example, to construct the matrix

use

To make it easy to make column vectors, a list of elements is considered to be a column vector.

Matrices are manipulated just like any other object in SymPy or Python.

One important thing to note about SymPy matrices is that, unlike every other object in SymPy, they are mutable. This means that they can be modified in place, as we will see below. The downside to this is that Matrix cannot be used in places that require immutability, such as inside other SymPy expressions or as keys to dictionaries. If you need an immutable version of Matrix, use ImmutableMatrix.

Operazioni di base
shape
Here are some basic operations on Matrix. To get the shape of a matrix use shape

accedere a righe e colonne
To get an individual row or column of a matrix, use row or col. For example, M.row(0) will get the first row. M.col(-1) will get the last column.

cancellare e inserire righe e colonne
To delete a row or column, use row_del or col_del. These operations will modify the Matrix in place.

To insert rows or columns, use row_insert or col_insert. These operations do not operate in place.

Unless explicitly stated, the methods mentioned below do not operate in place. In general, a method that does not operate in place will return a new Matrix and a method that does operate in place will return None.

Metodi base
As noted above, simple operations like addition and multiplication are done just by using +, *, and **. To find the inverse of a matrix, just raise it to the -1 power.

To take the transpose of a Matrix, use T.

Costruttori
Several constructors exist for creating common matrices. To create an identity matrix, use eye. eye(n) will create an n×n identity matrix.

To create a matrix of all zeros, use zeros. zeros(n, m) creates an n×m matrix of 0s.

Similarly, ones creates a matrix of ones.

To create diagonal matrices, use diag. The arguments to diag can be either numbers or matrices. A number is interpreted as a 1×1 matrix. The matrices are stacked diagonally. The remaining elements are filled with 0s.

SymPy – 20 – risolutori – 2

Monument Valley – Masterson

Continuo da qui, copio qui.

Risolvere equazioni differenziali
To solve differential equations, use dsolve. First, create an undefined function by passing cls=Function to the symbols function.

f and g are now undefined functions. We can call f(x), and it will represent an unknown function.

Derivatives of f(x) are unevaluated.

(see the Derivatives [prossimamente] section for more on derivatives).

To represent the differential equation

we would thus use

To solve the ODE, pass it and the function to solve for to dsolve.

dsolve returns an instance of Eq. This is because in general, solutions to differential equations cannot be solved explicitly for the function.

The arbitrary constants in the solutions from dsolve are symbols of the form C1, C2, C3, and so on.

SymPy – 19 – risolutori – 1

Continuo da qui, copio qui.

Una nota a riguardo le equazioni
Recall from the gotchas [qui] section of this tutorial that symbolic equations in SymPy are not represented by = or ==, but by Eq.

However, there is an even easier way. In SymPy, any expression not in an Eq is automatically assumed to equal 0 by the solving functions. Since a=b if and only if a−b=0, this means that instead of using x == y, you can just use x - y. For example

This is particularly useful if the equation you wish to solve is already equal to 0. Instead of typing solveset(Eq(expr, 0), x), you can just use solveset(expr, x).

Risolvere algebricamente equazioni
The main function for solving algebraic equations is solveset. The syntax for solveset is solveset(equation, variable=None, domain=S.Complexes) Where equations may be in the form of Eq instances or expressions that are assumed to be equal to zero.

Please note that there is another function called solve which can also be used to solve equations. The syntax is solve(equations, variables) However, it is recommended to use solveset instead.

When solving a single equation, the output of solveset is a FiniteSet or an Interval or ImageSet of the solutions.

If there are no solutions, an EmptySet is returned and if it is not able to find solutions then a ConditionSet is returned.

In the solveset module, the linear system of equations is solved using linsolve. In future we would be able to use linsolve directly from solveset. Following is an example of the syntax of linsolve.

1. List of Equations Form:

2. Augmented Matrix Form:

3. A*x = b Form

Note: The order of solution corresponds the order of given symbols.

In the solveset module, the non linear system of equations is solved using nonlinsolve. Following are examples of nonlinsolve.

No, nonlinsolve risulta non definita, da indagare.
Uh! risolto, con l’aggiornamento di versione come raccontato nel post precedente; si continua 😊

In the solveset module, the non linear system of equations is solved using nonlinsolve. Following are examples of nonlinsolve.

1. When only real solution is present:

2. When only complex solution is present:

3. When both real and complex solution is present:

4. If non linear system of equations is Positive dimensional system (A system with infinitely many solutions is said to be positive-dimensional):

Note:

1. The order of solution corresponds the order of given symbols.

2. Currently nonlinsolve doesn’t return solution in form of LambertW (if there is solution present in the form of LambertW). solve can be used for such cases(not all solution):

3. Currently nonlinsolve is not properly capable of solving the system of equations having trigonometric functions. solve can be used for such cases(not all solution):

solveset reports each solution only once. To get the solutions of a polynomial including multiplicity use roots.

The output {0: 1, 3: 2} of roots means that 0 is a root of multiplicity 1 and 3 is a root of multiplicity 2.

Currently solveset is not capable of solving the equations solvable by LambertW (Transcendental equation solver). solve can be used for such cases:

SymPy – 18 – calcolo infinitesimale – 4

Continuo da qui, copio qui.

Ho avuto un po’ di problemi legati alla versione non aggiornata; chissà se ora è OK 😊
L’aggiornamento l’ho effettuato via Anaconda, seguendo le indicazioni qui.

Inoltre riconsidero anche l’uso della Sympy Live Shell: è lenta e torno a IPython 😊

Differenze finite
So far we have looked at expressions with analytic derivatives and primitive functions respectively. But what if we want to have an expression to estimate a derivative of a curve for which we lack a closed form representation, or for which we don’t know the functional values for yet. One approach would be to use a finite difference approach.

The simplest way the differentiate using finite differences is to use the differentiate_finite function:

If we want to expand the intermediate derivative we may pass the flag evaluate=True:

This form however does not respect the product rule.

If you already have a Derivative instance, you can use the as_finite_difference method to generate approximations of the derivative to arbitrary order:

here the first order derivative was approximated around x using a minimum number of points (2 for 1st order derivative) evaluated equidistantly using a step-size of 1. We can use arbitrary steps (possibly containing symbolic expressions):

If you are just interested in evaluating the weights, you can do so manually:

note that we only need the last element in the last sublist returned from finite_diff_weights. The reason for this is that the function also generates weights for lower derivatives and using fewer points (see the documentation of finite_diff_weights for more details).

If using finite_diff_weights directly looks complicated, and the as_finite_difference method of Derivative instances is not flexible enough, you can use apply_finite_diff which takes order, x_list, y_list and x0 as parameters:

SymPy – 17 – calcolo infinitesimale – 3

Continuo da qui, copio qui.

Limiti
SymPy can compute symbolic limits with the limit function. The syntax to compute

is limit(f(x), x, x0).

limit should be used instead of subs whenever the point of evaluation is a singularity. Even though SymPy has objects to represent ∞, using them for evaluation is not reliable because they do not keep track of things like rate of growth. Also, things like ∞-∞ and ∞/∞ return NaN (not-a-number). For example

Like Derivative and Integral, limit has an unevaluated counterpart, Limit. To evaluate it, use doit.

To evaluate a limit at one side only, pass '+' or '-' as a third argument to limit. For example, to compute

do

Espansione di serie
SymPy can compute asymptotic series expansions of functions around a point. To compute the expansion of f(x) around the point x=x0 terms of order xn, use f(x).series(x, x0, n). x0 and n can be omitted, in which case the defaults x0=0 and n=6 will be used.

The O(x4) term at the end represents the Landau order term at x=0 (not to be confused with big O notation used in computer science, which generally represents the Landau order term at x=∞). It means that all x terms with power greater than or equal to x4 are omitted. Order terms can be created and manipulated outside of series. They automatically absorb higher order terms.

If you do not want the order term, use the removeO method.

The O notation supports arbitrary limit points (other than 0):

Aggiornamento: il capitolo continua con Finite differences ma semplicemente non funziona 😡
Le funzioni risultano non definite; Stack Overflow ha diversi post a riguardo, tutti segnalantio problemi. Salto 😜 Forse in futuro… 😉

Aggiornamento 2: era solo la versione da aggiornare; non leggete l’aggiornamento precedente, quello scancellato 😉