Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
Download

Chapter 2

44 views
unlisted
ubuntu2004-eol
Kernel: SageMath 9.7

Chapter 2 - Sage & basic mathematical structures - A warmup


2.1. Solving equations and inequalities
2.2. One-variable functions and 2D-plotting
2.3. Basics on sets and combinatorics
2.4. Classical & conditional probability
2.5. Plane geometry
2.6. Relations & mappings
2.7. Difference equations

In this chapter, we will explore various methods in SageMath for solving different mathematical tasks, starting with an overview of how to tackle equations and inequalities. Given the diverse types of equations—such as polynomial, linear, difference, or differential equations (and their systems)—we will encounter additional techniques for solving these equations with SageMath later in the course. The material in the first section of this chapter focuses on introducing the reader to fundamental methods in SageMath for solving relatively simple equations, rather than delving into specialized types. We will then move on to plotting functions and discuss methods for providing graphical confirmation of our findings.

2.1. Solving equations and inequalities

In Chapter 1, we learned that Sage can solve equations (and inequalities) using the command solve(){\tt{solve()}}. In fact, this command can be effectively applied to systems of equations involving more than one variable. It's important to note that the solve(){\tt{solve()}} command aims to express the solution of an equation without resorting to floating-point numbers. If this is not possible, it will return the solution in symbolic form. To obtain a numeric approximation of the solution, we can use the find_root(){\tt{find\_root()}} command, among other techniques.

Exercise

Solve the following equations or inequalities in SageMath:

  1. 3x2+2x5=03x^2+2x-\sqrt{5}=0;

  2. 3x3+2x25x+6=03x^3+2x^2-5x+\sqrt{6}=0;

  3. x2783 x^2 -78 \geq 3

  4. x31=0x^3-1=0;

  5. ex=1e^x=-1;

  6. exx2=0e^x - x^2 = 0;

  7. sin(x)=x\sin(x) =x;

  8. log(x2)=6/5\log(x^2)=6/5.

Solution:

show( solve(3*x^2+2*x-sqrt(5)==0, x))

[x=1335+113,x=1335+113]\displaystyle \left[x = -\frac{1}{3} \, \sqrt{3 \, \sqrt{5} + 1} - \frac{1}{3}, x = \frac{1}{3} \, \sqrt{3 \, \sqrt{5} + 1} - \frac{1}{3}\right]

Or we use the command roots{\tt{roots}}, as follows:

f=3*x^2+2*x-sqrt(5) show(f.roots())

[(1335+113,1),(1335+113,1)]\displaystyle \left[\left(-\frac{1}{3} \, \sqrt{3 \, \sqrt{5} + 1} - \frac{1}{3}, 1\right), \left(\frac{1}{3} \, \sqrt{3 \, \sqrt{5} + 1} - \frac{1}{3}, 1\right)\right]

The second number inside the parenthesis indicates the multiplicity of the root (1 in our case).

As for the second equation, which is of degree 3, by the Fundamental Theorem of Algebra we know that there exist 3 complex solutions.

SageMath computes them very quickly:

solve(3*x^3+2*x^2-5*x+sqrt(6)==0, x)
[x == -1/2*(-1/6*sqrt(6) + 1/54*sqrt(572/3*sqrt(6) - 142/3) - 143/729)^(1/3)*(I*sqrt(3) + 1) - 49/162*(-I*sqrt(3) + 1)/(-1/6*sqrt(6) + 1/54*sqrt(572/3*sqrt(6) - 142/3) - 143/729)^(1/3) - 2/9, x == -1/2*(-1/6*sqrt(6) + 1/54*sqrt(572/3*sqrt(6) - 142/3) - 143/729)^(1/3)*(-I*sqrt(3) + 1) - 49/162*(I*sqrt(3) + 1)/(-1/6*sqrt(6) + 1/54*sqrt(572/3*sqrt(6) - 142/3) - 143/729)^(1/3) - 2/9, x == (-1/6*sqrt(6) + 1/54*sqrt(572/3*sqrt(6) - 142/3) - 143/729)^(1/3) + 49/81/(-1/6*sqrt(6) + 1/54*sqrt(572/3*sqrt(6) - 142/3) - 143/729)^(1/3) - 2/9]

However, as it seems from Sage's output, extracting or copying these solutions can be messy.

To overpass this problem, or at least encode in a better way the given solutions, one can work as follows:

sols = solve([3*x^3+2*x^2-5*x+sqrt(6)==0],x)

