Tag Archives: Maths

Playing with Graphs

1 Reply

This twenty-fourth article of the mathematical journey through open source, demonstrates the concepts of graph theory for playing with graphs using the graphs package of Maxima.

<< Twenty-third Article

In our previous article, we familiarized ourselves with the notion of simple graphs, and how graphs package of Maxima allows us to create and visualize them. Assuming that knowledge, in this article, we are going to play with graphs and their properties, using the functionalities provided by Maxima’s graphs package.

Graph modifications

We have already created various graphs with create_graph() and make_graph() functions of the graphs package of Maxima. What if we want to modify some existing graphs, say by adding or removing some edges or vertices? Exactly for such operations, Maxima provides the following functions:

  • add_edge(<edge>, <g>) – Adds <edge> into the graph <g>
  • add_edges(<edge_list>, <g>) – Adds edges specified by <edge_list> into the graph <g>
  • add_vertex(<vertex>, <g>) – Adds <vertex> into the graph <g>
  • add_vertices(<vertex_list>, <g>) – Adds vertices specified by <vertex_list> into the graph <g>
  • connect_vertices(<u_list>, <v_list>, <g>) – Connects all vertices from <u_list> to all vertices in <v_list> in the graph <g>
  • contract_edge(<edge>, <g>) – Merges the vertices of the <edge> and the edges incident on those vertices, in the graph <g>
  • remove_edge(<edge>, <g>) – Removes the <edge> from the graph <g>
  • remove_vertex(<vertex>, <g>) – Removes the <vertex> and the associated edges from the graph <g>

Some of the above functions are demonstrated below:

$ maxima -q
(%i1) load(graphs)$ /* Loading the graphs package */
...
0 errors, 0 warnings
(%i2) g: create_graph(4, [[0, 1], [0, 2]]);
(%o2)                     GRAPH(4 vertices, 2 edges)
(%i3) print_graph(g)$

Graph on 4 vertices with 2 edges.
Adjacencies:
  3 :
  2 :  0
  1 :  0
  0 :  2  1
(%i4) add_edge([1, 2], g)$
(%i5) print_graph(g)$

Graph on 4 vertices with 3 edges.
Adjacencies:
  3 :
  2 :  1  0
  1 :  2  0
  0 :  2  1
(%i6) contract_edge([0, 1], g)$
(%i7) print_graph(g)$

Graph on 3 vertices with 1 edges.
Adjacencies:
  3 :
  2 :  0
  0 :  2

In the above examples, if we do not intend to modify the original graph, we may make a copy of it using copy_graph(), and then operate on the copy, as follows:

(%i8) h: copy_graph(g);
(%o8)                     GRAPH(3 vertices, 1 edges)
(%i9) add_vertex(1, h)$
(%i10) print_graph(h)$

Graph on 4 vertices with 1 edges.
Adjacencies:
  1 :
  0 :  2
  2 :  0
  3 :
(%i11) print_graph(g)$ /* Notice g is unmodified */

Graph on 3 vertices with 1 edges.
Adjacencies:
  3 :
  2 :  0
  0 :  2
(%i12) quit();

Advanced graph creations

New graphs can also be created based on existing graphs and their properties by various interesting operations. Few of them are:

  • underlying_graph(<dg>) – Returns the underlying graph of the directed graph <dg>
  • complement_graph(<g>) – Returns the complement graph of graph <g>
  • line_graph(<g>) – Returns a graph that represents the adjacencies between the edges of graph <g>
  • graph_union(<g1>, <g2>) – Returns a graph with edges and vertices of both graphs <g1> and <g2>
  • graph_product(<g1>, <g2>) – Returns the Cartesian product of graphs <g1> and <g2>

Here are some examples to demonstrate the simpler functions:

$ maxima -q
(%i1) load(graphs)$
...
0 errors, 0 warnings
(%i2) g: create_graph(4, [[0, 1], [0, 2], [0, 3]], directed = true);
(%o2)                     DIGRAPH(4 vertices, 3 arcs)
(%i3) print_graph(g)$

Digraph on 4 vertices with 3 arcs.
Adjacencies:
  3 :
  2 :
  1 :
  0 :  3  2  1
(%i4) h: underlying_graph(g);
(%o4)                     GRAPH(4 vertices, 3 edges)
(%i5) print_graph(h)$

Graph on 4 vertices with 3 edges.
Adjacencies:
  0 :  1  2  3
  1 :  0
  2 :  0
  3 :  0
(%i6) print_graph(complement_graph(h))$

Graph on 4 vertices with 3 edges.
Adjacencies:
  3 :  2  1
  2 :  3  1
  1 :  3  2
  0 :
(%i7) print_graph(graph_union(h, complement_graph(h)))$

Graph on 8 vertices with 6 edges.
Adjacencies:
  4 :
  5 :  6  7
  6 :  5  7
  7 :  5  6
  3 :  0
  2 :  0
  1 :  0
  0 :  3  2  1
(%i8) quit();

Basic graph properties

graph_order(<g>), vertices(<g>) returns the number of vertices and the list of vertices, respectively, in the graph <g>. graph_size(<g>), edges(<g>) returns the number of edges and the list of edges, respectively, in the graph <g>.

A graph is a collection of vertices and edges. Hence, most of its properties are centred around them. The following are graph related predicates provided by the graphs package of Maxima:

  • is_graph(<g>) – returns true if <g> is a graph, false otherwise
  • is_digraph(<g>) – returns true if <g> is a directed graph, false otherwise
  • is_graph_or_digraph(<g>) – returns true if <g> is a graph or a directed graph, false otherwise
  • is_connected(<g>) – returns true if graph <g> is connected, false otherwise
  • is_planar(<g>) – returns true if graph <g> can be placed on a plane without its edges crossing over each other, false otherwise
  • is_tree(<g>) – returns true if graph <g> has no simple cycles, false otherwise
  • is_biconnected(<g>) – returns true if graph <g> will remain connected even after removal of any one its vertices and the edges incident on that vertex, false otherwise
  • is_bipartite(<g>) – returns true if graph <g> is bipartite, i.e. 2-colourable, false otherwise
  • is_isomorphic(<g1>, <g2>) – returns true if graphs <g1> and <g2> have the same number of vertices and are connected in the same way, false otherwise. And, isomorphism(<g1>, <g2>) returns an isomorphism (that is a one-to-one onto mapping) between the graphs <g1> and <g2>, if it exists.
  • is_edge_in_graph(<edge>, <g>) – returns true if <edge> is in graph <g>, false otherwise
  • is_vertex_in_graph(<vertex>, <g>) – returns true if <vertex> is in graph <g>, false otherwise

The following example specifically demonstrates the isomorphism property, from the list above:

$ maxima -q
(%i1) load(graphs)$
...
0 errors, 0 warnings
(%i2) g1: create_graph(3, [[0, 1], [0, 2]]);
(%o2)                     GRAPH(3 vertices, 2 edges)
(%i3) g2: create_graph(3, [[1, 2], [0, 2]]);
(%o3)                     GRAPH(3 vertices, 2 edges)
(%i4) is_isomorphic(g1, g2);
(%o4)                                true
(%i5) isomorphism(g1, g2);
(%o5)                      [2 -> 0, 1 -> 1, 0 -> 2]
(%i6) quit();

Graph neighbourhoods

Lot of properties of graphs are to do with vertex and edge neighbourhoods, also referred as adjacencies.

For example, a graph itself could be represented by an adjacency list or matrix, which specifies the vertices adjacent to the various vertices in the graph. adjacency_matrix(<g>) returns the adjacency matrix of the graph <g>.

Number of edges incident on a vertex is called the valency or degree of the vertex, and could be obtained using vertex_degree(<v>, <g>). degree_sequence(<g>) returns the list of degrees of all the vertices of the graph <g>. In case of a directed graph, the degrees could be segregated as in-degree and out-degree, as per the edges incident into and out of the vertex, respectively. vertex_in_degree(<v>, <dg>) and vertex_out_degree(<v>, <dg>) respectively return the in-degree and out-degree for the vertex <v> of the directed graph <dg>.

neighbors(<v>, <g>), in_neighbors(<v>, <dg>), out_neighbors(<v>, <dg>) respectively return the list of adjacent vertices, adjacent in-vertices, adjacent out-vertices of the vertex <v> in the corresponding graphs.

average_degree(<g>) computes the average of degrees of all the vertices of the graph <g>. max_degree(<g>) finds the maximal degree of vertices of the graph <g>, and returns one such vertex alongwith. min_degree(<g>) finds the minimal degree of vertices of the graph <g>, and returns one such vertex alongwith.

Here follows a neighbourhood related demonstration:

$ maxima -q
(%i1) load(graphs)$
...
0 errors, 0 warnings
(%i2) g: create_graph(4, [[0, 1], [0, 2], [0, 3], [1, 2]]);
(%o2)                     GRAPH(4 vertices, 4 edges)
(%i3) string(adjacency_matrix(g)); /* string for output in single line */
(%o3)           matrix([0,0,0,1],[0,0,1,1],[0,1,0,1],[1,1,1,0])
(%i4) degree_sequence(g);
(%o4)                            [1, 2, 2, 3]
(%i5) average_degree(g);
(%o5)                                  2
(%i6) neighbors(0, g);
(%o6)                              [3, 2, 1]
(%i7) quit();

Graph connectivity

A graph is ultimately about connections, and hence lots of graph properties are centred around connectivity.

vertex_connectivity(<g>) returns the minimum number of vertices that need to be removed from the graph <g> to make the graph <g> disconnected. Similarly, edge_connectivity(<g>) returns the minimum number of edges that need to be removed from the graph <g> to make the graph <g> disconnected.

vertex_distance(<u>, <v>, <g>) returns the length of the shortest path between the vertices <u> and <v> in the graph <g>. The actual path could be obtained using shortest_path(<u>, <v>, <g>).

girth(<g>) returns the length of the shortest cycle in graph <g>.

vertex_eccentricity(<v>, <g>) returns the maximum of the vertex distances of vertex <v> with any other vertex in the connected graph <g>.

diameter(<g>) returns the maximum of the vertex eccentricities of all the vertices in the connected graph <g>.

radius(<g>) returns the minimum of the vertex eccentricities of all the vertices in the connected graph <g>.

graph_center(<g>) returns the list of vertices that have eccentricities equal to the radius of the connected graph <g>.

graph_periphery(<g>) is the list of vertices that have eccentricities equal to the diameter of the connected graph.

A minimal connectivity related demonstration for the graph shown in Figure 29 follows:

$ maxima -q
(%i1) load(graphs)$
...
0 errors, 0 warnings
(%i2) g: create_graph(9, [[0, 1], [0, 2], [1, 8], [8, 3], [2, 3], [3, 4], [4, 5],
        [3, 6], [3, 7]]);
(%o2)                     GRAPH(9 vertices, 9 edges)
(%i3) vertex_connectivity(g);
(%o3)                                  1
(%i4) edge_connectivity(g);
(%o4)                                 1
(%i5) shortest_path(0, 7, g);
(%o5)                           [0, 2, 3, 7]
(%i6) vertex_distance(0, 7, g);
(%o6)                                 3
(%i7) girth(g);
(%o7)                                 5
(%i8) diameter(g);
(%o8)                                 4
(%i9) radius(g);
(%o9)                                 2
(%i10) graph_center(g);
(%o10)                                [3]
(%i11) graph_periphery(g);
(%o11)                             [5, 1, 0]
Figure 29: Graph connectivities

Figure 29: Graph connectivities

Vertex 3 is the only centre of the graph and 0, 1, 5 are the peripheral vertices of the graph.

Graph colouring

Graph colouring has been a fascinating topic in graph theory, since its inception. It is all about colouring the vertices or edges of a graph in such a way that no adjacent elements (vertex or edge) have the same colour.

Smallest number of colours needed to colour the vertices of a graph, such that no two adjacent vertices have the same colour, is called the chromatic number of the graph. chromatic_number(<g>) computes the same. vertex_coloring(<g>) returns a list representing the colouring of the vertices of <g>, along with the chromatic number.

Smallest number of colours needed to colour the edges of a graph, such that no two adjacent edges have the same colour, is called the chromatic index of the graph. chromatic_index(<g>) computes the same. edge_coloring(<g>) returns a list representing the colouring of the edges of <g>, along with the chromatic index.

The following demonstration continues colouring the graph from the above demonstration:

(%i12) chromatic_number(g);
(%o12)                                 3
(%i13) vc: vertex_coloring(g);
(%o13) [3, [[0, 3], [1, 1], [2, 2], [3, 1], [4, 2], [5, 1], [6, 2], [7, 2], [8, 2]]]
(%i14) chromatic_index(g);
(%o14)                                 5
(%i15) ec: edge_coloring(g);
(%o15) [5, [[[0, 1], 1], [[0, 2], 2], [[1, 8], 2], [[3, 8], 5], [[2, 3], 1], 
                           [[3, 4], 2], [[4, 5], 1], [[3, 6], 3], [[3, 7], 4]]]
(%i16) draw_graph(g, vertex_coloring=vc, edge_coloring=ec, vertex_size=5,
        edge_width=3, show_id=true)$
(%i17) quit();

Figure 30 shows the coloured version of the graph, as obtained by %i16.

Figure 30: Graph colouring

Figure 30: Graph colouring

Bon voyage

With this article, we have completed a 2 year long mathematical journey through open source, starting from mathematics in Shell, covering Bench Calculator, Octave, and finally concluding with Maxima. I take this opportunity to thank my readers and wish them bon voyage with whatever they have gained through our interactions. However, this is not the end – get set for our next odyssey.

   Send article as PDF   

Visualizing Graph Theory

2 Replies

This twenty-third article of the mathematical journey through open source, introduces graph theory with visuals using the graphs package of Maxima.

<< Twenty-second Article

Graphs here refer to the structures formed using some points (or vertices), and some lines (or edges) connecting them. A simple example would be a graph with two vertices, say ‘0’ & ‘1’, and an edge between ‘0’ and ‘1’. If all the edges of a graph have a sense of direction from one vertex to another, typically represented by an arrow at their end(s), we call that a directed graph with directed edges. In such a case, we consider the edges to be not between two vertices but from one vertex to another vertex. Directed edges are also referred to as arcs. In the above example, we could have two directed arcs, one from ‘0’ to ‘1’, and another from ‘1’ to ‘0’. Figures 24 & 25 show an undirected and a directed graph, respectively.

Figure 24: Simple undirected graph

Figure 24: Simple undirected graph

Figure 25: Simple directed graph

Figure 25: Simple directed graph

Graph creation & visualization

Now, how did we get those figures? Using the graphs package of Maxima, which is loaded by invoking load(graphs) at the Maxima prompt. In the package, vertices are represented by whole numbers 0, 1, … and an edge is represented as a list of its two vertices. Using these notations, we can create graphs using the create_graph(<vertex_list>, <edge_list>[, directed]) function. And, then draw_graph(<graph>[, <option1>, <option2>, …]) would draw the graph pictorially. Code snippet below shows the same in action:

$ maxima -q
(%i1) load(graphs)$
...
0 errors, 0 warnings
(%i2) g: create_graph([0, 1], [[0, 1]])$
(%i3) dg: create_graph([0, 1], [[0, 1]], directed=true)$
(%i4) draw_graph(g, show_id=true, vertex_size=5, vertex_color=yellow)$
(%i5) draw_graph(dg, show_id=true, vertex_size=5, vertex_color=yellow)$

The “show_id=true” option draws the vertex numbers, and vertex_size and vertex_color draws the vertices with the corresponding size and colour.

Note that a graph without any duplicate edges and without any loops is called a simple graph. And, an edge from a vertex U to V is not duplicate of an edge from vertex V to U in a directed graph but is duplicate in an undirected graph. Maxima’s package supports only simple graphs, i.e. graphs without duplicate edges and loops.

A simple graph can also be equivalently represented by adjacency of vertices, meaning by lists of adjacent vertices for every vertex. print_graph(<graph>) exactly displays those lists. The following code demonstration, in continuation from the previous code, demonstrates that:

(%i6) print_graph(g)$

Graph on 2 vertices with 1 edges.
Adjacencies:
  1 :  0
  0 :  1
(%i7) print_graph(dg)$

Digraph on 2 vertices with 1 arcs.
Adjacencies:
  1 :
  0 :  1
(%i8) quit();

create_graph(<num_of_vertices>, <edge_list>[, directed]) is another way of creating a graph using create_graph(). Here, the vertices are created as 0, 1, …, <num_of_vertices> – 1. So, both the above graphs could equivalently be created as follows:

$ maxima -q
(%i1) load(graphs)$
...
0 errors, 0 warnings
(%i2) g: create_graph(2, [[0, 1]]);
(%o2)                     GRAPH(2 vertices, 1 edges)
(%i3) dg: create_graph(2, [[0, 1]], directed=true);
(%o3)                     DIGRAPH(2 vertices, 1 arcs)
(%i4) quit();

make_graph(<vertices>, <predicate>[, directed]) is another interesting way of creating a graph, based on vertex connectivity conditions specified by the <predicate> function. <vertices> could be an integer or a set/list of vertices. If it is a positive integer, then the vertices would be 1, 2, …, <vertices>. In any case, <predicate> should be a function taking two arguments of the vertex type and returning true or false. make_graph() creates a graph with the vertices specified as above, and with the edges between the vertices for which the <predicate> function returns true.

A trivial case would be, if the <predicate> always returns true, then it would create a complete graph, i.e. a simple graph where all vertices are connected to each other. Here are a couple of demonstrations of make_graph():

$ maxima -q
(%i1) load(graphs)$
...
0 errors, 0 warnings
(%i2) f(i, j) := true$
(%i3) g: make_graph(6, f);
(%o3)                   GRAPH(6 vertices, 15 edges)
(%i4) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i5) f(i, j) := is(mod(i, j)=0)$
(%i6) g: make_graph(10, f, directed = true);
(%o6)                  DIGRAPH(10 vertices, 17 arcs)
(%i7) draw_graph(g, show_id=true, vertex_color=yellow, vertex_size=4)$
(%i8) quit();

Figures 26 shows the output graphs from the above code.

Figure 26: More simple graphs

Figure 26: More simple graphs

Graph varieties

One aware of graphs, would have known or at least heard of a variety of them. Here’s a list of some of them, available in Maxima’s graphs package (through functions):

  • Empty graph (empty_graph(n)) – Graph with a given n vertices but no edges
  • Complete graph (complete_graph(n)) – Simple graph with all possible edges for a given n vertices
  • Complete bipartite graph (complete_bipartite_graph(m, n)) – Simple graph with two set of vertices, having all possible edges between the vertices from the two sets, but with no edge between the vertices of the same set.
  • Cube graph (cube_graph(n)) – Graph representing an n-dimensional cube
  • Dodecahedron graph (dodecahedron_graph()) – Graph forming a 3-D polyhedron with 12 pentagonal faces
  • Cuboctahedron graph (cuboctahedron_graph()) – Graph forming a 3-D polyhedron with 8 triangular faces and 12 square faces
  • Icosahedron graph (icosahedron_graph()) – Graph forming a 3-D polyhedron with 20 triangular faces
  • Icosidodecahedron graph (icosidodecahedron_graph()) – Graph forming a 3-D uniform star polyhedron with 12 star faces and 20 triangular faces

And here follows a demonstration of the above, along with the visuals (left to right, top to bottom) in Figure 27:

$ maxima -q
(%i1) load(graphs)$
...
0 errors, 0 warnings
(%i2) g: empty_graph(5);
(%o2)                     GRAPH(5 vertices, 0 edges)
(%i3) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i4) g: complete_graph(5);
(%o4)                     GRAPH(5 vertices, 10 edges)
(%i5) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i6) g: complete_bipartite_graph(5, 3);
(%o6)                     GRAPH(8 vertices, 15 edges)
(%i7) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i8) g: cube_graph(3);
(%o8)                    GRAPH(8 vertices, 12 edges)
(%i9) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i10) g: cube_graph(4);
(%o10)                   GRAPH(16 vertices, 32 edges)
(%i11) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i12) g: dodecahedron_graph();
(%o12)                   GRAPH(20 vertices, 30 edges)
(%i13) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i14) g: cuboctahedron_graph();
(%o14)                   GRAPH(12 vertices, 24 edges)
(%i15) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i16) g: icosahedron_graph();
(%o16)                   GRAPH(12 vertices, 30 edges)
(%i17) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i18) g: icosidodecahedron_graph();
(%o18)                   GRAPH(30 vertices, 60 edges)
(%i19) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i20) quit();
Figure 27: Graph varieties

Figure 27: Graph varieties

Graph beauties

Graphs are really beautiful to visualize. Some of the many beautiful graphs, available in Maxima’s graphs package (through functions), are listed below:

  • Circulant graph (circulant_graph(n, [x, y, …])) – Graph with vertices 0, …, n-1, where every vertex is adjacent to its xth, yth, … vertices. Visually, it has a cyclic group of symmetries./li>
  • Flower graph (flower_snark(n)) – Graph like a flower with n petals and 4*n vertices
  • Wheel graph (wheel_graph(n)) – Graph like a wheel with n vertices
  • Clebsch graph (clebsch_graph()) – An another symmetrical graph beauty, named by J J Seidel
  • Frucht graph (frucht_graph()) – A graph with 12 vertices, 18 edges, and no nontrivial symmetries, such that every vertex have 3 neighbours. It is named after Robert Frucht
  • Grötzsch graph (grotzch_graph()) – A triangle-free graph with 11 vertices and 20 edges, named after Herbert Grötzsch
  • Heawood graph (heawood_graph()) – A symmetrical graph with 14 vertices and 21 edges, named after Percy John Heawood
  • Petersen graph (petersen_graph()) – A symmetrical graph with 10 vertices and 15 edges, named after Julius Petersen
  • Tutte graph (tutte_graph()) – A graph with 46 vertices and 69 edges, such that every vertex have 3 neighbours. It is named after W T Tutte

And here follows a demonstration of some of the above, along with the visuals (left to right, top to bottom) in Figure 28:

$ maxima -q
(%i1) load(graphs)$
...
0 errors, 0 warnings
(%i2) g: circulant_graph(10, [1, 3]);
(%o2)                   GRAPH(10 vertices, 20 edges)
(%i3) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i4) g: circulant_graph(10, [1, 3, 4, 6]);
(%o4)                   GRAPH(10 vertices, 40 edges)
(%i5) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i6) g: flower_snark(3);
(%o6)                   GRAPH(12 vertices, 18 edges)
(%i7) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i8) g: flower_snark(5);
(%o8)                   GRAPH(20 vertices, 30 edges)
(%i9) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i10) g: flower_snark(8);
(%o10)                   GRAPH(32 vertices, 48 edges)
(%i11) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i12) g: flower_snark(10);
(%o12)                   GRAPH(40 vertices, 60 edges)
(%i13) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i14) g: wheel_graph(3);
(%o14)                    GRAPH(4 vertices, 6 edges)
(%i15) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i16) g: wheel_graph(4);
(%o16)                    GRAPH(5 vertices, 8 edges)
(%i17) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i18) g: wheel_graph(5);
(%o18)                    GRAPH(6 vertices, 10 edges)
(%i19) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i20) g: wheel_graph(10);
(%o20)                   GRAPH(11 vertices, 20 edges)
(%i21) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i22) g: wheel_graph(100);
(%o22)                  GRAPH(101 vertices, 200 edges)
(%i23) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i24) g: clebsch_graph();
(%o24)                   GRAPH(16 vertices, 40 edges)
(%i25) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i26) g: grotzch_graph();
(%o26)                   GRAPH(11 vertices, 20 edges)
(%i27) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i28) g: petersen_graph();
(%o28)                   GRAPH(10 vertices, 15 edges)
(%i29) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i30) g: tutte_graph();
(%o30)                   GRAPH(46 vertices, 69 edges)
(%i31) draw_graph(g, show_id=true, vertex_color=yellow)$
(%i32) quit();
Figure 28: Graph beauties

Figure 28: Graph beauties

What next?

With all these visualizations, you may be wondering, what to do with these apart from just staring at them. Visualizations were just to motivate you, visually. In fact, every graph has a particular set of properties, which distinguishes it from the other. And, there is a lot of beautiful mathematics involved with these properties. If you are motivated enough, watch out for playing around with these graphs and their properties.

Twenty-fourth Article >>

   Send article as PDF   

Manage your Finances with Maxima

2 Replies

This twenty-second article of the mathematical journey through open source, show cases managing our usual financial computations using the financial package of Maxima.

<< Twenty-first Article

In today’s credit-oriented world, loan, EMIs, the principal, interest rates, savings, etc are the common lingo. How well do we really understand these, or for that matter are able to compute the actual finances related to such things? Or, do we just accept what the other party has to provide us. One may say, either get into the gory details to understand and verify, or forget about understanding and just accept it if the offering suits you – you cannot not understand and also verify at the same time. However, with the finance package of Maxima, you can actually verify without getting much into the computational details. This article is going to exactly walk you through that.

Basic operations

To be able to use the finance package of Maxima, the first thing to do, after getting the Maxima prompt is to load the finance package, using load(). And, then go ahead to do the various computations. days360() is the simplest function to give the number of days between two dates, assuming 360 days per year, 30 days per month – one of the common interest calculation norm.

$ maxima -q
(%i1) load(finance);
(%o1)     /usr/share/maxima/5.24.0/share/contrib/finance/finance.mac
(%i2) days360(2014, 1, 1, 2014, 10, 1);
(%o2)                                 270
(%i3) days360(2014, 1, 1, 2014, 12, 31);
(%o3)                                 360
(%i4) days360(2014, 1, 1, 2014, 3, 1);
(%o4)                                 60
(%i5) days360(2014, 1, 1, 2015, 1, 1);
(%o5)                                 360
(%i6) quit();

Note that days360() takes from and to dates, each as a triplet of year, month, date.

One common computation, we often deal with, is the final amount we would get after applying a particular rate of compound interest on a particular principal amount – just use fv(<rate>, <principal>, <period>) for computing the future value of the <principal>, at the compound interest rate of <rate>, after <period> periods of the rate. As an example, what would be the future value of investing ₹10,000 for 10 years at the compound interest rate of 15% per year – just call fv(0.15, 10000, 10);

$ maxima -q
(%i1) load(finance)$ /* $ to suppress its output */
(%i2) fv(0.15, 10000, 10);
(%o2)                          40455.57735707907

If you are interested in, how exactly did it calculate, or what is the formula, that is where Maxima is fun to play with symbols. Just use some symbols, instead of actual numbers:

(%i3) string(fv(r, p, n));
(%o3)                              p*(r+1)^n
(%i4) quit();

And, there you go. p*(r+1)^n is the future value for investing p amount for n periods at the compounded interest rate of r per period.

How about doing an inverse computation? Suppose, I want to get ₹10,00,000 after 5 years from my investment today at the interest rate of 10.75%. Now for that, how much should I invest today? Don’t scratch your head, just call pv(<rate>, <future_val>, <period>);