Then we obtain any solutions seperately, by typing sols[0],sols[1],sols[2]{\tt{sols[0]}}, {\tt{sols[1]}}, {\tt{sols[2]}} (recall that Sage's labeling in lists starts by 0).

sols[0]
x == -1/2*(-1/6*sqrt(6) + 1/54*sqrt(572/3*sqrt(6) - 142/3) - 143/729)^(1/3)*(I*sqrt(3) + 1) - 49/162*(-I*sqrt(3) + 1)/(-1/6*sqrt(6) + 1/54*sqrt(572/3*sqrt(6) - 142/3) - 143/729)^(1/3) - 2/9
sols[1]
x == -1/2*(-1/6*sqrt(6) + 1/54*sqrt(572/3*sqrt(6) - 142/3) - 143/729)^(1/3)*(-I*sqrt(3) + 1) - 49/162*(I*sqrt(3) + 1)/(-1/6*sqrt(6) + 1/54*sqrt(572/3*sqrt(6) - 142/3) - 143/729)^(1/3) - 2/9

In case we want the real/imaginry part of some solution, we can use the commands real(){\tt{real()}} and imag(){\tt{imag()}}.

real(sols[1])
real_part(x == 1/2*I*sqrt(3)*(-1/6*sqrt(6) + 1/54*sqrt(572/3*sqrt(6) - 142/3) - 143/729)^(1/3) - 1/2*(-1/6*sqrt(6) + 1/54*sqrt(572/3*sqrt(6) - 142/3) - 143/729)^(1/3) - 49/162*I*sqrt(3)/(-1/6*sqrt(6) + 1/54*sqrt(572/3*sqrt(6) - 142/3) - 143/729)^(1/3) - 49/162/(-1/6*sqrt(6) + 1/54*sqrt(572/3*sqrt(6) - 142/3) - 143/729)^(1/3) - 2/9)
imag(sols[1])
imag_part(x == 1/2*I*sqrt(3)*(-1/6*sqrt(6) + 1/54*sqrt(572/3*sqrt(6) - 142/3) - 143/729)^(1/3) - 1/2*(-1/6*sqrt(6) + 1/54*sqrt(572/3*sqrt(6) - 142/3) - 143/729)^(1/3) - 49/162*I*sqrt(3)/(-1/6*sqrt(6) + 1/54*sqrt(572/3*sqrt(6) - 142/3) - 143/729)^(1/3) - 49/162/(-1/6*sqrt(6) + 1/54*sqrt(572/3*sqrt(6) - 142/3) - 143/729)^(1/3) - 2/9)
sols[2]
x == (-1/6*sqrt(6) + 1/54*sqrt(572/3*sqrt(6) - 142/3) - 143/729)^(1/3) + 49/81/(-1/6*sqrt(6) + 1/54*sqrt(572/3*sqrt(6) - 142/3) - 143/729)^(1/3) - 2/9
sols[3] #we expect that this will not give something, since the list index 3 is out of our range
--------------------------------------------------------------------------- IndexError Traceback (most recent call last) Input In [42], in <cell line: 1>() ----> 1 sols[Integer(3)]
File /ext/sage/9.7/src/sage/structure/sequence.py:527, in Sequence_generic.__getitem__(self, n) 521 return Sequence(list.__getitem__(self, n), 522 universe = self.__universe, 523 check = False, 524 immutable = False, 525 cr = self.__cr) 526 else: --> 527 return list.__getitem__(self,n)
IndexError: list index out of range

Recall that Sage allows us to compute the ``length of a collection or list'' , via the command len{\tt{len}}. Hence, typing len(sols){\tt{len(sols)}} we expect that Sage will print out the number 3.

len(sols)
3

Of course, this fits with the fact that the initial equation was of degree 3 (Fundamental Theorem of Algebra).

Let us now treat the given inequality:

solve(x**2-78>=3, x)
[[x <= -9], [x >= 9]]

As usual, we see that the solution set of certain inequalities consists of the union and intersection of open intervals.

The next is again a polynomial equation of degree 3, so we expect 3 roots over the complex numbers:

solve(x^3-1==0, x)
[x == 1/2*I*sqrt(3) - 1/2, x == -1/2*I*sqrt(3) - 1/2, x == 1]

Note however that the choice of the ring affects the result!

x = PolynomialRing(RationalField(), 'x').gen() #this creates a polynomial ring over the rational numbers (ℚ) with the variable x. f = x^3 - 1 f.roots()
[(1, 1)]

By creating the polynomial ring over the rational numbers (ℚ) using the command PolynomialRing(RationalField(),x){\tt{PolynomialRing(RationalField(), 'x')}}, you are working within a context where only rational numbers (fractions) are considered. So in this case the roots(){\tt{roots()}} method only provides roots that are rational numbers!

The command .gen(){\tt{.gen()}} used above returns the generator of this polynomial ring, which in this case is the variable xx. The comand f.roots(){\tt{f.roots()}} computes the roots of the polynomial f(x)=x31f(x)=x^3-1. It returns a list of pairs, where each pair consists of a root and its multiplicity.

solve(e^x+1==0, x)
[x == I*pi]
solve(exp(x)-x^2 == 0 , x)
[x == -e^(1/2*x), x == e^(1/2*x)]

As we pronounced , to find a numeric approximation we can use the command find_root(){\tt{find\_root()}}. This requires also a closed interval on which to search for a solution (however, this command will only return one solution on the specified interval, if one exists. It will not find the complete solution set over the entire real numbers)

find_root(exp(x) -x^2 == 0, -1 , 1)
-0.7034674224983917

Let us illustrate the solution (and see below for more examples in 2D-plotting )

plot(exp(x) -x^2, x, -1, 1)
Image in a Jupyter notebook
find_root(sin(x) == x, -pi/2 , pi/2)
0.0

Let us also solve the equation given in 8. In this point recall that log{\tt{log}} in Sage means the natural logarithm.

solve(log(x^2)==6/5, x)
[x == -e^(3/5), x == e^(3/5)]

Hence we get two solutions satisfying the given equation.

Remark

We mention that the command solve(){\tt{solve()}} can be very slow for large systems of equations. Below we present a few simple examples.

Notice that systems of linear equations will be treated later in terms of linear algebra functions.

reset('y')
solve( [3*x -2* y == 2, -2*x -5*y == 1 ], x,y) #this will give us an error
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Input In [10], in <cell line: 1>() ----> 1 solve( [Integer(3)*x -Integer(2)* y == Integer(2), -Integer(2)*x -Integer(5)*y == Integer(1) ], x,y) NameError: name 'y' is not defined

This mean that we forgot to introduce yy as a variable, since in fact we first typed reset(y){\tt{reset(y)}}, in order to ``clear'' some possible assignement to this variable.

This is because SageMath has been programmed to keep in its memory any assignment we have done, untill we decide to reset it.

Thus, to solve the system of equations given above we may proceed as follows:

y=var("y") solve( [3*x -2* y == 2, -2*x -5*y == 1 ], x,y)
[[x == (8/19), y == (-7/19)]]

Another example:

solve([2*x + y == -1 , -4*x - 2*y == 2],x,y)
[[x == -1/2*r1 - 1/2, y == r1]]

Above, r1r1 (or better written r1r_1) is a free variable (and this parametrizes the solution set).

When there is more than one free variable, SageMath enumerates them by r1,r2,,rkr1,r2, \ldots, rk (that is, by r1,,rkr_1, \ldots, r_k).

x, y, z=var("x, y, z") solve([ 2*x + 3*y + 5*z == 1, 4*x + 6*y + 10*z == 2, 6*x + 9*y + 15*z == 3], x,y,z)
[[x == -5/2*r2 - 3/2*r3 + 1/2, y == r3, z == r2]]

Finally, let us solve a system of inequalities in several variables (this can lead to complicated expressions).

solve([ x-y >=5, 2*x+y <=3], x,y)
[[x == (8/3), y == (-7/3)], [x == -1/2*y + 3/2, y < (-7/3)], [x == y + 5, y < (-7/3)], [y + 5 < x, x < -1/2*y + 3/2, y < (-7/3)]]

To explore a full list of options of the command solve{\tt{solve}} in SageMath, type the following:

solve?

2.2. One-variable functions and 2D-plotting

In SageMath the graph of a real-valued function ff defined in some subset II of the real line R\R, is given by the command plot(f(x),x,a,b){\tt{plot(f(x), x, a, b)}}, where a,ba, b are the end points of II.

Exercise

Present the graph of the folowing trigonometric functions for x[2π,2π]x\in[-2\pi, 2\pi]:

  1. f(x)=sin(x)f(x)=\sin(x),

  2. μ(x)=sin(2x)\mu(x)=\sin(2x),

  3. g(x)=cos(x)g(x)=\cos(x),

  4. ϕ(x)=tan(x)\phi(x)=\tan(x).

Solution:
  1. For the function sin(x)\sin(x), we get the graph

f(x)=sin(x) plot(f(x), x, -2*pi, 2*pi)
Image in a Jupyter notebook

Or, after introducing the function f(x)f(x) we could just type f.plot(2pi,2pi){\tt{f.plot(-2*pi, 2*pi)}}.

  1. Similarly, in this case we type

plot(sin(2*x), x, -2*pi, 2*pi)
Image in a Jupyter notebook
  1. Let us now plot the function cos(x)\cos(x):

g(x)=cos(x) plot(g(x), x, -2*pi, 2*pi)
Image in a Jupyter notebook
  1. The function ϕ(x)=tan(x)\phi(x)=\tan(x) has asymptotes. In SageMath to produce a meaningful plot for this function one can proceed with the syntax

fi(x)=tan(x) plot(fi(x), x, -2*pi, 2*pi, detect_poles='show', ymin=-50,ymax=50)
Image in a Jupyter notebook

Remark

The options xmin (or ymin) and xmax (or ymax) are interval bounds over the xx-axis (respectively over yy-axis) for which the Sage will display the function.

For instance, above we required the function tan(x)\tan(x) to be dispalyed in the yy-axis from 50-50 to 5050.

So this option changes the window range of the graph.

We also used the option detect_poles{\tt{detect\_poles}} since its defaul value is False.

This option often enables to draw a vertical asymptote at poles of the function.

Note that without these options the graph of tan(x)\tan(x) will not appear correctly.

Other options that we can often use are color, thickness, linestyle etc.

Thickness is about the thickness of the curves in the graph, and the defaul value is 1.

Linestyle has as default the solid, or "-", while we also have further options for this, as “--” or “dashed” , “.-.” or “dashdot”, “” or “dotted”, etc.

Later we will see that we can add text in a graph, plot a single point, give labels to the axes, etc.

Exercise

Plot in one figure the graph of f(x)=x24f(x)=\frac{x^2}{4} for 8x8-8\leq x\leq 8 and the line passing through the points [4,4][4, 4] and [2,2][2, -2], together with these two points.

Solution:

f(x) = x^2 / 4 a=plot(f, -8, 8, rgbcolor=(0.2,0.2,0.4)) b=line([(4, 4), (2, -2)]) c=point((4,4),pointsize=30) d=point((2,-2),pointsize=30) (a+b+c+d).show(ymin=-3, ymax=12, aspect_ratio=1)
Image in a Jupyter notebook

Here we fixed the option aspect_ratio=1{\tt{aspect\_ratio=1}}, to have equal scales for xx and yy.

More details on points, lines, etc, we will meet below when we will use Sage to study plane geometry.

Exercise

Plot the following functions for the given domain:

  1. f(x)=exp(x)ln(x)f(x)=\displaystyle\frac{\exp(x)}{\ln(x)}, x[0,4]x\in[0, 4];

  2. g(x)=exp(2x)xg(x)=\displaystyle\frac{\exp(2x)}{x}, x[2,4]x\in[2, 4];

  3. h(x)=xsin1xh(x)=\displaystyle x\sin\frac{1}{x}, for x(1,1)x\in(-1, 1).

  4. r(x)=2(x1)(x2)(x4)r(x)=\displaystyle \frac{2}{(x-1)(x-2)(x-4)}, for x[2π,2π]x\in[-2\pi, 2\pi].

Solution:

  1. Recall that the domain of ln(x)\ln(x) is the open set (0,+)(0, +\infty), with ln(1)=0\ln(1)=0, hence the given function ff is not defined at x=1x=1. Sage indicate this as follows:

ff(x)=exp(x)/ln(x) plot(ff(x), x, 0, 4, detect_poles='show', ymin=-100,ymax=100, rgbcolor=(0.2,0.2,0.4))
Image in a Jupyter notebook
  1. In the second case the given function is well-defined since for instance is required x[2,4]x\in[2, 4] and so x0x\neq 0.

g(x)=exp(2*x)/x plot(g(x), x, 2, 4, color="black", linestyle='dotted') # Note that in this case we asked Sage for black color and linestyle ``dotted''
Image in a Jupyter notebook
  1. In this case one arrives to the graph

h(x)=x*sin(1/x) plot(h(x), x, -1, 1, color="darkgreen")
Image in a Jupyter notebook
  1. In this case one can use the options detect_poles="True"{\tt{detect\_poles="True"}}, or detect_poles="show"{\tt{detect\_poles="show"}}, together with an arrangement of window size options. In particular, we may proceed with the syntax

r(x)=2/((x-1)*(x-2)*(x-4)) plot(r(x), x, -2*pi, 2*pi, detect_poles='show').show(ymin=-10,ymax=10)
Image in a Jupyter notebook

Or we could directly type: plot(r(x),x,2pi,2pi,detect_poles="show",ymin=10,ymax=10){\tt{plot(r(x), x, -2*pi, 2*pi, detect\_poles="show", ymin=-10,ymax=10)}}.

Tips

{\bullet} In SageMath it is very easy to join two (or more) graphs in one system of axes together, as in the following example:

a=plot(exp(x)/ln(x), x, 0, 4, detect_poles='show', ymin=-100,ymax=100, rgbcolor=(0.2,0.2,0.4)) b=plot(exp(2*x)/x, x, 0, 4, color="darkcyan", linestyle='dotted', thickness=2) show(a+b)
Image in a Jupyter notebook

Let us present another similar example:

P = plot(sin(x)/x, -10, 10, color='blue') + \ plot(x*cos(x), -10, 10, color='red') + \ plot(1/tan(x), -10, 10, color='green') P.show(ymin=-pi, ymax=pi)
Image in a Jupyter notebook

\bullet When one is not interested in the domain of a given function ff, it is still possible to plot ff as follows (Sage will consider itself a domain for ff).

f(x)=3*x+4 f.plot()
Image in a Jupyter notebook

{\bullet} In SageMath one can use the option plot_points{\tt{plot\_points}}, which is about the minimal number of computed points, with 200 as default value.

For instance:

plot(ln(x), (x,0,10), plot_points=20, linestyle='', marker='.')
Image in a Jupyter notebook
plot(ln(x), (x,0,10), plot_points=60, linestyle='', marker='.')
Image in a Jupyter notebook

Note that the marker can be a TeX symbol. As an example

plot(ln(1/x), (x, 0, 10), plot_points=20, linestyle='', marker=r'$\flat$')
Image in a Jupyter notebook

\bullet We can add a name/title to a plot (see also below for the command legend_label{\tt{legend\_label}}). Let us give an example.

plot(3*x^2+2*x+1, (x,-4,4), color="darkred", title='A plot of $3x^2+2x+1}$')
Image in a Jupyter notebook

\bullet In fact, we can also set the position of the title, including coordinates, via the option title_pos=(x,y){\tt{title\_pos=(x,y)}}. For instance:

plot(3*x^2+2*x+1, (x,-4,4), title='A plot of $3x^2+2x+1$', title_pos=(0.7,0.5))
Image in a Jupyter notebook

\bullet We can add labels to the axes:

plot(sqrt(3)*x^3/3, (x,-1,1), axes_labels=['$x$', '$y$'], color="brown")
Image in a Jupyter notebook

\bullet. In some cases we may want the axes to use values such as π/3,2π/3\pi/3, 2\pi/3, etc.

We may obtain this by using the option ticks=π/3,tick_formatter=π{\tt{ticks=\pi/3, tick\_formatter=\pi}}.

Let us present an example:

F=plot(cos(x), x, -2*pi, 2*pi, color="darkorange", ticks = pi/3, tick_formatter=pi) F.show()
Image in a Jupyter notebook

\bullet We may also like to add a legend. For instance

plot(2*x+1, x, -2, 4, color="darkblue", legend_label="the line $y=2x+1$")
Image in a Jupyter notebook

\bullet. Sage allows us to fix the color that we like in the legend, as below:

plot(sin(x)/x, x, -10, 10, color="darkblue", legend_label="the plot of $sin(x)/x$", legend_color="blue")
Image in a Jupyter notebook

\bullet Often it can be useful to plot together two graphs (or more) but using distinct system of axes, or even labeling. For instance

gr1=plot(cos(x), x, -3*pi, 3*pi, color="purple", ticks = pi/3, tick_formatter=pi) gr2=plot(cos(4*x), x, -3*pi, 3*pi, color="orange") show(graphics_array([gr1, gr2], 2, 1))
Image in a Jupyter notebook

Note that if we write the last command as show(graphics_array([gr1,gr2],1,2)){\tt{show(graphics\_array([gr1, gr2], 1,2))}} we will not obtain an optimal result:

show(graphics_array([gr1, gr2], 1, 2)) # in the command graphics_array the final option 1, 2 or 2, 1, #or more general a, b, splits the figure into a ``lines'' and b ``columns''.
Image in a Jupyter notebook

To explore further options of the command plot{\tt{plot}} in SageMath, type the following and scroll down the appearing file.

plot?

Exercise for practice

Plot together the graphical presentation of the functions f(x)=cos(x)f(x)=\cos(x) and g(x)=cos(2x)g(x)=\cos(2x) for x[4π,4π]x\in[-4\pi, 4\pi], label the axes, and for the xx-axis use values as π/3,2π/3\pi/3, 2\pi/3, etc.

Finally put a title to the figure.

Plot of piecewise functions

Suppose that f(x)f(x) is a real-valued function defined by f(x)={f1,if x(a,b),f2,if x[b,d). f(x)= \left\{ \begin{matrix} f_1, & if \ x\in(a, b), \\ f_2, & if \ x\in[b, d). \end{matrix}\right.

or with different kind of domains. In Sage to define such functions we can use the piecewise{\tt{piecewise}} command. Let us present examples.

F = piecewise([((0,10), x^3), ([-10,0], -x^2)]); F
piecewise(x|-->x^3 on (0, 10), x|-->-x^2 on [-10, 0]; x)

To obtain the collection of domains of the component functions we may type F.domains(){\tt{F.domains()}}. On the other hand, to otain the end points of the intervals, we may type F.end_points(){\tt{F.end\_points()}}. For our example this gives

F.domains()
((0, 10), [-10, 0])
F.end_points()
[-10, 0, 10]

Exercise

Plot the following functions in SageMath:

f(x)={1,if x(π,0),1,if x[0,π).f(x)= \left\{ \begin{matrix} -1, & if \ x\in(-\pi, 0), \\ 1, & if \ x\in[0, \pi). \end{matrix}\right.
g(x)={x3,if x(π,0),x2,if x(0,π).g(x)= \left\{ \begin{matrix} -x^3, & if \ x\in(-\pi, 0), \\ x^2, & if \ x\in(0, \pi). \end{matrix}\right.

Solution:

For the first case we use the cell

f2(x) = 1; f1(x) = -1 f = piecewise([[(-pi,0),f1],[(0,pi),f2]]) a=plot(f(x), x, -pi, pi) b=point((0, 1), size = 50) show(a+b)
Image in a Jupyter notebook

To correct the above graph (that is, to avoid the vertical line, which is not part of the graph), we should use the exclude{\tt{exclude}} command.

f2(x) = 1; f1(x) = -1 f = piecewise([[(-pi,0),f1],[(0,pi),f2]]) a=plot(f(x), x, -pi, pi, exclude=[0]) b=point((0, 1), size = 50) show(a+b)
Image in a Jupyter notebook

For the second case, we may type:

G = piecewise([[(-pi,0),-x^3],[(0,pi), 4*x]]) plot(G(x), x, -pi, pi, color="red")
Image in a Jupyter notebook

Another Example of 2D-graphics (lines Vs circles)

a=sum(line([(0,0),(cos(i),sin(i))], color="darkcyan", linestyle="dotted") for i in [0,pi/20,..,2*pi]) b=sum(circle((0,0), i, color = 'darkgreen', thickness = 1.5) for i in [0,0.05,..,1]) c=point((0, 0), size=30, color = 'red') show(a+b+c)
Image in a Jupyter notebook

Note that editing where ii is evaluated, we result to different figures:

a=sum(line([(0,0),(cos(i),sin(i))], color="darkcyan", linestyle="dotted") for i in [0,pi/90,..,2*pi]) b=sum(circle((0,0), i, color = 'darkgreen', thickness = 1.5) for i in [0,0.3,..,1]) c=point((0, 0), size=30, color = 'red') show(a+b+c)
Image in a Jupyter notebook

A colourful example (rainbow colors)

Ep = plot([]) # Define empty plot mcols = rainbow(20) # Define 20 evenly spaced colors for k in range(20): Ep = Ep + plot(sin(x+k)/(x+k),1,10, color=mcols[k],thickness=1.5) #one may fix another function Ep # Print graphs
Image in a Jupyter notebook

2.3. Basics on sets and combinatorics

Combinatorics is a branch of mathematics that deals with counting and arranging objects in a systematic way.

Combinatorics is used in a wide range of fields, including computer science, statistics, physics, and biology.

It is also used in various real-life situations, such as in the design of codes and ciphers, scheduling of events, and the arrangement of seating in a theater or stadium.

With Sage it's much easier to do combinatorics in a quick way, as we will see soon via examples.

Sets

We will be working with sets a lot, so let's see how do declare them and what you can do with them after.

set_4_values = Set([1,2,3,4]) # create the set from list of integers set_4_values_str = Set(['One','Two','Three','Four']) # also works with list of strings print('Set of integers:', set_4_values, '\nSet of strings:', set_4_values_str)
Set of integers: {1, 2, 3, 4} Set of strings: {'Four', 'Three', 'One', 'Two'}

Sets can be finite and infinite and we can check this with .is_finite() method

print('If Rational field set finite?:', Set(QQ).is_finite()) print('If set of interegers finite?:', set_4_values.is_finite())
If Rational field set finite?: False If set of interegers finite?: True

Basic operations with sets

We now describe some basic operations that we can apply when we use sets.

A) Intersection

The intersection of some given sets can be treated with .intersection() method. For example:

X = Set(ZZ).intersection(Primes()) # intersection between Integers Field and Primes Field print('Is 4 in set X?:', 4 in X) print('Is 3 in set X?:', 3 in X) print('Is 101 in set X?:', 101 in X) print('Is 6/2 in set X?:', 6/2 in X) print('Is 4/3 in set X?:', 4/3 in X)
Is 4 in set X?: False Is 3 in set X?: True Is 101 in set X?: True Is 6/2 in set X?: True Is 4/3 in set X?: False
B) Union

The union of some given sets can be treated with the .union() method or plus '+' sign. For instance:

X = Set(ZZ).union(Set([math.e, 4/3])); # intersection between Rational Field and Primes Field print('Is 4 in set X?:', 4 in X) print('Is e=2.7182... in set X?:', math.e in X) print('Is 101 in set X?:', 101 in X) print('Is 5/2 in set X?:', 5/2 in X) print('Is 4/3 in set X?:', 4/3 in X) Y = Set([1,2,3,4]) + Set([10,9]) # you can use plus '+' sign instead and perform with any set print('\nIs 4 in set Y?:', 4 in Y) print('Is 5 in set Y?:', 9 in Y) print('Is 9 in set Y?:', 9 in Y) print('Is 10 in set Y?:', 10 in Y)
Is 4 in set X?: True Is e=2.7182... in set X?: True Is 101 in set X?: True Is 5/2 in set X?: False Is 4/3 in set X?: True Is 4 in set Y?: True Is 5 in set Y?: True Is 9 in set Y?: True Is 10 in set Y?: True
C) Difference of sets

The difference of some given sets is treated with the .difference() method or minus '-' sign

X = Set(ZZ).difference(Primes()) # difference between Integers Field and Primes Field print('Is 4 in set X?:', 4 in X) print('Is 3 in set X?:', 3 in X) print('Is 101 in set X?:', 101 in X) print('Is 4/1 in set X?:', 4/1 in X) print('Is 4/3 in set X?:', 4/3 in X) Y = Set([1,2,3,4]) - [3,1] # you can use minus '-' sign instead and perform with any iterable print('\nIs 4 in set Y?:', 4 in Y) print('Is 3 in set Y?:', 3 in Y) print('Is 1 in set Y?:', 1 in Y)
Is 4 in set X?: True Is 3 in set X?: False Is 101 in set X?: False Is 4/1 in set X?: True Is 4/3 in set X?: False Is 4 in set Y?: True Is 3 in set Y?: False Is 1 in set Y?: False
D) Cartesian product of sets

We can also do Cartesian product with different sets with cartesian_product() function:

Hint: Cartesian_prod=Set1×Set2={(s1,s2)s1Set1,s2Set2}{\tt{Cartesian\_prod = Set1 × Set2 = \{(s1, s2) | s1 ∈ Set1, s2 ∈ Set2\}}}

fruits = Set(['Apple', 'Banana', 'Pineapple', 'Pear']) vegetables = Set(['Tomato', 'Corn', 'Letuse', 'Cucumber']) fruit_and_veg = cartesian_product([fruits, vegetables]) print(fruit_and_veg)
The Cartesian product of ({'Banana', 'Apple', 'Pineapple', 'Pear'}, {'Corn', 'Tomato', 'Cucumber', 'Letuse'})

We can check the cardinality of our sets and Cart. products with .cardinality() method

print('Size of set with vegetables:', vegetables.cardinality()) print('Size of cartesian product with fruits & vegetables:', fruit_and_veg.cardinality())
Size of set with vegetables: 4 Size of cartesian product with fruits & vegetables: 16

We can choose random elements from our sets with .random_element() method

print("\nRandom fruit:", fruits.random_element()) print("\nRandom fruit&veg:", fruit_and_veg.random_element())
Random fruit: Banana Random fruit&veg: ('Banana', 'Tomato')

Subsets

To create subsets from some give set, we can apply yhe Subsets() wrapper, as below.

Hint: If AA and BB are sets, then AA is a subset of BB, denoted by ABA \subseteq B, if and only if every element of AA is also an element of BB

AB(x)(xAxB)A \subseteq B \quad \Leftrightarrow \quad (\forall x) \left(x \in A \Rightarrow x \in B\right)

Breakfast = Subsets(fruit_and_veg, 3) # specify the set and number of elements in a desired subset print(Breakfast) print('Random breakfast consists of:', Breakfast.random_element())
Subsets of The Cartesian product of ({'Banana', 'Apple', 'Pineapple', 'Pear'}, {'Corn', 'Tomato', 'Cucumber', 'Letuse'}) of size 3 Random breakfast consists of: {('Pineapple', 'Letuse'), ('Apple', 'Tomato'), ('Pear', 'Letuse')}

Combinatorics

In this section we use natural numbers to describe some indivisible items located in real life space, and deal with questions in how many ways certain events may appear.

For instance, how to compute the number of their (pre)orderings, specific choices, and so on. In many of these problems, “common sense” is sufficient.

We just need to use the rules of product and sum in the right way, as we show in the following examples. We proceed with exercises on combinatorics.

Permutations (without repetition)

Suppose we have a set of nn (distinguishable) objects, and we wish to arrange them in some order. We can choose a first object in nn ways, then a second in n1n − 1 ways, a third in n2n − 2 ways, and so on, until we choose the last object for which there is only one choice. The total number of possible arrangements is the product of these, hence there are exactly n!=n(n1)(n2)...321n! = n\cdot(n − 1)\cdot(n − 2)\cdot. . .\cdot 3 \cdot 2 \cdot 1 distinct orders of the objects.

The number P(n)P(n) of distinct orderings of a finite set with nn elements is given by the factorial function: p(n)=n!p(n) = n!

Exercise

Find out how many permutations we can make with integers from 11 to 33 and from 11 to 44

Solution:

We can use built-in Permutations(num_of_int){\tt{Permutations(num\_of\_int)}} function for 3 and then for 4 permutations

perm_3 = Permutations(3) print('All permutations from 1 to 3:', perm_3.list(), '\nNumber of permutations:', perm_3.cardinality()) perm_4 = Permutations(4) print('\nAll permutations from 1 to 4:', perm_4.list(), '\nNumber of permutations:', perm_4.cardinality())
All permutations from 1 to 3: [[1, 2, 3], [1, 3, 2], [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1]] Number of permutations: 6 All permutations from 1 to 4: [[1, 2, 3, 4], [1, 2, 4, 3], [1, 3, 2, 4], [1, 3, 4, 2], [1, 4, 2, 3], [1, 4, 3, 2], [2, 1, 3, 4], [2, 1, 4, 3], [2, 3, 1, 4], [2, 3, 4, 1], [2, 4, 1, 3], [2, 4, 3, 1], [3, 1, 2, 4], [3, 1, 4, 2], [3, 2, 1, 4], [3, 2, 4, 1], [3, 4, 1, 2], [3, 4, 2, 1], [4, 1, 2, 3], [4, 1, 3, 2], [4, 2, 1, 3], [4, 2, 3, 1], [4, 3, 1, 2], [4, 3, 2, 1]] Number of permutations: 24

We can also use the factorial(n){\tt{factorial(n)}} function

Hint: factorial(n)=n!=n(n1)(n2)...321factorial(n) = n! =n\cdot(n − 1)\cdot(n − 2)\cdot. . .\cdot 3 \cdot 2 \cdot 1

num_of_3_perm = factorial(3) num_of_4_perm = factorial(4) print('Number of permutations for 3 values:', num_of_3_perm) print('Number of permutations for 4 values:', num_of_4_perm)
Number of permutations for 3 values: 6 Number of permutations for 4 values: 24

Exercise

In how many ways can 44 chemistry, 33 math, and 22 physics books be set on a shelf, if the chemistry books are set on the left, math books in the middle, and physics books on the right?

Solution:

In this case we will use the factorial(n){\tt{factorial(n)}} function for each group and then multiply them.

num_of_chemistry_perm = factorial(4) # number of ways to permutate chemistry books num_of_math_perm = factorial(3) # number of ways to permutate math books num_of_physics_perm = factorial(2) # number of ways to permutate physics books print('Number of permutations for chemistry books:', num_of_chemistry_perm) print('Number of permutations for math books:', num_of_math_perm) print('Number of permutations for physics books:', num_of_physics_perm) # As far as our book subjects are predefined (in the meaning of positions on the shelf) # we just need to multiply our permutations overall_num = num_of_chemistry_perm * num_of_math_perm * num_of_math_perm print('Number of all permutations:', overall_num)
Number of permutations for chemistry books: 24 Number of permutations for math books: 6 Number of permutations for physics books: 2 Number of all permutations: 864

Exercise

In how many ways can 44 chemistry, 33 math, and 22 physics books be set on a shelf, if they are grouped by subject?

Solution:

We will, again, use the factorial(n){\tt{factorial(n)}} function for each subject and then add them up, but as we can change the place of subjects on shelf, we will also multiply the sum by the number of subjects.

num_of_chemistry_perm = factorial(4) # number of ways to permutate chemistry books num_of_math_perm = factorial(3) # number of ways to permutate math books num_of_physics_perm = factorial(2) # number of ways to permutate physics books print('Number of permutations for chemistry books:', num_of_chemistry_perm) print('Number of permutations for math books:', num_of_math_perm) print('Number of permutations for physics books:', num_of_physics_perm) # Suppose that our book subjects are predefined (in the meaning of positions on the shelf), # in this case we just multiply our permutations if_fixed_overall_num = num_of_chemistry_perm * num_of_math_perm * num_of_math_perm print('Number of permutations if subjects are fixed:', if_fixed_overall_num) # But in our case, subjects can be placed in any order, # as far as we have 3 subjects, number of permutations is 3!, # so we need to multiply by 3! overall_num = factorial(3) * if_fixed_overall_num print('Number of all permutations:', overall_num)
Number of permutations for chemistry books: 24 Number of permutations for math books: 6 Number of permutations for physics books: 2 Number of permutations if subjects are fixed: 864 Number of all permutations: 5184

kk-Permutations (without repetitions)

Suppose SS is a set with nn elements. Suppose we wish to choose and arrange in order just kk of the members of SS, where 1kn1 ≤ k ≤ n.

This is called a k-permutation without repetition of the nn elements.

The same reasoning as above shows that this can be done in: vnk=n!(nk)!v_n^k=\frac{n!}{(n-k)!} ways.

The right side of this result also makes sense for k=0k = 0, (there is just one way of choosing nothing), and for k=nk = n, since 0!=10! = 1.

Exercise

How many 4-letter words(meaning a set where order of elements plays role) we can make with ABCDEFGABCDEFG where letters can't repeat?

Solution:

We need to find kk-letter words, which means that we can include only kk(4 in our case) elements from our group of nn(7 in our case) without repetitions, where the order counts, for that we can use this formula vnk=n!(nk)!v_n^k=\frac{n!}{(n-k)!}Where,

nn - is the number of all elements,

kk - is the number of allowed elements from all

In Sage for this we can define a function, as follows:

def u_perm (n, k) : a = factorial(n) / factorial(n-k) return a u_perm(7, 4)
840

Let's see how we can double-check our answer with the practical experiment.

letters = Set(['A','B','C','D','E','F','G']) # our letters set letters_subset = Subsets(# create a subsets letters, # from "letters" set k=4, # with 4 elements, which represent our 4-letter words submultiset=False # where elements don't repeat ) for _ in range(4): print(letters_subset.random_element()) # let's see an examples of such words # We received our subsets, but as we know, set doesn't count for order of elements, so we should account for it # So we have 35 subsets print('Number of sets (not including order of elements):', letters_subset.cardinality()) # And each of them has 4 elements, and, as we alredy know, the formula to permutate 4 elements is factorial(4)=4! # So, to permutate all 35 subsets with 4 elements each we need to multiply it by factorial of 4! print('Number of sets (including the order)', letters_subset.cardinality()*factorial(4))
{'A', 'D', 'C', 'B'} {'E', 'F', 'D', 'C'} {'F', 'E', 'D', 'A'} {'A', 'C', 'F', 'G'} Number of sets (not including order of elements): 35 Number of sets (including the order) 840

Exercise

We have 55 physics books and 66 math books to put on a shelf with 6 slots.

In how many ways can we put the books on the shelf, if the first two slots needs to be filled with physics books, and the next four with math books (order plays role)?

Solution:

We can include only 22 physics books from 55, and 44 math books from 6, the order counts, so we can use UnkU_n^k formula for both subjects and add the results

def v_perm (n, k) : a = factorial(n) / factorial(n-k) return a num_of_math_perm = v_perm(6,4) # number of ways to find permutations of 6 math books and take 4 num_of_physics_perm = v_perm(4,2) # number of ways to find permutations of 4 physics books and take 2 print('Number of permutations of 6 math books and take 4:', num_of_math_perm) print('Number of permutations of 4 physics books and take 2:', num_of_physics_perm) # Suppose that our book subjects are predefined (in the meaning of positions on the shelf), # in this case we just add our permutations # As far as our book subjects are predefined (in the meaning of slots in the shelf) # we just need to multiply our permutations overall_num = num_of_math_perm * num_of_math_perm print('Number of all permutations:', overall_num)
Number of ways to find permutations of 6 math books and take 4: 360 Number of ways to find permutations of 4 physics books and take 2: 12 Number of all permutations: 129600

Exercise

We have 55 physics books and 66 math books to put on a shelf with 44 slots. In how many ways can we put the books on the shelf (order plays role)?

Solution:

We can only put 4 books into shelf from 11, the order counts, so we can use UnkU_n^k formula for both subjects and add the results

def v_perm (n, k) : a = factorial(n) / factorial(n-k) return a num_of_math_perm = v_perm(11,4) # number of ways to find permutations of 6 math books and take 4 print('Number of ways to find permutations of 6 math books and 5:', num_of_math_perm)
Number of ways to find permutations of 6 math books and take 4: 7920

kk-Combinations (without repetitions)

Consider a set SS with n elements. A kk-combination of the elements of SS is a selection of kk elements of SS, 0kn0 ≤ k ≤ n, when order does not matter.

For k1k ≥ 1, the number of possible results of a subsequential choosing of our kk elements, is n(n1)(n2)...(nk+1)n(n − 1)\cdot(n − 2)\cdot ... \cdot(n − k + 1) (a kk-permutation).

We obtain the same kk-tuple in k!k! distinct orders. Hence the number of kk-combinations is (nk)=vnkk!=n!(nk)!k!{n\choose k}=\frac{v_{n}^{k}}{k!}=\frac{n!}{(n-k)!\cdot k!} If k=0k = 0, the same formula is still true, since 0!=10! = 1, and there is just one way to select all nn elements.

Exercise

How many 4-letter sets(meaning a regular set, where order of elements doesn't play role) we can make with ABCDEFG, where letters can't repeat?

Solution:

We need to find kk-letter sets, that means that we can include only kk(4 in our case) elements from our group of nn(7 in our case) without repetitions, where the order of elements doesn't count, for that we can use this formula (nk)=n!(nk)!k!{n\choose k}=\frac{n!}{(n-k)!\cdot k!}Where,

nn - is the number of all elements,

kk - is the number of allowed elements from all

Basically, what we've changed is that we divided the order our (represented by k!) from our kk-permutations formula

def combinations(n, k): a = factorial(n) / (factorial(n-k) * factorial(k)) return a combinations(7, 4)
35

Or we can use already defined function in SAGE called binomial(n,k)binomial(n,k)

binomial(7,4)
35

We can prove this with sets

letters = Set(['A','B','C','D','E','F','G']) # our letters set letters_subset = Subsets(# create a subset letters, # from "letters" set k=4, # with 5 elements, which represents our 5-letter set submultiset=False # where elements don't repeat ) for _ in range(4): print(letters_subset.random_element()) # let's see an examples of such a words print( '\nNumber of sets fitting our definition:', letters_subset.cardinality() ) #letters_subset.list() # uncomment the line to see all sets
{'A', 'D', 'C', 'B'} {'E', 'F', 'C', 'B'} {'A', 'E', 'C', 'D'} {'F', 'D', 'B', 'A'} Number of sets fitting our definition: 35

Exercise

If a coin is tossed seven times, in how many ways can it fall 55 heads and 22 tails?

Solution:

Suppose we have 7 spots to put the coins on.

If we choose any 22 spots for tails, the other 55 will automatically be tails. So we can use (nk){n\choose k} formula, where from 7 we choose 2:

def combinations(n, k): a = factorial(n) / (factorial(n-k) * factorial(k)) return a print('Number of combinations:', combinations(7, 2)) # 'or' binomial(7,2)
Number of combinations: 21 or 21

But we can also go from the other way, where from 7 spots we choose 5 spots for heads and 2 tails would be set automatically:

print('Number of combinations:', combinations(7, 5))
Number of combinations: 21

Exercise

The government consists of 5050 people: 1515 from green party, 77 from red party, and 2828 from blue party.

In how many ways could a 1212 person committee be selected if it contain's 22 red, 66 blue, and 44 green?

Solution:

We need to get 2 people from 15 red, 4 people from 30 green and 6 from 55 blue.

So we can use (nk){n\choose k} formula, for each party, and then multiply the results, as all three combinations are independent from each other and can be mixed

def combinations(n, k): a = factorial(n) / (factorial(n-k) * factorial(k)) return a comb_of_red = combinations(7, 2) comb_of_green = combinations(15, 4) comb_of_blue = combinations(28, 6) print('Number of red combinations:', combinations(7, 2)) print('Number of green combinations:', combinations(15, 4)) print('Number of blue combinations:', combinations(28, 6)) all_comb_result = comb_of_red*comb_of_green*comb_of_blue print('Possible commities from 12 people, where 4 green, 2 red and 6 blue:', all_comb_result)
Number of red combinations: 21 Number of green combinations: 1365 Number of blue combinations: 376740 Possible commities from 12 people, where 4 green, 2 red and 6 blue: 10799252100

kk-Permutations (with repeatitions)

Let SS be a set with nn distinct elements. We wish to select kk elements, 0kn0 ≤ k ≤ n from SS with repetition permitted.

This is called a kk-permutation with repetition.

Since the first selection can be done in nn ways, and similarly the second can also be done in nn ways etc.

The total number V(n,k)V(n, k) of kk-permutations with repetitions is nkn^k. Hence V(n,k)=nkV(n,k)=n^k

Exercise

Consider that we have an car plate, which consists of 55 digits. How many plates can be created?

Solution:

When we're dealing with the repeated permutations we can use the following formula V(n,k)=nkV(n,k)=n^kWhere,

nn - number of elements,

kk - number of wanted elements

We have 1010 digits (number of elements nn) and we want to create car plates with 55 digits (number of wanted elements kk)

print('Number of possible plates:', 10**5)
Number of possible plates: 100000

Exercise

Imagine you have a coin and you're throwing it 1010 times. What is the number of possible outcomes?

Solution:

Applying our formula V(n,k)V(n,k), where 2 is number of elements (n=2n=2) and 10 wanted elements (k=10k=10)

print('Number of possible outcomes:', 2**10)
Number of possible outcomes: 1024

kk-Combinations (with repeatitions)

If we are interested in a choice of kk elements without taking care of order, we speak of kk-combinations with repetitions.

At first sight, it does not seem to be easy to determine the number.

We reduce the problem to another problem we have already solved, namely combinations without repetitions:

The number of k-combinations with repetitions from nn elements equals for every k0k ≥ 0 and n1n ≥ 1 (nk)repetitions=(n+k1k){n\choose k}_{\rm repetitions}={n+k-1\choose k}

Exercise

From 55 types of beer (consider infinite amount of each), one of which is dark, you are putting 1010 of them on the table.

Determine the number of combination of beers you can put on the table.

Solution:

When we're dealing with the repeated permutations, we can use the following formula (nk)repetitions=(n+k1k){n\choose k}_{\rm repetitions}={n+k-1\choose k} Where,

nn - number of types of elements,

kk - number of wanted elements.

Applying to our task we have 5 types of beer (5 types of elements) and 10 bottles on the table (10 wanted elements)

def comb_repet(n, r): return binomial(n+r-1, r) print('Number of possible combinations of beers:', comb_repet(5, 10))
Number of possible combinations of beers: 1001

Exercise (continuation)

From 55 types of beer (consider infinite amount each), one of which is dark, you are putting 1010 of them on the table.

Determine the number of combination of beers if you must include at least one dark beer.

Solution:

As we have already selected one dark beer, we can still put 9 more beers of 5 types.

We can use our (nr)repetitions{n\choose r}_{\rm repetitions} formula, where our n=5n=5 and r=9r=9

def comb_repet(n, r): return binomial(n+r-1, r) print('Number of combinations of beer on the table with at least 1 dark beer:', comb_repet(5, 9))
Number of combinations of beer on the table with at least 1 dark beer: 715
f(x) = (exp(2/x)+1-(3/log(x)))/(sin(x+1/x)) f.limit(x=0, dir='+')
x |--> 0
quest = 140 right_ans = 43 prob = 1/4 def count (right_ans): return binomial(140, right_ans) * prob**right_ans * (1-prob)**(140-right_ans) result = 0 for r_ans in range(43, 140+1): result = result + count(r_ans) print(result.n(45))
0.07385206378684
from scipy.stats import binom # Set the parameters for the binomial distribution n = 140 p = 0.25 # Calculate the probability of getting at least 43 correct answers probability = 1 - binom.cdf(42, n, p) print(probability)
0.07385206378684228

2.4. Classical & conditional probability

Probability theory is the study of random events and the likelihood of their occurrence.

It has wide-ranging applications in fields such as statistics, engineering, physics, and finance.

With the help of Sage, we can perform various calculations and simulations to gain a better understanding of probability concepts and solve complex problems.

Exercise (rolling a die)

Suppose we roll a fair six-sided die. What is the probability of rolling a 3?

Solution:

First, we define a probability space for the experiment of rolling a die, which consists of the set of possible outcomes {1, 2, 3, 4, 5, 6}, each of which has equal probability:

outcomes = FiniteEnumeratedSet([1,2,3,4,5,6])

To calculate the probability of rolling a 3, we can use the following command in SAGE:

1 / outcomes.cardinality()
1/6

This indicates that the probability of rolling a 3 is 1/6, which is consistent with our intuition.

Exercise (flipping a coin)

Suppose we flip a fair coin twice. What is the probability of getting exactly one head?

Solution:

We can define a probability space for this experiment as follows: {HH, HT, TH, TT}, where H represents heads and T represents tails.

toss_outcomes = FiniteEnumeratedSet(['H', 'T']) series_outcomes = cartesian_product([toss_outcomes, toss_outcomes]) print("All possible outcomes:", series_outcomes.list()) print("Each outcome has probability", 1 / series_outcomes.cardinality())
All possible outcomes: [('H', 'H'), ('H', 'T'), ('T', 'H'), ('T', 'T')] Each outcome has probability 1/4

There are two possible outcomes with only one head: HT and TH. To calculate the probability of getting exactly one head, we can use the following code in Sage:

2 / series_outcomes.cardinality()
1/2

Exercise

What is the probability that when rolling two dice the sum is 7, if we know that neither of the rolls resulted in a 2?

Solution:

Let BB be the event that neither of the rolls results into 2, and let AA be the event “sum is 7”.

The set of all possible outcomes is again denoted by Ω\Omega. Thus:

P(AB)=P(AB)P(B)=AB/ΩB/Ω=ABBP(A\mid B) = \frac{P(A\cap B)}{P(B)} = \frac{|A\cap B| / |\Omega|}{|B| / |\Omega|} = \frac{|A\cap B|}{|B|}.

The number 7 can appear as a sum in four ways if there is no 2 ((1,6), (3,4), (4,3), (6,1)), that is, AB=4,B=55=25|A \cap B| = 4, |B| = 5\cdot 5 = 25. Thus P(AB)=4/25P(A|B)= 4/25.

PAcondB = 4/25 print(n(PAcondB))
0.160000000000000

We can also run the simulation to estimate the resulting probability:

num_rolls = 100_000 #Generate num_rolls of two dice #random_element generates random int between x and y-1 inclusive die1=[IntegerRing().random_element(x=1, y=7) for _ in range(num_rolls)] die2=[IntegerRing().random_element(x=1, y=7) for _ in range(num_rolls)]
#Let's select only rolls where neither dice resulted in a 2 good_rolls = [] for x, y in zip(die1, die2): if x != 2 and y != 2: good_rolls.append((x, y))
#Let's calculate number of rolls with sum 7 count_sum_7 = 0 for x, y in good_rolls: if x + y == 7: count_sum_7 += 1 #The final probability is the ratio: result = count_sum_7 / len(good_rolls) print(n(result))
0.160287666307084

Exercise

Michaela has two mailboxes, one at gmail.com and the other at hotmail.com.

Her username is the same at both servers, but the passwords are different.

She does not remember which password corresponds to which server.

When typing in the password for accessing his mailbox, she makes a typo with probability 5%

(that is, if she tries to type in a specific password, she types what she intended with probability 95%).

At the server hotmail.com, Michaela typed in the username and a password, but the server told her that something is wrong.

What is the probability that she chose the correct password but just mistyped?

(Assume that the username is always typed correctly and that making a typo cannot turn wrong password into a good one.)

Solution:

Let AA be the event that Michaela typed in a wrong password at hotmail.com. This event is the union of two disjoint events:

  • A1A_1: she wanted to type in the correct password and mistyped,

  • A2A_2: he wanted to type in the wrong password (the one from gmail.com) and either mistyped it or not.

We are looking for a conditional probability P(A1A)P(A_1\mid A) which, according to the formula for conditional probability, is:

P(A1A)P(A)=P(A1)P(A1A2)=P(A1)P(A1)+P(A2)\frac{P(A_1\cap A)}{P(A)} = \frac{P(A_1)}{P(A_1\cup A_2)} = \frac{P(A_1)}{P(A_1)+P(A_2)}

where we have used P(A1A2)=P(A1)+P(A2)P(A_1 \cup A_2) = P(A_1) + P(A_2) since A1A_1 and A2A_2 are disjoint.

We just need to determine the probabilities P(A1)P(A_1) and P(A2)P(A_2).

The event A1A_1 is the intersection of two independent events: Michaela wanted to type in a correct password and Michaela mistyped.

According to the problem statement, the probability of the first event is 1/2 and the probability of the second event is 1/20.

In total P(A1)=1/21/20P(A_1)=1/2 \cdot 1/20 (we multiply the probabiliies, since the events are independent).

Further we have (directly from the problem statement) P(A2)=1/2P(A_2) = 1/2. In total P(A)=P(A1)+P(A2)P(A)=P(A_1)+P(A_2). We can now evaluate

P(A1A)=P(A1)P(A)P(A_1\mid A)= \frac{P(A_1)}{P(A)}
PA1 = 1/2 * 1/20 PA2 = 1/2 PA = PA1 + PA2 PA1condA = PA1 / PA print("P(A1|A)=",PA1condA)
P(A1|A)= 1/21

Geometric probability (optional)

The method of geometric probability can be used in the case that the given sample space consists of some region of a line, region, space (where we can determine (respectively) length, area, volume, etc).

We assume that the probability, is equal to the ratio of the area of the subregion to the area of the sample space.

Exercise

From Edinburgh Waverley station trains depart every hour (in the direction to Aberdeen).

From Aberdeen to Edinburgh they also depart every hour.

Assume that the trains move between these two stations with a uniform speed 72 km/h and are 100 meters long.

The trip takes 2 hrs in either direction.

The trains meet each other somewhere along the route. After visiting an Edinburgh pub, John, who lives in Aberdeen,

takes the train home and falls asleep at the departure.

During the trip from Edinburgh to Aberdeen he wakes up and randomly sticks his head out of the train for five seconds, on the side of the train where the trains travel in the opposite direction.

What is the probability that he loses his head? (We are assuming that there are no other trains involved.)

Solution:

The mutual speed of the oncoming trains is

2 * (72 * 1000 / (60 * 60))
40

fourty metres per second, the oncoming train passes John’s window for

n(100 / 40)
2.50000000000000

two and a half seconds. During John’s trip two trains pass by John’s window in the opposite direction.

Any overlap of the 2.5 seconds of the passing time interval with the 5 second time interval when John’s head might be sticking out is fatal.

Thus, for each train, the space of “favourable” outcomes is an interval of length 7.5 seconds somewhere in the sample space.

For two trains, it is double this amount.

The sample space of all outcomes is the interval 0,2×60×60⟨0, 2\times60\times60⟩.

Thus the probability of losing the head is:

result = 2 * 7.5 / (2 * 60 * 60) print(result)
0.00208333333333333

Exercise

In a certain country, a bus departs from town A to town B once a day at a random time between 8 a.m. and 8 p.m.

Once a day in the same time interval another bus departs in the opposite direction.

The trip in either direction takes 5 hours. What is the probability that the buses meet, assuming they use the same route?

Solution:

The sample space is a square 12×1212 \times 12 (1212 hours for each bus).

If we denote the time of the departure of the buses as xx and yy respectively, then they meet on the trail if and only if xy5|x−y| \le 5.

This inequality determines the region in the square of “favourable events”. Let's plot it:

num_points = 100_000 # Generate num_points points x=[RealField().random_element(min=8, max=20) for _ in range(num_points)] y=[RealField().random_element(min=8, max=20) for _ in range(num_points)]
# Split points into two groups corresponding to the case # where one piece of rod is less than 20 cm, and the rest favourable_events = [] other_events = [] for i, j in zip(x, y): if abs(i-j) <= 5: favourable_events.append((i, j)) else: other_events.append((i, j))
# Plots these two sets sing red and blue color scatter1 = scatter_plot(favourable_events, facecolor='blue', marker='o', markersize=10, edgecolor='blue') scatter2 = scatter_plot(other_events, facecolor='red', marker='o', markersize=10, edgecolor='red') # Show the plot show(scatter1+scatter2, axes=true, aspect_ratio=1, ticks=[1, 1])
Image in a Jupyter notebook

This is a complement to the union of two (red) right-angled isosceles triangles with legs of length 77. Their total area is 4949, so the area of the “favourable part” is 14449=95144 − 49 = 95. The probability is p=95144p = \frac{95}{144}.

print("probability = ", (12^2-7^2) / 12^2)
probability = 95/144

Exercise

A rod of length two meters is randomly divided into three parts. Determine the probability that at least one part is at most 20 cm long.

Solution:

Random division of a rod into three parts is given by two points of the cut, xx and yy (we first cut the rod in the distance xx from the origin, we do not move it and again cut it in the distance yy from the origin). The sample space is a square CC with side 2 m.

If we place the square CC so that its two sides lie on axes in the plane, then the condition that at least one part is at most 20 cm determines in the square a subregion OO: O={(x,y)Cx20x180y20y180xy20}.O = \{(x, y) \in C \mid x \le 20 \lor x \ge 180 \lor y ≤ 20 \lor y ≥ 180 \lor |x − y| \le 20\}.

Let's plot how this area looks like:

num_points = 100_000 #Generate num_points points x=[RealField().random_element(min=0, max=2) for _ in range(num_points)] y=[RealField().random_element(min=0, max=2) for _ in range(num_points)]
# Split points into two groups corresponding to the case # where one piece of rod is less than 20 cm, and the rest short_rod_cases = [] other_cases = [] for i, j in zip(x, y): if (i <= 0.2 or i > 1.8 or j <= 0.2 or j >= 1.8 or abs(i-j) <= 0.2): short_rod_cases.append((i, j)) else: other_cases.append((i, j))
# Plots these two sets using colors scatter1 = scatter_plot(short_rod_cases, facecolor='antiquewhite', marker='o', markersize=10, edgecolor='antiquewhite') scatter2 = scatter_plot(other_cases, facecolor='grey', marker='o', markersize=10, edgecolor='grey') # Show the plot show(scatter1+scatter2, axes=true, aspect_ratio=1, ticks=[0.2,0.2])
Image in a Jupyter notebook

We are interested in the ratio of the blue area to the total area of the square. Let's calculate its approximate value:

len(short_rod_cases) / (len(short_rod_cases) + len(other_cases))
0.50878

Using some simple geometry we can determine that the side of the red triangles is 1.41.4. Knowing this, we can calculate the area of two red triangles exactly as 1.422\frac{1.4^2}{2}. The answer to the problem will be the ratio of the blue area to the total area of the square: (2221.422)/22=221.4222=41.964=204400=51100(2^2-2\frac{1.4^2}{2})/2^2=\frac{2^2-1.4^2}{2^2}=\frac{4-1.96}{4}=\frac{204}{400}=\frac{51}{100}

print("red area = ", 1.4^2, ", final probability", 1-1.4^2/2^2)
red area = 1.96000000000000 , final probability 0.510000000000000

Exercise

Two meter-long pole is randomly divided into three pieces.

Determine the probability that a triangle can be built of the pieces.

Solution:

num_points = 100_000 #Generate num_points breaking boints x=[RealField().random_element(min=0, max=1) for _ in range(num_points)] y=[RealField().random_element(min=0, max=1) for _ in range(num_points)]
# Split the time points into two groups corresponding to the case # where three pieces form a triangle, and the otherwise good_pairs = [] bad_pairs = [] for i, j in zip(x, y): #lengths of the 3 pieces a = min(i, j) b = abs(i - j) c = 1 - max(i, j) # a + b + c - max([a,b,c]) gives the sum of two smaller pieces if a + b + c - max([a,b,c]) >= max([a,b,c]): good_pairs.append((i, j)) else: bad_pairs.append((i, j))
# Plots these two sets using red and blue color scatter1 = scatter_plot(good_pairs, facecolor='blue', marker='o', markersize=10, edgecolor='blue') scatter2 = scatter_plot(bad_pairs, facecolor='red', marker='o', markersize=10, edgecolor='red') # Show the plot show(scatter1+scatter2, axes=true, aspect_ratio=1, ticks=[0.1, 0.1])
Image in a Jupyter notebook
print("Estimated probability:", n(len(good_pairs) / len(x)))
Estimated probability: 0.250710000000000
result = 1/2^2 print("The exact probability is the area of two blue triangles:", result, ", which is about {:.3}".format(n(result)))
The exact probability is the area of two blue triangles: 1/4 , which is about 0.250

Exercise (Alice & Bob)

Alice and Bob play the following game. First, Alice tosses four coins. After that, Bob tosses five coins.

Bob wins if she gets strictly more heads than Alice.

What is the probability that Bob wins?

Solution:

With given constraint, there are only 29=5122^9=512 possible outcomes in the game.

toss_outcomes = [0, 1]*9 # 1 for head and 0 for tail toss_sequence = Arrangements(toss_outcomes, 9).list() len(toss_sequence)
512

We can generate all possible outcomes of 9 tosses and count cases where Alice wins.

count = 0 for seq in toss_sequence: count += sum(seq[:4]) < sum(seq[4:]) print('there are', count, 'sequences where Bob wins, so the probability is', count / len(toss_sequence))
there are 256 sequences where Bob wins, so the probability is 1/2

2.5. Plane Geometry

Let us now discuss basic tasks from Plane Geometry.

Points, lines and basic operations

# Define a point with coordinates (2,3) p = point((2,3), size=70) # Define another point with coordinates (4,1) q = point((4,1), size=70) # Define a line passing through points p and q l = line([(2, 3), (4, 1)])
# plot points and the line x_min = -5 x_max = 5 y_min = -5 y_max = 5 # Create a plot object line_plot = plot(l) PL = p+q+line_plot show(PL)
Image in a Jupyter notebook

Exercise

Find the length of the line segment joining the points (-1, 2) and (3, -4).

Solution:

We can use the distance formula to find the length of the line segment: d=((x2x1)2+(y2y1)2)d = \sqrt{((x2 - x1)^2 + (y2 - y1)^2)}

Substituting the coordinates of the two points, we get:

d = sqrt((3 - (-1))^2 + (-4 - 2)^2) show(d)

213\displaystyle 2 \, \sqrt{13}

Relative position of lines

Once we have the equations of the two lines, we can set them equal to each other to find the point of intersection.

If the two lines are not parallel, they will intersect at a single point.

Given two lines, the intersection point is the point at which the two lines cross each other. In slope-intercept form, we can write the equation of each line using the slope-intercept form, y = mx + b, where m is the slope and b is the y-intercept (the value of y when x is zero).

# function to find intersection, given slopes and intercept values of each def find_intersection(m1, b1, m2, b2): """Find the intersection point of two lines in slope-intercept form.""" # Check if the lines are parallel (i.e. have the same slope) if m1 == m2: return None # Lines are parallel and do not intersect # Compute the x-coordinate of the intersection point x = (b2 - b1) / (m1 - m2) # Compute the y-coordinate of the intersection point y = m1 * x + b1 return (x, y)

Exercise

Find the intersection point of the lines: y=2x+1,y=3x+4.y = 2x + 1,\quad y = -3x + 4.

Solution:

We can use the previous function, to find the intersection:

# Use the above function to find intersection intersection = find_intersection(2, 1, -3, 4) print(intersection)
(3/5, 11/5)

Here is a simpler way to understand the relative possition of two lines in the plane, just by solving the induced system:

Exercise (intersection point of 2 lines)

Find the intersection point of the lines given by the equations:

1(x)=2x+1,2(x)=x+3.\ell_1(x)=2x+1\,,\quad \ell_2(x)=-x+3\,.

Solution:

# Define the variables x, y = var("x, y") # Define the equations of the lines line1 = y == 2*x + 1 line2 = y == -x + 3 # Solve the system of equations solution = solve([line1, line2], x, y) solution
[[x == (2/3), y == (7/3)]]

Distance between two lines

The distance between the two parallel lines is the absolute value of the difference between their y-intercepts divided by the square root of 1 plus the square of their common slope.

Therefore the formula for the distance between two parallel lines in slope-intercept form is:

d=b2b1/(1+m12).d = |b2 - b1| / \sqrt{(1 + m1^2)}\,.
def find_distance_parallel(m1, b1, m2, b2): """Find the distance between two parallel lines in slope-intercept form.""" # Check if the lines are parallel (i.e. have the same slope) if m1 != m2: return None # Lines are not parallel # Compute the distance between the lines distance = abs(b2 - b1) / sqrt(1 + m1^2) return distance

Exercise (distance between lines)

Find the distance between the lines: y=2x+1,y=2x+5.y = 2x + 1, \quad y = 2x + 5.

Solution:

# Use the above function to find distance between two given lines distance = find_distance_parallel(2, 1, 2, 5) show(distance)

455\displaystyle \frac{4}{5} \, \sqrt{5}

Exercise

Perpendicular lines.

Two lines are perpendicular if and only if their slopes are negative reciprocals of each other. That is, m1 * m2 = -1

def are_perpendicular(m1, b1, m2, b2): """Check if two lines in slope-intercept form are perpendicular.""" # Check if the slopes are negative reciprocals of each other if m1 * m2 == -1: return True return False
# Check if y = 2x + 1 and y = -1/2 x + 5 are perpendicular perpendicular = are_perpendicular(2, 1, -1/2, 5) print(perpendicular)
True

Exercise (intersection of parametric lines)

Determine the intersection of the lines

p:x+y4=0,q:x=1+2t,y=2+t,tR.p : x + y − 4 = 0,\quad q : x = −1 + 2t, y = 2 + t, \quad t ∈\mathbb{R}.

Solution:

# Declare t as a symbol var('t') # Define the equations of the two lines p_eq = x + y - 4 == 0 q_eq1 = x == -1 + 2 * t q_eq2 = y == 2 + t # Substitute the second equation of line q into the first equation of line q and solve for t t_sol = y-2 # equation of line 2 in slope-intercept form q_subs = x == -1 + 2*(y-2) # use the function to find intersection find_intersection(-1, 4, 1/2, 5/2)
(1, 3)
# plot x_min = -5 x_max = 5 y_min = -5 y_max = 5 # Create a plot object p = plot([], xmin=x_min, xmax=x_max, ymin=y_min, ymax=y_max) # Add the first line to the plot p += plot(4-x, (x, x_min, x_max), color='blue', legend_label='Line 1') # Add the second line to the plot p += parametric_plot((t*2-1, t+2), (t, -5, 5), color='red', legend_label='Line 2') # Show the plot p.show()
Image in a Jupyter notebook

Exercise (slope of a line)

Find which of the below is the parametric equation P0+tuP_{0}+tu of the line having slope 53-\frac{5}{3} and passes through the point P0=[5,9]P_{0}=[5, -9]:

(i)[5,9]+t(3,5),(ii)[5,9]t(5,3),(iii)[3,5]+t(5,9),(iv)[5,9]+t(3,5).\begin{matrix} (i) & [5, -9]+t(-3, 5), & (ii) & [5, -9]-t(5, 3), \\ (iii) & [3, 5]+t(5, -9), & (iv) & [5, -9]+t(3, -5). \end{matrix}

Solution:

Given a line with slope mm its direction vector uu is given by u=(1,m)u=(1, m).

Hence, if the slope is a rational, i.e., m=pqm=\frac{p}{q} for some p,qZp, q\in\Z, then u=(q,p)u=(q, p).

It follows that the right answer is (i)(i) and (iv)(iv), which represent the same line.

It is interesting to present the figure of the lines in Sage, all together in one plot, with different colours and thickness/style, in order to be able to visualize them.

To do so, since the lines are given in parametric form we can use the command parametric_plot{\tt{parametric\_plot}}. We present the cell here:

t = var("t") P1= parametric_plot((5-3*t, -9+5*t), (t, -12, 12), thickness=4, color="green", legend_label="The right line") P2= parametric_plot((5-5*t, -9+3*t), (t, -10, 10), color="red", legend_label="A wrong line") P3= parametric_plot((3+5*t, 5-9*t), (t, -7, 9), color="blue", legend_label="A wrong line") P4= parametric_plot((5+3*t, -9-5*t), (t, -14, 14), linestyle="--", color="orange", legend_label="Also the right line") Q= point((5,-9), size=70) PL=P1+P2+P3+P4+Q show(PL)
Image in a Jupyter notebook

Exercise (Sage interactive environment)

Create a Sage interactive cell that allows the user to use two sliders m,bm, b and then plot the line y=ax+by=ax+b on R2\mathbb{R}^2.

Solution:

A Sage cell doing this operation has the form given below:

@interact def Lines(a = slider(-10, 10, 1), b = slider(-10, 10, 1, 0)): plot(a*x+b, -10,10).show()
Interactive function <function Lines at 0x7f05fcb16e60> with 2 widgets a: TransformIntSlider(value=0, descri…

Exercise

Find a parametric equation for a line in R3\mathbb{R}^3 given by the equations x2y+z=2,2x+yz=5.x − 2y + z = 2,\quad 2x + y − z = 5\,.

Solution:

To find a parametric equation for the line in R3, we need to first find a vector that lies on the line and a direction vector for the line. We can do this by solving the system of equations:

var("y") var("z") # Define the system of equations eq1= x-2*y+z-2 eq2= 2*x+y-z-5 eqns = [eq1, eq2] # Solve the system of equations sol = solve(eqns, x, y, z) sol
[[x == 1/5*r1 + 12/5, y == 3/5*r1 + 1/5, z == r1]]

This tells us that the line is given by the parametric equations:

x = 5/3t - 1/3s + 4/3 y = t z = s

where t and s are arbitrary parameters.

Exercise (area of triangle)

Find the area of the triangle with vertices at the points A=[0,0]A=[0, 0], B=[4,0]B=[4, 0] AND C=[4,3]C=[4, 3], using the Heron’s formula:

Area=s(sa)(sb)(sc),wheres=a+b+c2.\text{Area} = \sqrt{s(s - a)(s - b)(s - c)}\,,\quad \text{where}\quad s = \frac{a + b + c}{2}\,.

Solution:

On should first compute the side lengths:

a=distance between B(4,0) and C(4,3)=3,b=distance between A(0,0) and C(4,3)=5,c=distance between A(0,0) and B(4,0)=4.a = \text{distance between } B(4, 0) \text{ and } C(4, 3) = 3,\quad b = \text{distance between } A(0, 0) \text{ and } C(4, 3) = 5,\quad c = \text{distance between } A(0, 0) \text{ and } B(4, 0) = 4.

It follows that the area equals 6. Here is a verification

# Define the side lengths a = 3 b = 5 c = 4 # Semi-perimeter s = (a + b + c) / 2 # Heron's formula for the area area = sqrt(s * (s - a) * (s - b) * (s - c)) area
6

Exercise (barycenter)

If the vertices of A triangle are A(x1,y1)A(x_1, y_1), B(x2,y2)B(x_2, y_2) and C(x3,y3)C(x_3, y_3), the coordinates of the barycenter GG can be calculated using the following formula: G=(x1+x2+x33,y1+y2+y33). G=\left( \frac{x_1 + x_2 + x_3}{3}, \frac{y_1 + y_2 + y_3}{3} \right)\,. This formula essentially finds the average of the coordinates of the vertices.

Calculate the barycenter of a triangle with vertices A=[0,0]A=[0, 0], B=[4,0]B=[4, 0] and C=[2,3]C=[2, 3].

Solution:

# Define the vertices of the triangle A = (0, 0) # Vertex A B = (4, 0) # Vertex B C = (2, 3) # Vertex C # Calculate the coordinates of the barycenter (centroid) G = ((A[0] + B[0] + C[0]) / 3, (A[1] + B[1] + C[1]) / 3) # Output the barycenter G
(2, 1)

And here is the illustration:

# Define the vertices of the triangle A = (0, 0) # Vertex A B = (4, 0) # Vertex B C = (2, 3) # Vertex C # Barycenter computed previously G = (2, 1) # Barycenter # Create a plot of the triangle triangle_plot = polygon([A, B, C], color='lightblue', alpha=0.5) + point(G, color='red', size=30) # Add vertices triangle_plot += point(A, color='black', size=30) + point(B, color='black', size=30) + point(C, color='black', size=30) # Show the plot triangle_plot.show()
Image in a Jupyter notebook

Exercise (orthocenter)

The orthocenter of a triangle is one of its notable points, defined as the point where the three altitudes of the triangle intersect. An altitude is a line segment from a vertex perpendicular to the opposite side (or the line that extends that side).

Find the orthocenter of a triangle defined by the vertices A=[0,0]A=[0, 0], B=[4,0]B=[4, 0] and C=[2,3]C=[2, 3].

Moreover, plot the triangle and the orthocenter.

Solution:

First we compute the orthocenter, that we denote by HH: Here how it goes:

# Define the vertices of the triangle A = (0, 0) # Vertex A B = (4, 0) # Vertex B C = (2, 3) # Vertex C # Calculate the slopes of the sides m_AB = (B[1] - A[1]) / (B[0] - A[0]) # Slope of side AB m_BC = (C[1] - B[1]) / (C[0] - B[0]) # Slope of side BC m_AC = (C[1] - A[1]) / (C[0] - A[0]) # Slope of side AC # Calculate the equations of the altitudes # Altitude from vertex C to side AB (vertical line since AB is horizontal) altitude_C = (2, 0) # This is the x-coordinate of the altitude line # Altitude from vertex A to side BC # Slope of altitude from A (negative reciprocal of m_BC) slope_altitude_A = -1/m_BC # Equation of altitude from A: y - y_A = slope * (x - x_A) # Using A(0,0) to find the equation def altitude_A(x): return slope_altitude_A * (x - A[0]) + A[1] # Altitude from vertex B to side AC # Slope of altitude from B (negative reciprocal of m_AC) slope_altitude_B = -1/m_AC # Equation of altitude from B: y - y_B = slope * (x - x_B) # Using B(4,0) to find the equation def altitude_B(x): return slope_altitude_B * (x - B[0]) + B[1] # Solve for the intersection of altitude A and C # We can set y = altitude_A(x) equal to y = altitude_C and solve for x # This would normally involve solving equations, but we'll use the known points x_H = altitude_C[0] # x-coordinate of the orthocenter y_H = altitude_A(x_H) # Solve for y using altitude_A # The orthocenter H H = (x_H, y_H) # Output the orthocenter H
(2, 4/3)

To obtain the illustration is easy:

# Define the vertices of the triangle A = (0, 0) B = (4, 0) C = (2, 3) # The orthocenter H = (2, 4/3) # Create the triangle plot triangle_plot = polygon([A, B, C], color='lightblue', alpha=0.5, legend_label='Triangle ABC') point_A = point(A, color='red', size=30, legend_label='A(0, 0)') point_B = point(B, color='green', size=30, legend_label='B(4, 0)') point_C = point(C, color='blue', size=30, legend_label='C(2, 3)') point_H = point(H, color='purple', size=30, legend_label='Orthocenter H(2, 4/3)') # Add the points and triangle to the plot plot = triangle_plot + point_A + point_B + point_C + point_H # Show the plot with legends plot.show()
Image in a Jupyter notebook

Exercise

Determine whether or not the points [0, 2, 1], [−1, 2, 0], [−2, 5, 2] and [0, 5, 4] in R3\mathbb{R}^3 all lie in the same plane.

Solution:

To determine whether or not the given points lie in the same plane, we can check whether they are linearly dependent.

If the points are linearly dependent, then they all lie in the same plane.

# Define the points as vectors v1 = vector([0, 2, 1]) v2 = vector([-1, 2, 0]) v3 = vector([-2, 5, 2]) v4 = vector([0, 5, 4]) # Create a matrix with the vectors as rows M = matrix([v1, v2, v3, v4]) # Check if the rank of the matrix is less than 4 if M.rank() < 4: print("The points lie in the same plane.") else: print("The points do not lie in the same plane.")
The points lie in the same plane.

Exercise for practice:

Determine whether or not the points [2, 2, 1], [−1, 2, 2], [−2, 4, 1] and [0, 1, 3] in R3\mathbb{R}^3 all lie in the same plane.

Exercise

Find the distance between the point [1, 2, -1, 3] and the subspace U defined by the equations

{2x1+x2x3+3x4=0,  4x13x2+2x3x4=0,  x1x2+3x32x4=0}.\Big\{2x_1 + x_2 - x_3 + 3x_4 = 0, \ \ 4x_1 - 3x_2 + 2x_3 - x_4 = 0, \ \ x_1 - x_2 + 3x_3 - 2x_4 = 0\Big\}.

Solution:

To find the distance between a point and a subspace, we need to find the projection of the point onto the subspace and then calculate the distance between the original point and the projected point.

# define point and subspace point = vector([1, 2, -1, 3]) U = [vector([2, 1, -1, 3]), vector([4, -3, 2, -1]), vector([1, -1, 3, -2])] A = matrix(QQ, [[2, 1, -1, 3], [4, -3, 2, -1], [1, -1, 3, -2]]) b = vector([0, 0, 0]) # find orthogonal projection of the point onto the subspace: projection = A.right_kernel().basis()[0] # calculate distance between the point and the subspace distance = norm(point - projection) print(distance)
1/17*sqrt(5654)

2.6. Relations & Mappings

In this paragraph we will devote some space to describe how one can treat via SageMath relations and mappings.

Let us begin with relations. A relation can be represented as a set of ordered pairs:

R = {(1, 2), (2, 3), (3, 4), (4, 5)}

Exercise

Give the domain DD and the codomain II of the relation

R={(a,v),(b,x),(c,x),(c,u),(d,v),(f,y)}R = \{(a, v), (b, x), (c, x), (c, u), (d, v), (f, y)\}

between the sets A={a,b,c,d,e,f} A = \{a, b, c, d, e, f\} and B={x,y,u,v,w}B = \{x, y, u, v, w\}. Is the relation RR a mapping?

Solution:

To check if a relation is a mapping, we need to ensure that each element in the domain is connected to exactly one element in the codomain. Domain and codomain can be found like this:

A = set(['a', 'b', 'c', 'd', 'e','f']) B = set(['x','y', 'u', 'v', 'w']) Rel = set( [('a','v'), ('b','x'), ('c','x'),('c','u'), ('d','v'), ('f','y')]) # to compute the domain we iterate over R, taking the first element of each pair domain = set(x[0] for x in Rel) print(f'domain = {domain}') # to compute the codomain, we take the second element codomain = set(x[1] for x in Rel) print(f'codomain = {codomain}')
domain = {'d', 'a', 'f', 'b', 'c'} codomain = {'x', 'u', 'y', 'v'}
is_mapping = len(Rel) == len(domain) if is_mapping: print('Rel is a mapping') else: print('Rel is not a mapping')
Rel is not a mapping

Exercise

Determine for each of the following relations over the set {a,b,c,d}\{a, b, c, d\} whether it is an ordering and whether it is complete:

R1={(a,a),(b,b),(c,c),(d,d),(b,a),(b,c),(b,d)},R_1 = \{ (a,a), (b,b), (c,c), (d,d), (b,a), (b,c), (b,d) \},R2={(a,a),(b,b),(c,c),(d,d),(d,a),(a,d)},R_2 = \{ (a,a), (b,b), (c,c), (d,d), (d,a), (a,d) \},R3={(a,a),(b,b),(c,c),(d,d),(a,b),(b,c),(b,d)},R_3 = \{ (a,a), (b,b), (c,c), (d,d), (a,b), (b,c), (b,d) \},R4={(a,a),(b,b),(c,c),(a,b),(a,c),(a,d),(b,c),(b,d),(c,d)},R_4 = \{ (a,a), (b,b), (c,c), (a,b), (a,c), (a,d), (b,c), (b,d), (c,d) \},R5={(a,a),(b,b),(c,c),(d,d),(a,b),(a,c),(a,d),(b,c),(b,d),(c,d)}.R_5 = \{ (a,a), (b,b), (c,c), (d,d), (a,b), (a,c), (a,d), (b,c), (b,d), (c,d) \}.

Solution:

def is_reflexive(Rel, A=None): if not A: A = set(x[0] for x in Rel) return all({(a,a) in Rel for a in A}) def is_symmetric(Rel): return all({(b,a) in Rel for (a,b) in Rel}) def is_antisymmetric(Rel): # if a <= b and a != b then not (b <= a) return all({ (b,a) not in Rel for (a,b) in Rel if a != b}) def is_transitive(Rel): return all( { (a,d) in Rel for (a,b) in Rel for (c,d) in Rel if b == c}) def is_total(Rel, A=None): if not A: A = set(x[0] for x in Rel) return all( { (a,b) in Rel or (b,a) in Rel for a in A for b in A}) def is_linear_order(Rel): '''This function checks whether Rel is a linear order. Returns a bool value. ''' # A is the domain of Rel A = set(x[0] for x in Rel) # checking the properties of linear order reflexive = is_reflexive(Rel, A) antisymmetric = is_antisymmetric(Rel) transitive = is_transitive(Rel) total = is_total(Rel, A) linear_order = reflexive and antisymmetric and transitive and total return linear_order def relation_summary(Rel): A = set(x[0] for x in Rel) # here we simply use all the functions written above reflexive = is_reflexive(Rel, A) antisymmetric = is_antisymmetric(Rel) transitive = is_transitive(Rel) total = is_total(Rel, A) print(f'Reflexive: {reflexive}') print(f'Antisymmetric: {antisymmetric}') print(f'Transitive: {transitive}') print(f'Total: {total}')
Rel1 = {('a','a'),('b','b'),('c','c'), ('d','d'), ('b','a'), ('b','c'), ('b','d')} relation_summary(Rel1)
Reflexive: True Antisymmetric: True Transitive: True Total: False

Exercise

Give all the elements in SRS \circ R if R={(2,4),(4,4),(4,5)}, R = \{ (2,4), (4,4), (4,5) \}, S={(3,1),(3,2),(3,5),(4,1),(4,4)}. S = \{ (3,1), (3,2), (3,5), (4,1), (4,4) \}.

Solution:

def RelCompose(Rel1, Rel2): '''Returns the composition of two relations''' RS = set() for (a,b) in Rel1: for (c,d) in Rel2: if b == c: RS.add( (a,d)) return RS
R = {(2,4), (4,4), (4,5)} S = {(3,1), (3,2), (3,5), (4,1), (4,4)}
print(RelCompose(R,S))
{(4, 4), (2, 4), (4, 1), (2, 1)}

Exercise

Write down all the relations over a two-element set {1,2}\{1, 2\}, which are symmetric but are neither reflexive nor transitive.

Solution:

for rel in Subsets([(1,1),(2,2),(1,2),(2,1)]): if is_symmetric(rel) and (not is_transitive(rel)) and (not is_reflexive(rel)): print(set(rel))
{(1, 2), (2, 1)} {(1, 1), (1, 2), (2, 1)} {(1, 2), (2, 1), (2, 2)}

Exercise

Consider the set of numbers that have three digits in the ternary notation and a relation such that two numbers are in the relation whenever they

a) begin with the same two digits in this notation,

b) end with the same two digits in this notation.

Write down the corresponding equivalence classes.

Solution:

def equivalence_classes(s, eq): classes = {} for x in s: eq_class = None for c in classes: if eq(x, c[0]): eq_class = c break if eq_class is None: eq_class = (x,) classes[eq_class] = [] classes[eq_class].append(x) return list(classes.values())
from itertools import product letters = '012' strings = [''.join(p) for p in product(letters, repeat=3) if p[0] != '0'] strings
['100', '101', '102', '110', '111', '112', '120', '121', '122', '200', '201', '202', '210', '211', '212', '220', '221', '222']
# this is an alternative way to define relations in python: as a function from set of pairs to bool def have_common_23(s1, s2): ''' tests if two strings have the same 2nd and 3rd characters''' return s1[1:3] == s2[1:3] print(have_common_23('122','222')) print(have_common_23('212','210'))
True False
equivalence_classes(strings, have_common_23)
[['100', '200'], ['101', '201'], ['102', '202'], ['110', '210'], ['111', '211'], ['112', '212'], ['120', '220'], ['121', '221'], ['122', '222']]

Define function needed for the part (a) of this exercise. It should check if given strings have the same two first characters. Then call equivalence_classes with your function

2.7. Difference equations

In this final section we will treat difference equations, also known as recurrences equations.

A homogeneous linear difference equation (or homogeneous linear recurrence) of order k is given by the expression a0xn+a1xn1++akxnk=0,a00,akk, a_0 x_n + a_1 x_{n-1} + \dots + a_k x_{n-k} = 0, \quad a_0 \ne 0, a_k \ne k, the right hand side can be non-zero - in this case we have a nonhomogeneous linear difference equation: a0xn+a1xn1++akxnk=f(n),a00,akk. a_0 x_n + a_1 x_{n-1} + \dots + a_k x_{n-k} = f(n), \quad a_0 \ne 0, a_k \ne k.

# this is an example of solving nonhomogeneous difference equation x(n+2) - 5*x(n+1) + 6*x(n) = 7*n from sympy import Function, rsolve from sympy.abc import n x = Function('x') eq = x(n+2) - 5*x(n+1) + 6*x(n) - 7*n
rsolve(eq, x(n))

2nC0+3nC1+7C0n\displaystyle 2^{n} C_{0} + 3^{n} C_{1} + 7 C_{0} n

This is the general solution with two parameters C0,C1C_0, C_1. We can find a solution for some initial conditions like this:

rsolve(eq, x(n), {x(0): 1, x(1): 3})

3n\displaystyle 3^{n}

Exercise

Solve the linear difference equation: a(n)2a(n1)+a(n2)=0a(n) - 2a(n-1) + a(n-2) = 0, with initial conditions a(0)=1,a(1)=1a(0) = 1, a(1) = 1.

Solution:

a = Function('a') rsolve(a(n) - 2*a(n-1) + a(n-2), a(n), {0: 1, 1: 1} )

1\displaystyle 1

Exercise

Solve the nonlinear difference equation: a(n)=2a(n1)a(n2)+1a(n) = 2a(n-1) - a(n-2) + 1, with initial conditions a(0)=1,a(1)=3a(0) = 1, a(1) = 3.

Solution:

rsolve(a(n) - 2*a(n-1) + a(n-2) - 1, a(n), {0: 1, 1: 3} )

n22+3n2+1\displaystyle \frac{n^{2}}{2} + \frac{3 n}{2} + 1

Exercise

Find the formula for Fibonacci numbers. Check that F10=55F_{10} = 55 with Sage.

Solution:

fibs = rsolve(a(n+2) - a(n+1) - a(n), a(n), {0: 0, 1: 1}) fibs

5(1252)n5+5(12+52)n5\displaystyle - \frac{\sqrt{5} \left(\frac{1}{2} - \frac{\sqrt{5}}{2}\right)^{n}}{5} + \frac{\sqrt{5} \left(\frac{1}{2} + \frac{\sqrt{5}}{2}\right)^{n}}{5}

fibs.subs({n: 10}).simplify()

55\displaystyle 55

Remark

It can be useful to present a code in Sage which will return the first nn terms FnF_{n} of the Fibonacci sequence for a cerntain positive integer nn.

We can do this by defining first a function as below (see also the last section of Chapter 1).

def fib(n): if n <= 1: return 1 else: return fib(n-1) + fib(n-2)
[fib(i) for i in range(20)]
[1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765]

Or we can use the function fibonacci_sequence{\tt{fibonacci\_sequence}} that Sage provides, and allows us to list Fibonacci terms FnF_n, when nn runs to a certain region anba\leq n\leq b for a,bNa, b\in\mathbb{N}, as follows:

fibs = [i for i in fibonacci_sequence(1, 21)] fibs
[1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765]
fibs2 = [i for i in fibonacci_sequence(5, 10)] fibs2
[5, 8, 13, 21, 34]

Let us also mention how we can add Fibonacci terms:

sum([i for i in fibonacci_sequence(1, 5)]) # this will add the terms F_1, ..., F_4
7
sum([i for i in fibonacci_sequence(3, 5)])
5
fib(n)?

Exercise

Apply the command rsolve{\tt{rsolve}} to the nonlinear difference equation: a(n)=a(n1)(1a(n2))a(n) = a(n-1) * (1 - a(n-2)), with initial conditions a(0)=1,a(1)=1/2a(0) = 1, a(1) = 1/2.

Next, print the first 10 terms of resulting formula.

Solution:

sol = rsolve(a(n+2) - a(n+1) + a(n+1)*a(n), a(n), {a(0): 1, a(1): 1/2}) sol

in(12i4)+(i)n(12+i4)\displaystyle i^{n} \left(\frac{1}{2} - \frac{i}{4}\right) + \left(- i\right)^{n} \left(\frac{1}{2} + \frac{i}{4}\right)

# a(2) should be equal to 0 for i in range(10): print(f'a({i}) = {sol.subs(n, i).simplify()}') # rsolve does not work for nonlinear eqations
a(0) = 1 a(1) = 1/2 a(2) = -1 a(3) = -1/2 a(4) = 1 a(5) = 1/2 a(6) = -1 a(7) = -1/2 a(8) = 1 a(9) = 1/2