(%i1) load(finance)$
(%i2) pv(0.1075, 1000000, 5);
(%o2)                          600179.7323625274
(%i3) fv(0.1075, 600180, 5);
(%o3)                          1000000.445928875
(%i4) quit();

That requires an investment of ₹6,00,180. Yes, so if you invest that amount, the future value would be ₹10,00,000 – that’s the check done above using fv().

Loans and EMIs

Fundamental thought behind today’s culture of buying a house, a car, or even smaller items is “let’s buy it on loan and pay it off in EMIs (equated monthly installments)”. These terms would be familiar to most of you, so no need to explain them. But, how do you compute these? One would say, there are some complicated formulas to do so. Yes, you are right. But, you don’t need to worry about any of them. Just tell Maxima to give you the complete schedule for your loan, at a given rate of interest, for a given period of time, using amortization(). The first example below provides the schedule for a home loan of ₹20 lakhs at an interest rate of 9.25% per annum (p.a.) for 5 years. The various columns in the output schedule provide the following information:

  • n: installment payment time – year for our example
  • Balance: principal + interest left over to be paid out, after the current installment is paid out
  • Interest: interest part being paid out in the current installment
  • Amortization: principal part being paid out in the current installment
  • Payment: current installment to be paid out, i.e. the EYI (equated yearly installment)

What is this EYI? We were supposed to talk only about EMI, right? Okay, in that case, we need to convert the rate of interest and the period, both in terms of months. So, we need to divide the per annum rate of interest by 12 to get it per month, and multiply the number of years by 12 to get that in number of months. That is exactly what the second example below shows. Note the 60 EMIs to be paid out over the period of 5 years.

$ maxima -q
(%i1) load(finance)$
(%i2) amortization(0.0925, 2000000, 5)$
           "n"      "Balance"    "Interest"  "Amortization"  "Payment"      
          0.000   2000000.000         0.000         0.000         0.000  
          1.000   1667475.420    185000.000    332524.580    517524.580  
          2.000   1304192.317    154241.476    363283.103    517524.580  
          3.000    907305.527    120637.789    396886.790    517524.580  
          4.000    473706.709     83925.761    433598.818    517524.580  
          5.000         0.000     43817.871    473706.709    517524.580  

(%i3) amortization(0.0925/12, 2000000, 5*12)$
           "n"      "Balance"    "Interest"  "Amortization"  "Payment"      
          0.000   2000000.000         0.000         0.000         0.000  
          1.000   1973656.870     15416.667     26343.130     41759.797  
          2.000   1947110.679     15213.605     26546.192     41759.797  
          3.000   1920359.860     15008.978     26750.818     41759.797  
          4.000   1893402.837     14802.774     26957.023     41759.797  
          5.000   1866238.021     14594.980     27164.816     41759.797  
          6.000   1838863.809     14385.585     27374.212     41759.797  
          7.000   1811278.588     14174.575     27585.221     41759.797  
          8.000   1783480.730     13961.939     27797.857     41759.797  
          9.000   1755468.598     13747.664     28012.133     41759.797  
         10.000   1727240.538     13531.737     28228.059     41759.797  
         11.000   1698794.887     13314.146     28445.651     41759.797  
         12.000   1670129.968     13094.877     28664.919     41759.797  
         13.000   1641244.090     12873.919     28885.878     41759.797  
         14.000   1612135.550     12651.257     29108.540     41759.797  
         15.000   1582802.632     12426.878     29332.918     41759.797  
         16.000   1553243.605     12200.770     29559.026     41759.797  
         17.000   1523456.728     11972.919     29786.877     41759.797  
         18.000   1493440.244     11743.312     30016.484     41759.797  
         19.000   1463192.382     11511.935     30247.861     41759.797  
         20.000   1432711.360     11278.775     30481.022     41759.797  
         21.000   1401995.381     11043.817     30715.980     41759.797  
         22.000   1371042.632     10807.048     30952.749     41759.797  
         23.000   1339851.289     10568.454     31191.343     41759.797  
         24.000   1308419.513     10328.020     31431.776     41759.797  
         25.000   1276745.450     10085.734     31674.063     41759.797  
         26.000   1244827.233      9841.580     31918.217     41759.797  
         27.000   1212662.979      9595.543     32164.253     41759.797  
         28.000   1180250.793      9347.610     32412.186     41759.797  
         29.000   1147588.763      9097.767     32662.030     41759.797  
         30.000   1114674.963      8845.997     32913.800     41759.797  
         31.000   1081507.453      8592.286     33167.510     41759.797  
         32.000   1048084.276      8336.620     33423.177     41759.797  
         33.000   1014403.463      8078.983     33680.814     41759.797  
         34.000    980463.026      7819.360     33940.437     41759.797  
         35.000    946260.965      7557.736     34202.061     41759.797  
         36.000    911795.264      7294.095     34465.702     41759.797  
         37.000    877063.889      7028.422     34731.375     41759.797  
         38.000    842064.793      6760.701     34999.096     41759.797  
         39.000    806795.913      6490.916     35268.880     41759.797  
         40.000    771255.168      6219.052     35540.745     41759.797  
         41.000    735440.463      5945.092     35814.705     41759.797  
         42.000    699349.687      5669.020     36090.776     41759.797  
         43.000    662980.711      5390.821     36368.976     41759.797  
         44.000    626331.390      5110.476     36649.320     41759.797  
         45.000    589399.565      4827.971     36931.825     41759.797  
         46.000    552183.057      4543.288     37216.508     41759.797  
         47.000    514679.671      4256.411     37503.386     41759.797  
         48.000    476887.197      3967.322     37792.474     41759.797  
         49.000    438803.406      3676.005     38083.791     41759.797  
         50.000    400426.052      3382.443     38377.354     41759.797  
         51.000    361752.873      3086.617     38673.179     41759.797  
         52.000    322781.588      2788.512     38971.285     41759.797  
         53.000    283509.900      2488.108     39271.689     41759.797  
         54.000    243935.492      2185.389     39574.408     41759.797  
         55.000    204056.031      1880.336     39879.461     41759.797  
         56.000    163869.167      1572.932     40186.865     41759.797  
         57.000    123372.528      1263.158     40496.638     41759.797  
         58.000     82563.728       950.997     40808.800     41759.797  
         59.000     41440.360       636.429     41123.368     41759.797  
         60.000         0.000       319.436     41440.360     41759.797  

(%i4) quit();

If you want to be a little stylish in your monthly payments, that is, instead of equal payments, you want to do increasing payments, Maxima could help you with that as well. arit_amortization() and geo_amortization() are two such functions, which provides the schedule for increasing payments with fixed amount increments and fixed rate increments, respectively. Here’s a small demo of the same:

$ maxima -q
(%i1) load(finance)$
(%i2) amortization(0.10, 100, 5)$
           "n"       "Balance"     "Interest"   "Amortization"  "Payment"      
          0.000       100.000         0.000         0.000         0.000  
          1.000        83.620        10.000        16.380        26.380  
          2.000        65.603         8.362        18.018        26.380  
          3.000        45.783         6.560        19.819        26.380  
          4.000        23.982         4.578        21.801        26.380  
          5.000         0.000         2.398        23.982        26.380  

(%i3) arit_amortization(0.10, 10, 100, 5)$
           "n"       "Balance"     "Interest"   "Amortization"  "Payment"      
          0.000       100.000         0.000         0.000         0.000  
          1.000       101.722        10.000        -1.722         8.278  
          2.000        93.615        10.172         8.106        18.278  
          3.000        74.698         9.362        18.917        28.278  
          4.000        43.890         7.470        30.809        38.278  
          5.000         0.000         4.389        43.890        48.278  

(%i4) geo_amortization(0.10, 0.05, 100, 5)$
           "n"       "Balance"     "Interest"   "Amortization"  "Payment"      
          0.000       100.000         0.000         0.000         0.000  
          1.000        85.907        10.000        14.093        24.093  
          2.000        69.200         8.591        16.707        25.298  
          3.000        49.558         6.920        19.642        26.562  
          4.000        26.623         4.956        22.935        27.891  
          5.000         0.000         2.662        26.623        29.285  

(%i5) quit();

%i2 has been provided for comparative analysis. %i3 shows the incremental payout with increments of ₹10 (the second parameter of arit_amortization()). %i4 shows the incremental payout with increments at the rate of 5% (the second parameter of geo_amortization()). Both these computations could be done in decrements as well – just pass the second parameter negative.

Plan your savings

Say, you have a savings account like PPF (public provident fund), giving you interest at a rate of 8% p.a., and at the end of 5 years, you want to have save ₹1,00,000. So, what should be your minimum yearly deposit into your account. It is not just divide by 5, as interest would be also added to your savings. saving() shows you the complete schedule as follows:

$ maxima -q
(%i1) load(finance)$
(%i2) saving(0.08 /* interest rate */, 100000 /* final savings */, 5 /* years */)$
           "n"       "Balance"     "Interest"    "Payment"      
          0.000          0.000         0.000         0.000  
          1.000      17045.645         0.000     17045.645  
          2.000      35454.943      1363.652     17045.645  
          3.000      55336.983      2836.395     17045.645  
          4.000      76809.588      4426.959     17045.645  
          5.000     100000.000      6144.767     17045.645

And, the minimum yearly deposit to be done is ₹17,046. The “Balance” and “Interest” columns, respectively, tell you about the balance and the interest accumulated in the corresponding year. If you are only interested in knowing the minimum amount to be deposited, you may just use the annuity_fv() function – basically computing the annuity of ₹17,046 every year for 5 years to have a future saving of ₹1,00,000 after 5 years.

(%i3) annuity_fv(0.08, 100000, 5);
(%o3)                         17045.64545668365
(%i4) quit();

Project planning

Finance management is a key ingredient of any project planning, whether it be a professional or personal project. Assume a project would take n years, with given yearly expenses, and say the available interest rate is r p.a. Then, a common question, every project manager needs to answer is what is the net present value (NPV) of the project, which needs to be invested for the project. The answer to this basic question is more often than not, one of the important factors for deciding whether to take up this project or not. Maxima provides npv() to compute the same. As an example, if a project needs ₹100, ₹200, ₹150, ₹600 over 4 years, respectively, and the current interest rate is 7%, what is the NPV? It would be ₹848 as shown below:

$ maxima -q
(%i1) load(finance)$
(%i2) npv(0.07, [100, 200, 150, 600]);
(%o2)                         848.3274983420189
(%i3) quit();

Another common practice to select between various projects is to compute the benefit to cost ratio of the various projects, and select the one with best ratio. The benefit to cost ratio for a given interest rate r could be computed using benefit_cost(). Here’s an example to demonstrate the same, assuming 18% as the rate of interest:

  • Project P1 (2 years): Yearly Benefits (100, 200), Yearly Costs (150, 50)
  • Project P2 (3 years): Yearly Benefits (0, 200, 100), Yearly Costs (100, 100, 0)
  • Project P3 (4 years): Yearly Benefits (0, 200, 200, 100), Yearly Costs (100, 100, 50, 50)
$ maxima -q
(%i1) load(finance)$
(%i2) benefit_cost(0.18, [100, 200], [150, 50]);
(%o2)                         1.400881057268722
(%i3) benefit_cost(0.18, [0, 200, 100], [100, 100, 0]);
(%o3)                         1.306173223448919
(%i4) benefit_cost(0.18, [0, 200, 200, 100], [100, 100, 50, 50]);
(%o4)                         1.489492494361802
(%i5) quit();

And clearly, over the long run the 4-year project P3 has better benefit cost ratio. But if only looking for shorter term benefits, then one may even go for project P1 as well.

Twenty-third Article >>

   Send article as PDF   

Laplace Transforms: Solving Engineering Problems

2 Replies

This twenty-first article of the mathematical journey through open source, demonstrates the Laplace transforms through Maxima.

<< Twentieth Article

In higher mathematics, transforms play an important role. A transform is mathematical logic to transform or convert a mathematical expression into another mathematical expression, typically from one domain to another. Laplace and Fourier are two very common examples, transforming from time domain to frequency domain. In general, such transforms have their corresponding inverse transforms. And this combination of direct and inverse transforms is very powerful in solving many real life engineering problems. Focus of this article is Laplace and its inverse transform, along with some problem solving highlights.

Laplace transform

Mathematically, Laplace transform F(s) of a function f(t) is defined as follows:

F(s) = \int_0^\infty{f(t) * exp(-s*t) dt},

where t represents time and s represents complex angular frequency.

To demonstrate it, let’s take a simple example of f(t) = 1. Substituting and integrating, we get F(s) = 1/s. Maxima has the function laplace() to do the same. In fact, with that, we can choose to let our variables t and s be anything else as well. But, as per our mathematical notations, preserving them as t and s would be the most appropriate. Let’s start with some basic Laplace transforms:
(Note that string() has been used to just flatten the expression)

$ maxima -q
(%i1) string(laplace(1, t, s));
(%o1)                                 1/s
(%i2) string(laplace(t, t, s));
(%o2)                                1/s^2
(%i3) string(laplace(t^2, t, s));
(%o3)                                2/s^3
(%i4) string(laplace(t+1, t, s));
(%o4)                              1/s+1/s^2
(%i5) string(laplace(t^n, t, s));
Is  n + 1  positive, negative, or zero?

p; /* Our input */
(%o5)                         gamma(n+1)*s^(-n-1)
(%i6) string(laplace(t^n, t, s));
Is  n + 1  positive, negative, or zero?

n; /* Our input */
(%o6)                  gamma_incomplete(n+1,0)*s^(-n-1)
(%i7) string(laplace(t^n, t, s));
Is  n + 1  positive, negative, or zero?

z; /* Our input, making it non-solvable */
(%o7)                          'laplace(t^n,t,s)
(%i8) string(laplace(1/t, t, s)); /* Non-solvable */
(%o8)                          'laplace(1/t,t,s)
(%i9) string(laplace(1/t^2, t, s)); /* Non-solvable */
(%o9)                         'laplace(1/t^2,t,s)
(%i10) quit();

Note that, in the above examples, the expression is preserved as is, in case of non-solvability.

laplace() is designed to understand various symbolic functions, such as sin(), cos(), sinh(), cosh(), log(), exp(), delta(), erf(). delta() is Dirac delta function, and erf() is the error function – others being the usual mathematical functions.

$ maxima -q
(%i1) string(laplace(sin(t), t, s));
(%o1)                              1/(s^2+1)
(%i2) string(laplace(sin(w*t), t, s));
(%o2)                             w/(w^2+s^2)
(%i3) string(laplace(cos(t), t, s));
(%o3)                              s/(s^2+1)
(%i4) string(laplace(cos(w*t), t, s));
(%o4)                             s/(w^2+s^2)
(%i5) string(laplace(sinh(t), t, s));
(%o5)                              1/(s^2-1)
(%i6) string(laplace(sinh(w*t), t, s));
(%o6)                            -w/(w^2-s^2)
(%i7) string(laplace(cosh(t), t, s));
(%o7)                              s/(s^2-1)
(%i8) string(laplace(cosh(w*t), t, s));
(%o8)                            -s/(w^2-s^2)
(%i9) string(laplace(log(t), t, s));
(%o9)                         (-log(s)-%gamma)/s
(%i10) string(laplace(exp(t), t, s));
(%o10)                              1/(s-1)
(%i11) string(laplace(delta(t), t, s));
(%o11)                                 1
(%i12) string(laplace(erf(t), t, s));
(%o12)                     %e^(s^2/4)*(1-erf(s/2))/s
(%i13) quit();

Interpreting the transform

A Laplace transform is typically a fractional expression consisting of a numerator and a denominator. Solving the denominator by equating it to zero, gives the various complex frequencies associated with the original function. These are called the poles of the function. For example, the Laplace transform of sin(w * t) is w/(s^2 + w^2), where the denominator is s^2 + w^2. Equating that to zero and solving it, gives complex frequency s = +iw, -iw; thus indicating that the frequency of the original expression sin(w * t) is w, which indeed is w. Here goes few demonstrations for the same:

$ maxima -q
(%i1) string(laplace(sin(w*t), t, s));
(%o1)                             w/(w^2+s^2)
(%i2) string(denom(laplace(sin(w*t), t, s))); /* The Denominator */
(%o2)                               w^2+s^2
(%i3) string(solve(denom(laplace(sin(w*t), t, s)), s)); /* The Poles */
(%o3)                        [s = -%i*w,s = %i*w]
(%i4) string(solve(denom(laplace(sinh(w*t), t, s)), s));
(%o4)                           [s = -w,s = w]
(%i5) string(solve(denom(laplace(cos(w*t), t, s)), s));
(%o5)                        [s = -%i*w,s = %i*w]
(%i6) string(solve(denom(laplace(cosh(w*t), t, s)), s));
(%o6)                           [s = -w,s = w]
(%i7) string(solve(denom(laplace(exp(w*t), t, s)), s));
(%o7)                               [s = w]
(%i8) string(solve(denom(laplace(log(w*t), t, s)), s));
(%o8)                               [s = 0]
(%i9) string(solve(denom(laplace(delta(w*t), t, s)), s));
(%o9)                                 []
(%i10) string(solve(denom(laplace(erf(w*t), t, s)), s));
(%o10)                              [s = 0]
(%i11) quit();

Involved Laplace transforms

laplace() also understands derivative() / diff(), integrate(), sum(), and ilt() – the inverse Laplace transform. Here are some interesting transforms showing the same:

$ maxima -q
(%i1) laplace(f(t), t, s);
(%o1)                         laplace(f(t), t, s)
(%i2) string(laplace(derivative(f(t), t), t, s));
(%o2)                      s*'laplace(f(t),t,s)-f(0)
(%i3) string(laplace(integrate(f(x), x, 0, t), t, s));
(%o3)                        'laplace(f(t),t,s)/s
(%i4) string(laplace(derivative(sin(t), t), t, s));
(%o4)                              s/(s^2+1)
(%i5) string(laplace(integrate(sin(t), t), t, s));
(%o5)                             -s/(s^2+1)
(%i6) string(sum(t^i, i, 0, 5));
(%o6)                         t^5+t^4+t^3+t^2+t+1
(%i7) string(laplace(sum(t^i, i, 0, 5), t, s));
(%o7)                1/s+1/s^2+2/s^3+6/s^4+24/s^5+120/s^6
(%i8) string(laplace(ilt(1/s, s, t), t, s));
(%o8)                                 1/s
(%i9) quit();

Note the usage of ilt() – inverse Laplace transform in the %i8 of above example. Calling laplace() and ilt() one after the other cancels their effect – that is what is meant by inverse. Let’s look into some common inverse Laplace transforms.

Inverse Laplace transforms

$ maxima -q
(%i1) string(ilt(1/s, s, t));
(%o1)                                  1
(%i2) string(ilt(1/s^2, s, t));
(%o2)                                  t
(%i3) string(ilt(1/s^3, s, t));
(%o3)                                t^2/2
(%i4) string(ilt(1/s^4, s, t));
(%o4)                                t^3/6
(%i5) string(ilt(1/s^5, s, t));
(%o5)                               t^4/24
(%i6) string(ilt(1/s^10, s, t));
(%o6)                             t^9/362880
(%i7) string(ilt(1/s^100, s, t));
(%o7) t^99/9332621544394415268169923885626670049071596826438162146859296389521759999
	3229915608941463976156518286253697920827223758251185210916864000000000000000
	0000000
(%i8) string(ilt(1/(s-a), s, t));
(%o8)                              %e^(a*t)
(%i9) string(ilt(1/(s^2-a^2), s, t));
(%o9)                   %e^(a*t)/(2*a)-%e^-(a*t)/(2*a)
(%i10) string(ilt(s/(s^2-a^2), s, t));
(%o10)                      %e^(a*t)/2+%e^-(a*t)/2
(%i11) string(ilt(1/(s^2+a^2), s, t));
Is  a  zero or nonzero?

n; /* Our input */
(%o11)                            sin(a*t)/a
(%i12) string(ilt(s/(s^2+a^2), s, t));
Is  a  zero or nonzero?

n; /* Our input */
(%o12)                             cos(a*t)
(%i13) assume(a < 0) or assume(a > 0)$
(%i14) string(ilt(1/(s^2+a^2), s, t));
(%o14)                             sin(a*t)/a
(%i15) string(ilt(s/(s^2+a^2), s, t));
(%o15)                              cos(a*t)
(%i16) string(ilt((s^2+s+1)/(s^3+s^2+s+1), s, t));
(%o16)                      sin(t)/2+cos(t)/2+%e^-t/2
(%i17) string(laplace(sin(t)/2+cos(t)/2+%e^-t/2, t, s));
(%o17)               s/(2*(s^2+1))+1/(2*(s^2+1))+1/(2*(s+1))
(%i18) string(rat(laplace(sin(t)/2+cos(t)/2+%e^-t/2, t, s)));
(%o18)                       (s^2+s+1)/(s^3+s^2+s+1)
(%i19) quit();

Observe that if we take the Laplace transform of the above %o outputs, they would give back the expressions, which are input to ilt() of the corresponding %i’s. %i18 specifically shows one such example. It does laplace() of the output at %o16, giving back the expression, which was input to ilt() of %i16.

Solving differential and integral equations

Now, with these insights, we can easily solve many interesting and otherwise complex problems. One of them is solving differential equations. Let’s take a simple example of solving f'(t) + f(t) = e^t, where f(0) = 0. First we take the Laplace transform of the equation. Then substitute the value for f(0), and simplify to obtain the Laplace of f(t), i.e. F(s). And, then finally compute the inverse Laplace transform of F(s), to get the solution for f(t).

$ maxima -q
(%i1) string(laplace(diff(f(t), t) + f(t) = exp(t), t, s));
(%o1)      s*'laplace(f(t),t,s)+'laplace(f(t),t,s)-f(0) = 1/(s-1)

Substituting f(0) as 0, and simplifying we get, laplace(f(t),t,s) = 1/((s-1)*(s+1)), for which we do an inverse Laplace transform:

(%i2) string(ilt(1/((s-1)*(s+1)), s, t));
(%o2)                          %e^t/2-%e^-t/2
(%i3) quit();

Thus, giving f(t) = (e^t – e^-t) / 2, i.e. sinh(t), which definitely satisfies the given differential equation.

Similarly, we can solve equations with integrals. And, why only integrals – rather equations with both differentials and integrals. Such equations come very often in solving for electrical circuits with resistors, capacitors, and inductors. Let’s again take a simple example to demonstrate the fact. Say, we have 1 ohm resistor, 1 farad capacitor, and 1 henry inductor in series being powered by a sinusoidal voltage source of frequency w. What would be the current in the circuit, assuming it to be zero at t = 0? It would yield the following equation: R * i(t) + 1/C * ∫ i(t) dt + L * di(t)/dt = sin(w*t), where R = 1, C = 1, L =1.

So, the equation simplifies to i(t) + ∫ i(t) dt + di(t)/dt = sin(w*t). Now, following the procedure as described above, we do the following steps:

$ maxima -q
(%i1) string(laplace(i(t) + integrate(i(x), x, 0, t) + diff(i(t), t) = sin(w*t),
	t, s));
(%o1) s*'laplace(i(t),t,s)+'laplace(i(t),t,s)/s+'laplace(i(t),t,s)-i(0) = w/(w^2+s^2)

Substituting i(0) as 0, and simplifying, we get laplace(i(t), t, s) = w/((w^2+s^2)*(s+1/s+1)). Solving that by inverse Laplace transform, we very easily get the complex expression for i(t) as follows:

(%i2) string(ilt(w/((w^2+s^2)*(s+1/s+1)), s, t));
Is  w  zero or nonzero?

n; /* Our input: Non-zero frequency */
(%o2) w^2*sin(t*w)/(w^4-w^2+1)-(w^3-w)*cos(t*w)/(w^4-w^2+1)+%e^-(t/2)*
	(sin(sqrt(3)*t/2)*(-(w^3-w)/(w^4-w^2+1)-2*w/(w^4-w^2+1))/sqrt(3)+
	cos(sqrt(3)*t/2)*(w^3-w)/(w^4-w^2+1))
(%i3) quit();

Twenty-second Article >>

   Send article as PDF   

Lists: The Building Blocks of Maxima

7 Replies

This twentieth article of the mathematical journey through open source, showcases the list manipulations in Maxima.

<< Nineteenth Article

Lists are the basic building blocks of Maxima. The fundamental reason being that Maxima is implemented in Lisp, the building blocks of which, are also lists.

Creation of lists

To begin with, let us walk through the ways to create a list. Simplest way to have a list in Maxima is to just define it using [ ]. So, [x, 5, 3, 2*y] is a list consisting of 4 members. However, Maxima provides two powerful functions for automatically generating lists: makelist(), create_list().

makelist() can take two forms. makelist(e, x, x0, xn) creates and returns a list using the expression e, evaluated for x using the values ranging from x0 to xn. makelist(e, x, L) – creates and returns a list using the expression e, evaluated for x using the members of the list L. Check out the example below for better clarity.

$ maxima -q
(%i1) makelist(2 * i, i, 1, 5);
(%o1)                        [2, 4, 6, 8, 10]
(%i2) makelist(concat(x, 2 * i - 1), i, 1, 5);
(%o2)                        [x1, x3, x5, x7, x9]
(%i3) makelist(concat(x, 2), x, [a, b, c, d]);
(%o3)                          [a2, b2, c2, d2]
(%i4) quit();

Note the interesting usage of concat() to just concatenate its arguments. Note that, makelist() is limited by the variation it can have – to be specific, just one – i in the first two examples, and x in the last one. If we want more, that is where the create_list() function comes into play.

create_list(f, x1, L1, …, xn, Ln) – creates and returns a list with members of the form f, evaluated for the variables x1, …, xn using the values from the corresponding lists L1, …, Ln. Here is a glimpse of its power:

$ maxima -q
(%i1) create_list(concat(x, y), x, [p, q], y, [1, 2]);
(%o1)                          [p1, p2, q1, q2]
(%i2) create_list(concat(x, y, z), x, [p, q], y, [1, 2], z, [a, b]);
(%o2)              [p1a, p1b, p2a, p2b, q1a, q1b, q2a, q2b]
(%i3) create_list(concat(x, y, z), x, [p, q], y, [1, 2, 3], z, [a, b]);
(%o3)    [p1a, p1b, p2a, p2b, p3a, p3b, q1a, q1b, q2a, q2b, q3a, q3b]
(%i4) quit();

Note the “all possible combinations” being created using the values for the variables x, y, z.

Once we have lists created, Maxima provides a host of functions to play around with them. Let’s take a walk through them.

Testing the lists

The following set of functions demonstrates the various checks on lists:

  • atom(v) – returns true if v is an atomic element, false otherwise
  • listp(L) – returns true if L is a list, false otherwise
  • member(v, L) – returns true if v is a member of list L, false otherwise
  • some(p, L) – returns true if predicate p is true for at least one member of list L, false otherwise
  • every(p, L) – returns true if predicate p is true for all members of list L, false otherwise
$ maxima -q
(%i1) atom(5);
(%o1)                                true
(%i2) atom([5]);
(%o2)                                false
(%i3) listp(x);
(%o3)                                false
(%i4) listp([x]);
(%o4)                                true
(%i5) listp([x, 5]);
(%o5)                                true
(%i6) member(x, [a, b, c]);
(%o6)                                false
(%i7) member(x, [a, x, c]);
(%o7)                                true
(%i8) some(primep, [1, 4, 9]);
(%o8)                                false
(%i9) some(primep, [1, 2, 4, 9]);
(%o9)                                true
(%i10) every(integerp, [1, 2, 4, 9]);
(%o10)                               true
(%i11) every(integerp, [1, 2, 4, x]);
(%o11)                               false
(%i12) quit();

List recreations

Next, is a set of functions operating on list(s) to create and return new lists:

  • cons(v, L) – returns a list with v followed by members of L
  • endcons(v, L) – returns a list with members of L followed by v
  • rest(L, n) – returns a list with members of L, except the first n members (if n is non-negative), otherwise except the last -n members. n is optional, in which case, it is taken as 1.
  • join(L1, L2) – returns a list with members of L1 and L2 interspersed
  • delete(v, L, n) – returns a list like L but with the first n occurences of v, deleted from it. n is optional, in which case all occurences of v are deleted
  • append(L1, …, Ln) – returns a list with members of L1, …, Ln, one after the other
  • unique(L) – returns a list obtained by removing the duplicate members in the list L
  • reverse(L) – returns a list with members of the list L in reverse order
$ maxima -q
(%i1) L: makelist(i, i, 1, 10);
(%o1)                   [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
(%i2) cons(0, L);
(%o2)                 [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
(%i3) endcons(11, L);
(%o3)                 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
(%i4) rest(L);
(%o4)                    [2, 3, 4, 5, 6, 7, 8, 9, 10]
(%i5) rest(L, 3);
(%o5)                       [4, 5, 6, 7, 8, 9, 10]
(%i6) rest(L, -3);
(%o6)                        [1, 2, 3, 4, 5, 6, 7]
(%i7) join(L, [a, b, c, d]);
(%o7)                      [1, a, 2, b, 3, c, 4, d]
(%i8) delete(6, L);
(%o8)                    [1, 2, 3, 4, 5, 7, 8, 9, 10]
(%i9) delete(4, delete(6, L));
(%o9)                      [1, 2, 3, 5, 7, 8, 9, 10]
(%i10) delete(4, delete(6, join(L, L)));
(%o10)        [1, 1, 2, 2, 3, 3, 5, 5, 7, 7, 8, 8, 9, 9, 10, 10]
(%i11) L1: rest(L, 7);
(%o11)                            [8, 9, 10]
(%i12) L2: rest(rest(L, -3), 3);
(%o12)                           [4, 5, 6, 7]
(%i13) L3: rest(L, -7);
(%o13)                             [1, 2, 3]
(%i14) append(L1, L2, L3);
(%o14)                  [8, 9, 10, 4, 5, 6, 7, 1, 2, 3]
(%i15) reverse(L);
(%o15)                   [10, 9, 8, 7, 6, 5, 4, 3, 2, 1]
(%i16) join(reverse(L), L); 
(%o16)   [10, 1, 9, 2, 8, 3, 7, 4, 6, 5, 5, 6, 4, 7, 3, 8, 2, 9, 1, 10]
(%i17) unique(join(reverse(L), L)); 
(%o17)                   [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
(%i18) L;
(%o18)                   [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
(%i19) quit();

Note that the list L is still not modified. For that matter, even L1, L2, L3 are not modified. In fact, that is what is meant by that all these functions recreate new modified lists, rather than modifying the existing ones.

List extractions

Here goes a set of functions extracting the various members of a list. first(L), second(L), third(L), fourth(L), fifth(L), sixth(L), seventh(L), eight(L), ninth(L), tenth(L) – respectively return the first, second, … member of the list L. last(L) – returns the last member of the list L

$ maxima -q
(%i1) L: create_list(i * x, x, [a, b, c], i, [1, 2, 3, 4]);
(%o1)       [a, 2 a, 3 a, 4 a, b, 2 b, 3 b, 4 b, c, 2 c, 3 c, 4 c]
(%i2) first(L);
(%o2)                                  a
(%i3) seventh(L);
(%o3)                                 3 b
(%i4) last(L);
(%o4)                                 4 c
(%i5) third(L); last(L);
(%o5)                                 3 a
(%o6)                                 4 c
(%i7) L;
(%o7)       [a, 2 a, 3 a, 4 a, b, 2 b, 3 b, 4 b, c, 2 c, 3 c, 4 c]
(%i8) quit();

Again, note that the list L is still not modified. But, why have we been talking of that? To bring out the fact, that we may need to modify the existing lists, and none of the above functions would do that. Note that, we may achieve that by assigning the return values of the various list recreation functions back to the original list. However, there are few functions, which does that right away.

List manipulations

The following are the two list manipulating functions provided by Maxima:

  • push(v, L) – inserts v at the beginning of the list L
  • pop(L) – removes and returns the first element from list L

Note that L must be a symbol bound to a list, not the list itself, in both the above functions, for them to modify it. Also, these functionalities are not available by default, so we need to load the basic Maxima file. Check out the demonstration below.

We may display L after doing these operations, or even check the length of L to verify the actual modification of L. And, in case we need to preserve a copy of the list, function copylist() can be used, as such.

$ maxima -q
(%i1) L: makelist(2 * x, x, 1, 10);
(%o1)                [2, 4, 6, 8, 10, 12, 14, 16, 18, 20]
(%i2) push(0, L); /* This doesn't work */
(%o2)            push(0, [2, 4, 6, 8, 10, 12, 14, 16, 18, 20])
(%i3) pop(L); /* Nor does this work */
(%o3)              pop([2, 4, 6, 8, 10, 12, 14, 16, 18, 20])
(%i4) load(basic); /* Loading the basic Maxima file */
(%o4)           /usr/share/maxima/5.24.0/share/macro/basic.mac
(%i5) push(0, L); /* Now, this works */
(%o5)               [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20]
(%i6) L;
(%o6)               [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20]
(%i7) pop(L); /* Even this works */
(%o7)                                  0
(%i8) L;
(%o8)                [2, 4, 6, 8, 10, 12, 14, 16, 18, 20]
(%i9) K: copylist(L);
(%o9)                [2, 4, 6, 8, 10, 12, 14, 16, 18, 20]
(%i10) length(L);
(%o10)                                10
(%i11) pop(L);
(%o11)                                 2
(%i12) length(L);
(%o12)                                 9
(%i13) K;
(%o13)               [2, 4, 6, 8, 10, 12, 14, 16, 18, 20]
(%i14) L;
(%o14)                 [4, 6, 8, 10, 12, 14, 16, 18, 20]
(%i15) pop([1, 2, 3]); /* Actual list is not allowed */
arg must be a symbol [1, 2, 3]
#0: symbolcheck(x=[1,2,3])(basic.mac line 22)
#1: pop(l=[1,2,3])(basic.mac line 26)
 -- an error. To debug this try: debugmode(true);
(%i16) quit();

Advanced list operations

And finally, if you are still with me, here is a bonus of two sophisticated list operations:

  • sublist_indices(L, p) – returns the list indices for the members of the list L, for which predicate p is true.
  • assoc(k, L, d) – L must have all its members in the form of x op y, where op is some binary operator. Then, assoc() searches for k in the left operand of the members of L. If found, it returns the corresponding right operand, otherwise d, or false, if d is missing.

Check out the demonstration below for both the above operations.

$ maxima -q
(%i1) sublist_indices([12, 23, 57, 37, 64, 67], primep);
(%o1)                              [2, 4, 6]
(%i2) sublist_indices([12, 23, 57, 37, 64, 67], evenp);
(%o2)                               [1, 5]
(%i3) sublist_indices([12, 23, 57, 37, 64, 67], oddp);
(%o3)                            [2, 3, 4, 6]
(%i4) sublist_indices([2 > 0, -2 > 0, 1 = 1, x = y], identity);
(%o4)                               [1, 3]
(%i5) assoc(2, [2^r, x+y, 2=4, 5/6]);
(%o5)                                  r
(%i6) assoc(6, [2^r, x+y, 2=4, 5/6]);
(%o6)                               false
(%i7) assoc(6, [2^r, x+y, 2=4, 5/6], na);
(%o7)                                na
(%i8) quit();

Twenty-first Article >>

   Send article as PDF   

Advanced Set Theory

4 Replies

This nineteenth article of the mathematical journey through open source, explores the advanced set theory concepts through Maxima.

<< Eighteenth Article

With the introduction to set theory fundamentals in the previous article, we are all set to explore the advanced realms of set theory through Maxima.

More set operations

We have already worked out the basic set creation techniques and some basic set operations provided by Maxima. Here are some next-level ones provided by Maxima:

  • makeset(expr, vars, varslist) – Sophisticated set creation using expressions
  • adjoin(x, S) – Returns a set with all elements of S and the element x
  • disjoin(x, S) – Returns a set with all elements S but without element x
  • powerset(S) – Returns the set of all subsets of S
  • subset(S, p) – Returns the subset of S, elements of which satisfy the predicate p
  • symmdifference(S1, S2) – Returns the symmetric difference between the sets S1 and S2, i.e. the elements in S1 or S2 but not in both

And here goes a demonstration of each one of them:

$ maxima -q
(%i1) makeset(a+b, [a, b], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6]]);
(%o1)                          {3, 5, 7, 9, 11}
(%i2) makeset(a-b, [a, b], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6]]);
(%o2)                                {- 1}
(%i3) makeset(a*b, [a, b], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6]]);
(%o3)                         {2, 6, 12, 20, 30}
(%i4) makeset(a + 2*a*b + b, [a, b], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6]]);
(%o4)                         {7, 17, 31, 49, 71}
(%i5) quit();
$ maxima -q
(%i1) S: {-4, 6, 7, 32, 0};
(%o1)                         {- 4, 0, 6, 7, 32}
(%i2) adjoin(3, S);
(%o2)                        {- 4, 0, 3, 6, 7, 32}
(%i3) adjoin(7, S);
(%o3)                        {- 4, 0, 6, 7, 32}
(%i4) S: adjoin(3, S); /* Updating S */
(%o4)                       {- 4, 0, 3, 6, 7, 32}
(%i5) adjoin(7, S);
(%o5)                       {- 4, 0, 3, 6, 7, 32}
(%i6) disjoin(7, S);
(%o6)                        {- 4, 0, 3, 6, 32}
(%i7) disjoin(5, S);
(%o7)                       {- 4, 0, 3, 6, 7, 32}
(%i8) quit();
$ maxima -q
(%i1) S: {-4, 0, 3, 6, 7, 32};
(%o1)                        {- 4, 0, 3, 6, 7, 32}
(%i2) S1: subset(S, evenp);
(%o2)                           {- 4, 0, 6, 32}
(%i3) powerset(S1);
(%o3) {{}, {- 4}, {- 4, 0}, {- 4, 0, 6}, {- 4, 0, 6, 32}, {- 4, 0, 32}, {- 4, 6},
	{- 4, 6, 32}, {- 4, 32}, {0}, {0, 6}, {0, 6, 32}, {0, 32}, {6}, {6, 32},
	{32}}
(%i4) S2: {-35, -26, 0, 7, 32, 100};
(%o4)                     {- 35, - 26, 0, 7, 32, 100}
(%i5) symmdifference(S1, S2);
(%o5)                    {- 35, - 26, - 4, 6, 7, 100}
(%i6) symmdifference(S, S2);
(%o6)                    {- 35, - 26, - 4, 3, 6, 100}
(%i7) quit();

Advanced set operations

With Maxima, much more than this can be done with sets, just by the functionalities provided by it. So now, let’s journey through such interesting advance functionalities:

Cartesian Product:

Given n sets, the function cartesian_product() returns a set of lists formed by the Cartesian product of the n sets. The following demonstration explains what it means:

$ maxima -q
(%i1) cartesian_product({0, 1, 2}, {a, b, c});
(%o1) {[0, a], [0, b], [0, c], [1, a], [1, b], [1, c], [2, a], [2, b], [2, c]}
(%i2) cartesian_product({0, 1}, {a, b}, {X, Y});
(%o2) {[0, a, X], [0, a, Y], [0, b, X], [0, b, Y], [1, a, X], [1, a, Y], [1, b, X],
	[1, b, Y]}
(%i3) cartesian_product({0, 1}, {a, b, c});
(%o3)          {[0, a], [0, b], [0, c], [1, a], [1, b], [1, c]}
(%i4) quit();

Set Partitions:

Given a set S, it could be partitioned into various subsets, based on various mathematical principles. Maxima provides a host of functions for such partitioning. The basic one being set_partitions(). It returns a set of all possible partitions of the given set. With a number as the second argument, it gives only the partitions with that exact number of sets in each partition. Examples follow to make sense out of this:

$ maxima -q
(%i1) S: {a, b, c};
(%o1)                              {a, b, c}
(%i2) set_partitions(S);
(%o2) {{{a}, {b}, {c}}, {{a}, {b, c}}, {{a, b}, {c}}, {{a, b, c}}, {{a, c}, {b}}}
(%i3) set_partitions(S, 1);
(%o3)                            {{{a, b, c}}}
(%i4) set_partitions(S, 2);
(%o4)            {{{a}, {b, c}}, {{a, b}, {c}}, {{a, c}, {b}}}
(%i5) set_partitions(S, 3);
(%o5)                          {{{a}, {b}, {c}}}
(%i6) set_partitions(S, 4);
(%o6)                                 {}
(%i7) belln(3);
(%o7)                                  5
(%i8) cardinality(set_partitions(S)); /* Number of elements */
(%o8)                                  5
(%i9) belln(4);
(%o9)                                 15
(%i10) belln(5);
(%o10)                                 52
(%i11) belln(6);
(%o11)                                 203
(%i12) quit();

In the above examples, belln() – the nth Bell number is the number of partitions of a set with n members.

integer_partitions(n) is a specific function, which partitions a given positive integer n into set of positive integers, sum of which adds up to the original integer. num_partitions(n) returns the number of such partitions returned by integer_partitions(n). Examples follow:

$ maxima -q
(%i1) integer_partitions(1);
(%o1)                                {[1]}
(%i2) num_partitions(1);
(%o2)                                 1
(%i3) integer_partitions(2);
(%o3)                            {[1, 1], [2]}
(%i4) num_partitions(2);
(%o4)                                 2
(%i5) integer_partitions(3);
(%o5)                      {[1, 1, 1], [2, 1], [3]}
(%i6) num_partitions(3);
(%o6)                                 3
(%i7) integer_partitions(4);
(%o7)           {[1, 1, 1, 1], [2, 1, 1], [2, 2], [3, 1], [4]}
(%i8) num_partitions(4);
(%o8)                                 5
(%i9) integer_partitions(0);
(%o9)                                {[]}
(%i10) num_partitions(0);
(%o10)                                 1
(%i11) integer_partitions(5, 1);
(%o11)                                {[5]}
(%i12) integer_partitions(5, 2);
(%o12)                      {[3, 2], [4, 1], [5, 0]}
(%i13) integer_partitions(5, 3);
(%o13)       {[2, 2, 1], [3, 1, 1], [3, 2, 0], [4, 1, 0], [5, 0, 0]}
(%i14) integer_partitions(5, 4);
(%o14) {[2, 1, 1, 1], [2, 2, 1, 0], [3, 1, 1, 0], [3, 2, 0, 0], [4, 1, 0, 0],
	[5, 0, 0, 0]}
(%i15) integer_partitions(5, 5);
(%o15) {[1, 1, 1, 1, 1], [2, 1, 1, 1, 0], [2, 2, 1, 0, 0], [3, 1, 1, 0, 0],
	[3, 2, 0, 0, 0], [4, 1, 0, 0, 0], [5, 0, 0, 0, 0]}
(%i16) num_partitions(5);
(%o16)                                 7
(%i17) num_distinct_partitions(5);
(%o17)                                 3
(%i18) quit();

Note that like set_partitions(), integer_partitions() also takes an optional second argument, limiting the partitions to partitions of cardinality equal to that number. However, note that all smaller size partitions are made equal to the corresponding size by adding the required number of zeroes.

Also, num_distinct_partitions(n) returns number of distinct integer partitions of n, i.e. integer partitions of n with only distinct integers.

Another powerful partitioning function is the function equiv_classes(S, r), which returns a partition of S, elements of which satisfy the binary relation r. Here goes a few examples:

$ maxima -q
(%i1) equiv_classes({0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, lambda([x, y],
	remainder(x - y, 2) = 0));
(%o1)               {{0, 2, 4, 6, 8, 10}, {1, 3, 5, 7, 9}}
(%i2) equiv_classes({0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, lambda([x, y],
	remainder(x - y, 3) = 0));
(%o2)              {{0, 3, 6, 9}, {1, 4, 7, 10}, {2, 5, 8}}
(%i3) equiv_classes({0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, lambda([x, y],
	remainder(x - y, 5) = 0));
(%o3)            {{0, 5, 10}, {1, 6}, {2, 7}, {3, 8}, {4, 9}}
(%i4) equiv_classes({0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, lambda([x, y],
	remainder(x - y, 6) = 0));
(%o4)           {{0, 6}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {5}}
(%i5) equiv_classes({1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, lambda([x, y],
	remainder(x, y) = 0));
(%o5)             {{1, 2, 4, 8}, {3, 6}, {5, 10}, {7}, {9}}
(%i6) quit();

Notice the relation being defined using lamda() for the property of divisibility by 2, 3, 5, 6, and among the set elements themselves, respectively.

A closely related function partition_set(S, p) partitions S into 2 sets, one with elements satisfying the predicate p, and the other not satisfying the predicate p. Small demonstration follows:

$ maxima -q
(%i1) partition_set({-1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 19, 26, 37, 100},
	primep);
(%o1) [{- 1, 0, 1, 4, 6, 8, 9, 10, 26, 100}, {2, 3, 5, 7, 11, 19, 37}]
(%i2) quit();

Miscellaneous:

And finally, lets look at some general but mathematically interesting operations:

  • divisors(n) – returns the set of positive divisors of n
  • permutations(S) – returns the set of all permutations of the elements of S
  • random_permutation(S) – returns one of the elements of permutations(S), randomly
  • extremal_subset(S, f, max | min) – returns the subset of S, for which the value of the function f is maximum or minimum

An all-together demonstration follows:

$ maxima -q
(%i1) divisors(9);
(%o1)                             {1, 3, 9}
(%i2) divisors(28);
(%o2)                       {1, 2, 4, 7, 14, 28}
(%i3) permutations({a, b, c});
(%o3) {[a, b, c], [a, c, b], [b, a, c], [b, c, a], [c, a, b], [c, b, a]}
(%i4) random_permutation({a, b, c});
(%o4)                             [c, b, a]
(%i5) random_permutation({a, b, c});
(%o5)                              [c, a, b]
(%i6) random_permutation({a, b, c});
(%o6)                              [b, c, a]
(%i7) extremal_subset({-5, -3, -1, 0, 1, 2, 3, 4, 5}, lambda([x], x*x), max);
(%o7)                            {- 5, 5}
(%i8) extremal_subset({-5, -3, -1, 0, 1, 2, 3, 4, 5}, lambda([x], x*x), min);
(%o8)                               {0}
(%i9) quit();

Twentieth Article >>

   Send article as PDF   

Fundamentals of Set Theory

2 Replies

This eighteenth article of the mathematical journey through open source, digs into the fundamentals of the set theory through Maxima.

<< Seventeenth Article

Most of you would have definitely heard about the word “set” as in “set theory”, and many of you may already know quite a bit about the sets. Then what’s so special about it. The specialty lies in the whole idea of getting our familiar maths being done by the computer. In this series, we have been mostly picking up familiar maths and then figuring out how that can be done by computer, and that also with no or very little programming knowledge. The same holds true for the sets as well. So, let’s get started with its fundamentals, starting from their creation.

Creation of sets

Just for the beginners of set theory, a set is an unordered collection of distinct items – any item, in any order, but unique. An item is commonly referred to as an element, as well. If an item is contained in a set, it is commonly referred as a member of the set. A set is typically represented by its members enclosed in braces {} and separated by commas. {6, -5, 9, 0}, {dog, cat, donkey, cow, buffalo}, {Kiran, Kasturi, Karan}, {6, horse, sapphire} are some examples. Notice that the first three sets have related items in them, but the last one doesn’t – and that’s perfectly fine. However, if the items in a set have relation(s) or condition(s), the set can also be expressed with that relation(s) or condition(s) mentioned within braces {}. For example, {All human beings younger than 35 years}, {All positive even numbers, All multiples of 3}. In Maxima, we can straight away represent the sets in the first notation, as follows:

$ maxima -q
(%i1) {6, -5, 9, 0};
(%o1)                           {- 5, 0, 6, 9}
(%i2) {dog, cat, donkey, cow, buffalo};
(%o2)                  {buffalo, cat, cow, dog, donkey}
(%i3) {Kiran, Kasturi, Karan};
(%o3)                       {Karan, Kasturi, Kiran}
(%i4) {6, horse, sapphire};
(%o4)                        {6, horse, sapphire}
(%i5) {axe, knife, spear, axe, scissor};
(%o5)                    {axe, knife, scissor, spear}
(%i6) quit();

Note that as the order of items in the set doesn’t matter, Maxima internally keeps them sorted, and hence displayed accordingly, in the above examples. Also note the last example – the duplicates are treated as a single item.

Sets can also be created from ordered lists using setify(). Items of sets could be expressions, but may not be automatically simplified. Check out the following:

$ maxima -q
(%i1) setify([x, y, z]);
(%o1)                              {x, y, z}
(%i2) [x, y, z, x]; /* Ordered list */
(%o2)                            [x, y, z, x]
(%i3) setify([x, y, z, x]);
(%o3)                              {x, y, z}
(%i4) string({x^2 - 1, (x + 1) * (x -1)});
(%o4)                         {(x-1)*(x+1),x^2-1}

string() has been used in %i4, to just have the output on a single line. But the important thing to note is that though the two items of the list are mathematically identical, they have been preserved as two distinct items – and thus not a set in real sense. Such cases can be actually setified by simplifying the individual items of the set using corresponding simplification function, e.g. rat() for rational expressions. And operating any function on every item of a set can be achieved using map(). Here’s an to example to set all those straight, continuing from the above one:

(%i5) string(map(rat, {x^2 - 1, (x + 1) * (x -1)}));
(%o5)                               {x^2-1}
(%i6) string(rat((x + 1) * (x -1)));
(%o6)                                x^2-1
(%i7) quit();

%i6 and %o6 above is just to demonstrate how rat() works. I know, you are still wondering what this wierd map() is and how it is working. So, here are a few more examples for the same:

$ maxima -q
(%i1) trigreduce(2 * sin(x) * cos(x));
(%o1)                              sin(2 x)
(%i2) {sin(2 * x), 2 * sin(x) * cos(x)};  /* Identical items */
(%o2)                     {2 cos(x) sin(x), sin(2 x)}
(%i3) map(trigreduce, {sin(2 * x), 2 * sin(x) * cos(x)});
(%o3)                             {sin(2 x)}
(%i4) string({apple / fruit + mango / fruit, (apple + mango) / fruit});
(%o4)            {(mango+apple)/fruit,mango/fruit+apple/fruit}
(%i5) string(map(rat, {apple / fruit + mango / fruit, (apple + mango) / fruit}));
(%o5)                        {(mango+apple)/fruit}
(%i6) quit();

In fact, the power of map() lies in its ability that it can take a function created on the fly, meaning using the lambda notation. Here goes few examples to demonstrate the lamda() first, and then map() using lambda():

$ maxima -q
(%i1) f: lambda([x], x^3)$
(%i2) f(5);
(%o2)                                 125
(%i3) lambda([x], x^3)(5);
(%o3)                                 125
(%i4) lambda([x, y], x+y)(4, 6);
(%o4)                                 10
(%i5) map(f, {0, 1, 2, 3});
(%o5)                            {0, 1, 8, 27}
(%i6) map(lambda([x], x^3), {0, 1, 2, 3});
(%o6)                            {0, 1, 8, 27}
(%i7) map(lambda([x, y], x+y), {a}, {3});
(%o7)                               {a + 3}
(%i8) map(lambda([x], x^4), {-2, -1, 0, 1, 2});
(%o8)                            {0, 1, 16}
(%i9) map(g, {-2, -1, 0, 1, 2});
(%o9)                {g(- 2), g(- 1), g(0), g(1), g(2)}
(%i10) quit();

lambda() takes two arguments. First a list of arguments of the function being defined, and second the expression for the return value of the function using those arguments. %i1 defines a function f with one argument, returning its cube. %i2 calls f(). However, the whole point of using lambda is using it without defining an explicit function like f(). So, %i3 & %i4 demonstrate exactly that. %i6, %i7, %i8 shows using lambda() with map(). Note the elimination of duplicates in %o8. %i9 is an another example of map().

Basic set operations

Enough with varieties of set creation. Let’s do some set operations. For the starters: union of sets is defined as a set with items of all the sets; intersection of sets is defined as a set with items common in all the sets; difference of two sets is defined as a set with items from the first set, not in the second set. And here goes the demonstration:

$ maxima -q
(%i1) union({1, 2}, {1, 3, 4}, {1, 2, 6, 7});
(%o1)                         {1, 2, 3, 4, 6, 7}
(%i2) intersection({1, 2}, {1, 3, 4}, {1, 2, 6, 7});
(%o2)                                 {1}
(%i3) setdifference({1, 2}, {1, 3, 4});
(%o3)                                 {2}
(%i4) quit();

Other basic set operations provided by Maxima are:-

  • cardinality() – returns the number of distinct items in a set
  • elementp() – checks for an item to be a member of a set
  • emptyp() – checks for the emptiness of a set
  • setequalp() – compares two sets for equality
  • disjointp() – checks for no common items in two sets
  • subsetp() – checks for the first set to be a subset of the second set

 

The following walk through demonstrates all of these:

$ maxima -q
(%i1) S1: {}$
(%i2) S2: {1, 2, 3}$
(%i3) S3: {3, 1, 5-3}$ /* Same as S2 */
(%i4) S4: {a, b, c}$
(%i5) S5: {2, 1, 2}$
(%i6) cardinality(S1);
(%o6)                                  0
(%i7) cardinality(S2);
(%o7)                                  3
(%i8) cardinality(S3);
(%o8)                                  3
(%i9) cardinality(S4);
(%o9)                                  3
(%i10) cardinality(S5);
(%o10)                                 2
(%i11) elementp(b, S3);
(%o11)                               false
(%i12) elementp(b, S4);
(%o12)                               true
(%i13) emptyp(S1);
(%o13)                               true
(%i14) emptyp(S2);
(%o14)                               false
(%i15) setequalp(S1, S2);
(%o15)                               false
(%i16) setequalp(S2, S3);
(%o16)                               true
(%i17) disjointp(S1, S2);
(%o17)                               true
(%i18) disjointp(S2, S3);
(%o18)                               false
(%i19) disjointp(S3, S4);
(%o19)                               true
(%i20) disjointp(S3, S5);
(%o20)                               false
(%i21) subsetp(S1, S2);
(%o21)                               true
(%i22) subsetp(S2, S3);
(%o22)                               true
(%i23) subsetp(S3, S2);
(%o23)                               true
(%i24) subsetp(S3, S4);
(%o24)                               false
(%i25) subsetp(S5, S3);
(%o25)                               true
(%i26) subsetp(S3, S5);
(%o26)                               false
(%i27) quit();

Playing with set elements

After clearing the fundamentals, mostly through numerical examples, it is now time to have some fun with symbol substitution of Maxima. And here’s some playing around:

$ maxima -q
(%i1) S: {a, b, c, a};
(%o1)                             {a, b, c}
(%i2) S: {a+b, b+c, c+d, d+a};
(%o2)                   {b + a, c + b, d + a, d + c}
(%i3) subst(a=c, S);
(%o3)                          {c + b, d + c}
(%i4) subst([a=c, b=d], S);
(%o4)                              {d + c}
(%i5) subst([a=c, b=d, c=-d], S);
(%o5)                                {0}
(%i6) subst([a=1, b=2, c=-3], S);
(%o6)                       {- 1, 3, d - 3, d + 1}
(%i7) T: {S, {S}};
(%o7)   {{b + a, c + b, d + a, d + c}, {{b + a, c + b, d + a, d + c}}}
(%i8) subst([a=c, b=d, c=-d], T);
(%o8)                            {{0}, {{0}}}
(%i9) subst([a=1, b=2, c=-3], T);
(%o9)         {{- 1, 3, d - 3, d + 1}, {{- 1, 3, d - 3, d + 1}}}
(%i10) quit();

Nineteenth Article >>

   Send article as PDF   

High School Trigo with Maxima

2 Replies

This seventeenth article of the mathematical journey through open source, takes through a tour of the high school trigonometry using Maxima.

<< Sixteenth Article

Trigonometry first gets introduced to students of standard IX, through triangles. And, then it takes its own journey through the jungle of formulae and tables. And one knows, being good at instant recall of various formulae, makes her good at trigo. The idea here is not to be good at mugging up the formulae but rather applying them to get the various end results. It assumes that you possibly already know the formulae.

Fundamental trigonometric functions

Maxima provides all the familiar fundamental trigonometric functions, including the hyperbolic ones. They can be tabulated as follows:

Mathematical Names

Normal

Hyperbolic

Functions

Inv. Functions

Functions

Inv. Functions
sine (sin)

sin()

asin()

sinh()

asinh()
cosine (cos)

cos()

acos()

cosh()

acosh()
tangent (tan)

tan()

atan()

tanh()

atanh()
cosecant (cosec)

csc()

acsc()

csch()

acsch()
secant (sec)

sec()

asec()

sech()

asech()
cotangent (cot)

cot()

acot()

coth()

acoth()

Note that all of their arguments are values in radians. And here follows a demonstration of a small subset of those:

$ maxima -q
(%i1) cos(0);
(%o1)                                  1
(%i2) cos(%pi/2);
(%o2)                                  0
(%i3) cot(0);
The number 0 isn't in the domain of cot
 -- an error. To debug this try: debugmode(true);
(%i4) tan(%pi/4);
(%o4)                                  1
(%i5) string(asin(1));
(%o5)                                %pi/2
(%i6) csch(0);
The number 0 isn't in the domain of csch
 -- an error. To debug this try: debugmode(true);
(%i7) csch(1);
(%o7)                               csch(1)
(%i8) asinh(0);
(%o8)                                  0
(%i9) string(%i * sin(%pi / 3)^2 + cos(5 * %pi / 6));
(%o9)                         3*%i/4-sqrt(3)/2
(%i10) quit();

Simplifications with special angles like %pi / 10 and its multiples could be enabled by loading the ntrig package. Check out the difference below, before and after loading the package:

$ maxima -q
(%i1) string(sin(%pi/10));
(%o1)                             sin(%pi/10)
(%i2) string(cos(2*%pi/10));
(%o2)                             cos(%pi/5)
(%i3) string(tan(3*%pi/10));
(%o3)                            tan(3*%pi/10)
(%i4) load(ntrig);
(%o4)        /usr/share/maxima/5.24.0/share/trigonometry/ntrig.mac
(%i5) string(sin(%pi/10));
(%o5)                            (sqrt(5)-1)/4
(%i6) string(cos(2*%pi/10));
(%o6)                            (sqrt(5)+1)/4
(%i7) string(tan(3*%pi/10));
(%o7)          sqrt(2)*(sqrt(5)+1)/((sqrt(5)-1)*sqrt(sqrt(5)+5))
(%i8) quit();

A very common trigonometric problem is, given a tangent value find the corresponding angle. Now, a common challenge is for every value, the angle could lie in two quadrants. For a positive tangent, the angle could be in the first or the third quadrant, and for a negative value, the angle could be in the second or the fourth quadrant. So, atan() cannot always calculate the correct quadrant of the angle. Now, how do we know it exactly. Obviously, we need some extra information, say the actual values of the perpendicular (p) and the base (b) of the tangent, rather than just the tangent value. With that, the angle location could be tabulated as follows:

Perpendicular (p)

Base (b)

Tangent (p/b)

Angle Quadrant

Positive

Positive

Positive

First

Positive

Negative

Negative

Second

Negative

Negative

Positive

Third

Negative

Positive

Negative

Fourth

And this functionality is captured in the atan2() function, which takes 2 arguments, the p and the b, and thus does provide the angle in the correct quadrant, as per the table above. Along with this, the infinities of tangent are also taken care. Here’s a demo:

$ maxima -q
(%i1) atan2(0, 1); /* Zero */
(%o1)                                  0
(%i2) atan2(0, -1); /* Zero */
(%o2)                                 %pi
(%i3) string(atan2(1, -1)); /* -1 */
(%o3)                               3*%pi/4
(%i4) string(atan2(-1, -1)); /* 1 */
(%o4)                              -3*%pi/4
(%i5) string(atan2(-1, 0)); /* - Infinity */
(%o5)                               -%pi/2
(%i6) string(atan2(5, 0)); /* + Infinity */
(%o6)                                %pi/2
(%i7) quit();

Trigonometric Identities

Maxima supports many built-in trigonometric identities, and one can add his own as well. First one to look at would be the set dealing with integral multiples and factors of %pi. We’ll declare a few integers and then play around with them.

$ maxima -q
(%i1) declare(m, integer, n, integer);
(%o1)                                done
(%i2) properties(m);
(%o2)                  [database info, kind(m, integer)]
(%i3) sin(m * %pi);
(%o3)                                  0
(%i4) string(cos(n * %pi));
(%o4)                              (-1)^n
(%i5) string(cos(m * %pi / 2)); /* No simplification */
(%o5)                            cos(%pi*m/2)
(%i6) declare(m, even); /* Will lead to simplification */
(%o6)                               done
(%i7) declare(n, odd);
(%o7)                               done
(%i8) cos(m * %pi);
(%o8)                                 1
(%i9) cos(n * %pi);
(%o9)                                - 1
(%i10) string(cos(m * %pi / 2));
(%o10)                             (-1)^(m/2)
(%i11) string(cos(n * %pi / 2));
(%o11)                            cos(%pi*n/2)
(%i12) quit();

Next is the relation between the normal & the hyperbolic trigo functions.

$ maxima -q
(%i1) sin(%i * x);
(%o1)                            %i sinh(x)
(%i2) cos(%i * x);
(%o2)                              cosh(x)
(%i3) tan(%i * x);
(%o3)                            %i tanh(x)
(%i4) quit();

By enabling the option variable halfangles, many half angle identities, come into play. To be specific, sin(x/2) gets further simplified in (0, 2 * %pi) range, cos(x/2) gets further simplified in (-%pi/2, %pi/2) range. Check out the differences, before and after enabling the option variable, along with the range modifications, in the examples below:

$ maxima -q
(%i1) string(2*cos(x/2)^2 - 1); /* No effect */
(%o1)                           2*cos(x/2)^2-1
(%i2) string(cos(x/2)); /* No effect */
(%o2)                              cos(x/2)
(%i3) halfangles:true; /* Enabling half angles */
(%o3)                                true
(%i4) string(2*cos(x/2)^2 - 1); /* Simplified */
(%o4)                               cos(x)
(%i5) string(cos(x/2)); /* Complex expansion for all x */
(%o5)         (-1)^floor((x+%pi)/(2*%pi))*sqrt(cos(x)+1)/sqrt(2)
(%i6) assume(-%pi < x, x < %pi); /* Limiting x values */
(%o6)                        [x > - %pi, x < %pi]
(%i7) string(cos(x/2)); /* Further simplified */
(%o7)                       sqrt(cos(x)+1)/sqrt(2)
(%i8) quit();

Trigonometric Expansions and Simplifications

Trigonometry is full of multiples of angles, sums of angles, products & powers of trigo functions, and the long list of relations between them. Multiples and sums of angles fall into one category. Products and powers of trigo functions in an another category. And its very useful to do conversions from one of these categories to the other one, to crack a range of simple to complex problems, catering to basic hobby science to quantum mechanics. trigexpand() does the conversion from “multiples & sums of angles” to “products & powers of trigo functions”. trigreduce() does exactly the opposite. Here’s goes a small demo:

$ maxima -q
(%i1) trigexpand(sin(2*x));
(%o1)                           2 cos(x) sin(x)
(%i2) trigexpand(sin(x+y)-sin(x-y));
(%o2)                           2 cos(x) sin(y)
(%i3) trigexpand(cos(2*x+y)-cos(2*x-y));
(%o3)                         - 2 sin(2 x) sin(y)
(%i4) trigexpand(%o3);
(%o4)                      - 4 cos(x) sin(x) sin(y)
(%i5) string(trigreduce(%o4));
(%o5)                  -2*(cos(y-2*x)/2-cos(y+2*x)/2)
(%i6) string(trigsimp(%o5));
(%o6)                       cos(y+2*x)-cos(y-2*x)
(%i7) string(trigexpand(cos(2*x)));
(%o7)                         cos(x)^2-sin(x)^2
(%i8) string(trigexpand(cos(2*x) + 2*sin(x)^2));
(%o8)                         sin(x)^2+cos(x)^2
(%i9) trigsimp(trigexpand(cos(2*x) + 2*sin(x)^2));
(%o9)                                 1
(%i10) quit();

In %o5 above, you might have noted that the 2’s could have been cancelled for further simplification. But that is not the job of trigreduce(), and for that we have to apply the trigsimp() function as shown in %i6. In fact, many other trigonometric identities based simplification are achieved using trigsimp(). Check out the %i7 to %o9 sequences for another such example.

Eighteenth Article >>

   Send article as PDF   

Polynomials in Maxima

4 Replies

This sixteenth article of the mathematical journey through open source, demonstrates the polynomial manipulations using Maxima.

<< Fifteenth Article

Polynomials have fascinated mathematicians since ages. Moreover, because of wide variety of their applications, ranging from basic algebra to puzzles to various sciences. Here, we are going to look at some of the polynomial manipulation functions provided by Maxima, and a couple of real world applications using some of them.

Fundamental polynomial operations

Let’s start with a demonstration of the fundamental polynomial operations, like addition, subtraction, multiplication, division. In all these, whenever needed,we’ll use expand() to expand the polynomials, and string() to display the polynomials in a flattened notation.

$ maxima -q
(%i1) p1: x^2 - y^2 + y$
(%i2) p2: -x^2 - y^2 + x$
(%i3) p3: (x + y)^2$
(%i4) string(p1 + p2);
(%o4)                             -2*y^2+y+x
(%i5) string(p1 + p2 + p3);
(%o5)                          (y+x)^2-2*y^2+y+x
(%i6) string(expand(p1 + p2 + p3));
(%o6)                         -y^2+2*x*y+y+x^2+x
(%i7) string(expand(p3 - p1));
(%o7)                            2*y^2+2*x*y-y
(%i8) string(expand(p1 * p2));
(%o8)                   y^4-y^3-x*y^2-x^2*y+x*y-x^4+x^3
(%i9) string(p1 / p2);
(%o9)                     (-y^2+y+x^2)/(-y^2-x^2+x)
(%i10) string(divide(p1, p2));
(%o10)                           [1,y+2*x^2-x]
(%i11) string(divide(x^2 - y^2, x + y));
(%o11)                              [x-y,0]
(%i12) quit();

Note that the / operator just places the polynomials as a fraction, rather then dividing them. And the function divide() actually divides the first polynomial by the second one, yielding a list with the quotient and the remainder of the division. If the division needs to be with respect to a particular variable, that can be passed as the third argument to divide. Check out the variation below, to understand what it means:

$ maxima -q
(%i1) string(divide(x^2 - y^2 + y, x + y, x));
(%o1)                              [x-y,y]
(%i2) string(divide(x^2 - y^2 + y, x + y, y));
(%o2)                            [-y+x+1,-x]
(%i3) quit();

Coefficients of a polynomial

Extracting the coefficients of a polynomial is another basic and common requirement for polynomial manipulations. Maxima provides two slightly different mechanisms. First one, just finds the coefficient of a given variable or its some power, using coeff(). Second one, segregates a polynomial into the coefficient of a given variable or its power, and the remaining terms, using bothcoef().

$ maxima -q
(%i1) string(bothcoef(expand(x^2 - y^2 + (x + y)^2), x^2));
(%o1)                             [2,2*x*y]
(%i2) string(bothcoef(expand(x^2 - y^2 + (x + y)^2), x));
(%o2)                            [2*y,2*x^2]
(%i3) string(bothcoef(expand(x^2 - y^2 + (x + y)^2), y^2));
(%o3)                          [0,2*x*y+2*x^2]
(%i4) string(bothcoef(expand(x^2 - y^2 + (x + y)^2), y));
(%o4)                            [2*x,2*x^2]
(%i5) string(coeff(expand((x + 2 * y)^50), x^20));
(%o5)                   50604606318512743383040*y^30
(%i6) string(coeff(expand((a + b + c + d)^4), a^3));
(%o6)                            4*d+4*c+4*b
(%i7) quit();

Polynomial fractions

Calculating the greatest common divisor (GCD) is one of the very useful operations to simply fractions of polynomials. Other common requirements are extracting the numerator, the denominator, and the highest power. Here goes the function demonstrations for all of these:

$ maxima -q
(%i1) gcd(x^3 + 3*x^2 + 3*x + 1, x^2 + 3*x + 2);
(%o1)                              x + 1
(%i2) string(ezgcd(x^3 + 3*x^2 + 3*x + 1, x^2 + 3*x + 2));
(%o2)                       [x+1,x^2+2*x+1,x+2]
(%i3) string(denom((x + 1)^-3 * (1 - x)^2));
(%o3)                             (x+1)^3
(%i4) string(num((x + 1)^-3 * (1 - x)^2));
(%o4)                             (1-x)^2
(%i5) hipow(expand((x + 1)^3 + (1 - x)^3), x);
(%o5)                                2
(%i6) quit();

Note that the ezgcd() function lists out the remainder polynomials, along with the GCD.

Polynomial fractions could be differentiated using the powerful ratdiff().

$ maxima -q
(%i1) string(ratdiff((x + 1)^-1 * (1 - x)^2, x));
(%o1)                     (x^2+2*x-3)/(x^2+2*x+1)
(%i2) string(ratdiff(1 / (x + 1), x));
(%o2)                         -1/(x^2+2*x+1)
(%i3) string(ratdiff((x^2 - 1) / (x + 1), x));
(%o3)                                1
(%i4) quit();

And, ratsubst() is a powerful expression substitution function, with intelligence. It would dig into the expression to simplify complicated expressions, including trigonometric ones. Check out the %i5, for one of its powerful application. ratsubst(<new>, <old>, <expr>) replaces the <old> expression by the <new> expression in the complete expression <expr>.

$ maxima -q
(%i1) string(ratsubst(u, x^2, x^3 + 3*x^2 + 3*x + 1));
(%o1)                          (u+3)*x+3*u+1
(%i2) string(ratsubst(u, x^2, (x+1)^3));
(%o2)                          (u+3)*x+3*u+1
(%i3) string(ratsubst(u, x^2, (x+1)^4));
(%o3)                       (4*u+4)*x+u^2+6*u+1
(%i4) string(ratsubst(u, x - 1, x^4 - 2*x^2 + 1));
(%o4)                         u^4+4*u^3+4*u^2
(%i5) string(ratsubst(sin(x)^2, 1 - cos(x)^2, cos(x)^4 - 2*cos(x)^2 + 1));
(%o5)                            sin(x)^4
(%i5) quit();

Variable eliminations & equation solving

Very often, we come across N sets of equations in M sets of unknowns, where M >= N. If M = N, then most likely there exists a unique solution. However, if M > N, then there may be many possible solutions, but with some constraints. And in such case, it would be helpful to deduce some simpler set of expressions. This can be achieved using eliminate() of Maxima. Let’s have two polynomial expressions in 3 variables x, y, z to demonstrate the same.

$ maxima -q
(%i1) p1: x^2 + 2*x*y + z^2$
(%i2) p2: x + y + z$
(%i3) string(eliminate([p1, p2], [x]));
(%o3)                            [2*z^2-y^2]
(%i4) string(eliminate([p1, p2], [y]));
(%o4)                         [-z^2+2*x*z+x^2]
(%i5) string(eliminate([p1, p2], [z]));
(%o5)                         [y^2+4*x*y+2*x^2]
(%i6) quit();

Note that in all the above polynomial expressions, the expressions are assumed to be equated to zero. A beautiful application of the above is solving recurrence relations. Let’s say we have a non-linear equation given by Yn+1 = r * Yn * (1 – Yn). And, we want to solve it for a scenario that as n tends to infinity, the value of Y oscillates between two distinct values. This means that Yn+2 = Yn. Hence, we would have two expressions with three unknowns to solve with, namely:

  1. Yn+1 = r * Yn * (1 – Yn)
  2. Yn = r * Yn+1 * (1 – Yn+1)

So, representing Yn+1 with yn1 and Yn by yn, we may use solve() to solve for yn & yn1 in terms of r. But, if rather we want to get a feel of the equation which pops up, in terms of yn, we would have to use eliminate() to eliminate yn1.

$ maxima -q
(%i1) p1: yn1 = r * yn * (1 - yn)$
(%i2) p2: yn = r * yn1 * (1 - yn1)$
(%i3) string(eliminate([p1, p2], [yn1]));
(%o3)          [yn*(r^3*yn^3-2*r^3*yn^2+(r^3+r^2)*yn-r^2+1)]
(%i4) string(solve([p1, p2], [yn, yn1])); /* Just to show the solution */
(%o4) [[yn = (r-1)/r,yn1 = (r-1)/r],
    [yn = -(sqrt(r^2-2*r-3)-r-1)/(2*r),yn1 = (sqrt(r-3)*sqrt(r+1)+r+1)/(2*r)],
    [yn = (sqrt(r^2-2*r-3)+r+1)/(2*r),yn1 = -(sqrt(r-3)*sqrt(r+1)-r-1)/(2*r)],
    [yn = 0,yn1 = 0]]
(%i5) quit()

Seventeenth Article >>

   Send article as PDF   

Expression Simplification using Maxima

2 Replies

This fifteenth article of the mathematical journey through open source, demonstrates the various techniques of expression simplification using Maxima.

<< Fourteenth Article

Expression simplification can be done in a variety of ways. Let’s start with the simple ones, and then move onto more powerful ones.

Real number simplification

A simple simplification is to convert a given number into a rational number using ratsimp(). Below follows the demonstration.

$ maxima -q
(%i1) ratsimp(9);
(%o1)                                  9
(%i2) ratsimp(10.0);
rat: replaced 10.0 by 10/1 = 10.0
(%o2)                                 10
(%i3) string(ratsimp(4.2)); /* string to print it on one line */
rat: replaced 4.2 by 21/5 = 4.2
(%o3)                               21/5
(%i4) string(factor(3.2)); /* string to print it on one line */
rat: replaced 3.2 by 16/5 = 3.2
(%o4)                                2^4/5
(%i5) string(ratsimp(4.3333)); /* string to print it on one line */
rat: replaced 4.3333 by 43333/10000 = 4.3333
(%o5)                            43333/10000
(%i6) string(ratsimp(4.3333333)); /* string to print it on one line */
rat: replaced 4.3333333 by 13/3 = 4.333333333333333
(%o6)                               13/3
(%i7) quit();

Another one is to check whether a number is an integer using askinteger(). And if yes, is it an even or odd number, again using askinteger(). Moreover, asksign() checks for the sign. In case of trying these with unknown variables, Maxima would ask the user the necessary question(s) to deduce the response, and store it for future analysis. For example, askinteger(x) would ask us if x is an integer. And, if we say yes, it can then deduce many more information by itself. Below are some examples:

$ maxima -q
(%i1) askinteger(1);
(%o1)                                 yes
(%i2) askinteger(1.0);
rat: replaced 1.0 by 1/1 = 1.0
(%o2)                                 yes
(%i3) askinteger(1.2);
rat: replaced 1.2 by 6/5 = 1.2
(%o3)                                 no
(%i4) askinteger(-9);
(%o4)                                 yes
(%i5) askinteger(2/3 + 3/4 + 1/6 + 5/12);
(%o5)                                 yes
(%i6) askinteger(-9, even);
(%o6)                                 no
(%i7) askinteger(0, even);
(%o7)                                 yes
(%i8) properties(x);
(%o8)                                []
(%i9) askinteger(x);
Is x an integer?
yes; /* This is our response */
(%o9)                                 yes
(%i10) askinteger(x + 9);
(%o10)                                 yes
(%i11) askinteger(2 * x, even);
(%o11)                                yes
(%i12) askinteger(2 * x + 1, even);
(%o12)                                no
(%i13) askinteger(x, even);
Is x an even number?
n; /* This is our response */
(%o13)                                no
(%i14) askinteger(x, odd);
(%o14)                                yes
(%i15) askinteger(x^3 - (x + 1)^2, even);
(%o15)                                no
(%i16) properties(x);
(%o16)          [database info, kind(x, integer), kind(x, odd)]
(%i17) asksign((-1)^x);
(%o17)                                neg
(%i18) asksign((-1)^(x+1));
(%o18)                                pos
(%i19) asksign((-1)^x+1);
(%o19)                               zero
(%i20) quit();

Maxima simplification uses the concept of properties of symbols. As an example, note the properties of x in the above demonstration at %i8 and %i16.

Complex number simplification

Complex numbers have two common useful forms, namely the exponential and the circular (with sine & cosine) forms. demoivre() converts the exponential to circular form and exponentialize() does the other way round. Using expand() along with them, can simplify complicated looking expressions. Here goes few examples:

$ maxima -q
(%i1) string(demoivre(exp(%i*x^2))); /* %i is sqrt(-1) */
(%o1)                       %i*sin(x^2)+cos(x^2)
(%i2) string(exponentialize(a*(cos(t)))); /* %i is sqrt(-1) */
(%o2)                    a*(%e^(%i*t)+%e^-(%i*t))/2
(%i3) string(expand(exponentialize(a*(cos(t)+%i*sin(t))))); /* %i is sqrt(-1) */
(%o3)                            a*%e^(%i*t)
(%i4) quit();

Expansions and Reductions

As already seen in the previous article, expand() expands an expression completely, by default. However, it can be controlled by specifying the maximum power to which to expand for both the numerator and the denominator, respectively. Using factor() can compact the expanded expressions. Moreover, in many cases, we would like to expand it only with respect to only some variable(s). Say, (x + a + b)^2 should be expanded with respect to x. expandwrt() is meant exactly for that. One example of each of these is shown below.

$ maxima -q
(%i1) string(expand(((x+1)^2-x^2)/(x+1)^2, 2, 0));
(%o1)                      2*x/(x+1)^2+1/(x+1)^2
(%i2) string(factor(((x+1)^2-x^2)/(x+1)^2));
(%o2)                           (2*x+1)/(x+1)^2
(%i3) string(expandwrt((x+a+b)^2,x));
(%o3)                      x^2+2*(b+a)*x+(b+a)^2
(%i4) quit();

Expressions containing logs, exponentials, radicals (powers) can be simplified using radcan(). Rule based simplifications can be achieved using sequential comparative simplification function scsim(). Both of these call for a few examples.

$ maxima -q
(%i1) string(radcan(exp(5 * log(x) + log(3 * exp(log(y) / 4)))));
(%o1)                          3*x^5*y^(1/4)
(%i2) radcan((log(2*x+2*x^2)-log(x))/(log(1+1/x)+log(2*x)));
(%o2)                                1
(%i3) expr: a^2 + b^2 + c^2$
(%i4) eq1: a^2 + 2*a*b + b^2 = 4$
(%i5) eq2: a * b = 6$
(%i6) string(scsimp(expr, eq1, eq2));
(%o6)                              c^2-8
(%i7) quit();

Unlike Octave, Maxima by default doesn’t evaluate its expressions, it only simplifies. What it means is that expressions with integers like cos(1), exp(2), sqrt(3), etc. may remain as is in the most simplified form, instead of evaluating to their respective float numerical values. In such cases, we may force the evaluation by passing the option numer. Similar evaluation can be achieved for predicates, using pred.

$ maxima -q
(%i1) cos(1);
(%o1)                             cos(1)
(%i2) cos(1), numer;
(%o2)                        .5403023058681398
(%i3) sqrt(7);
(%o3)                             sqrt(7)
(%i4) sqrt(7), numer;
(%o4)                        2.645751311064591
(%i4) 1 + 2 > 9;
(%o4)                              3 > 9
(%i5) 1 + 2 > 9, pred;
(%o5)                              false
(%i6) string(%e^%pi < %pi^%e); /* string to print it on one line */
(%o6)                         %e^%pi < %pi^%e
(%i7) %e^%pi < %pi^%e, pred;
(%o7)                              false
(%i8) quit();

Summation Simplifications

Symbolical representation and manipulations of summations can be beautifully achieved using sum() and the option simpsum. Check out the code execution below:

$ maxima -q
(%i1) sum(a * i, i, 1, 5);
(%o1)                                15 a
(%i2) sum(a, i, 1, 5);
(%o2)                                 5 a
(%i3) sum(a * i, i, 1, n);
                                      n
                                     ====
                                     \
(%o3)                              a  >    i
                                     /
                                     ====
                                     i = 1
(%i4) simpsum: true$ /* Enable simplification of summations */
(%i5) string(sum(a * i, i, 1, n)); /* string to print it on one line */
(%o5)                            a*(n^2+n)/2
(%i6) string(sum(i^2 - i, i, 1, n) + sum(j, j, 1, n));
(%o6)                         (2*n^3+3*n^2+n)/6
(%i7) string(factor(sum(i^2 - i, i, 1, n) + sum(j, j, 1, n)));
(%o7)                         n*(n+1)*(2*n+1)/6
(%i8) quit();

Sixteenth Article >>

   Send article as PDF