Introduction to SageMath

Institute of Mathematical Sciences, Chennai

by Ajit Kumar (ICT Mumbai)

Sage Days 126 (2024)

This is the SageMath worksheet used in the introductry lecture and sagemath in education.

What is SageMath???

  • SageMath is a free and open source Computer Algebra System (CAS).
  • It is build on top of many existing open-source packages.

    https://www.sagemath.org/links-components.html

  • SageMath distribution comprises of about 300 packages.

    https://github.com/sagemath/sage/tree/develop/build/pkgs

  • Mission:

    Creating a viable free open source alternative to

    Magma, Maple, Mathematica and Matlab

  • Native language of SageMath is PYTHON.
  • SageMath is a very nice mathematical tool for learning, teaching and researh in mathematics.
  • Sage has nice interface with $\LaTeX$
  • It is available for all platforms: Linux, MAC and Windows (through WSL)
  • Sage can be used online:

    CoCalc: CoCalc is a cloud-based collaborative software oriented towards research, teaching, and scientific publishing purposes.

    https://cocalc.com/

    SageCell: an easy-to-use web interface

    https://sagecell.sagemath.org

Why to use SageMath???

  • It is FREE with amazing capability.
  • Being Python based, it is easy to learn and use.
  • SageMath is a very nice pedagogical tool.
  • Useful for research in Mathematics.
  • SageMath can be used to explore most of the standard topics at UG and PG levels

Documentations/References??

https://doc.sagemath.org/

Getting Started with SageMath??

SageMath as an Advanced Calculator??

In??[1]:
5225+52152
Out[1]:
57377
In??[2]:
62358325*52537625
Out[2]:
3276158294478125
In??[3]:
527^52 # Raising power
Out[3]:
3420991468019301396709099078165893954774460027305174670899339207166277742791692606935376919696170202891188906522026349072469298089736788771521
In??[4]:
527**52 # Raising power
Out[4]:
3420991468019301396709099078165893954774460027305174670899339207166277742791692606935376919696170202891188906522026349072469298089736788771521
In??[5]:
527/52  # Returns the rational number
Out[5]:
527/52
In??[6]:
527/52.0
Out[6]:
10.1346153846154
In??[7]:
527//52 # Interger division or quotient
Out[7]:
10
In??[9]:
52
Out[9]:
52
In??[10]:
527%52 # Returns the remainder
Out[10]:
7
In??[11]:
sqrt(571)
Out[11]:
sqrt(571)
In??[12]:
sqrt(571.)
Out[12]:
23.8956062906970
In??[13]:
sqrt(571).n(digits=500)
Out[13]:
23.895606290697041027024598960419273218234209318776886205748348863252658749564111975288374701754911041957796822801793174362475921982754149890134288365991303256504735148720811921350913969247629375354202958995614760631834894382046297206842813343160719440830782731958797857238216948878657687573210702547242694745564339564388089074931202763974440326442884108750502375336744715224689568464946869015083114030741845870443872957300490205091149474834952141044306667324835262194979130225710792412199518967280387

Scientific functions

In??[1]:
sin(43.9),cos(-42.0),asin(0.25),exp(1.2),log(5242,10),ln(2.2)
Out[1]:
(-0.0822042844004434,
 -0.399985314988351,
 0.252680255142079,
 3.32011692273655,
 log(5242)/log(10),
 0.788457360364270)

Standard Mathematical Constants

In??[16]:
pi.n()
Out[16]:
3.14159265358979
In??[17]:
e.n()
Out[17]:
2.71828182845905
In??[18]:
golden_ratio.n()
Out[18]:
1.61803398874989
In??[19]:
euler_gamma.n()
Out[19]:
0.577215664901533
In??[20]:
R = RealField(200);R
Out[20]:
Real Field with 200 bits of precision
In??[21]:
R(e)
Out[21]:
2.7182818284590452353602874713526624977572470936999595749670
In??[22]:
R(pi)
Out[22]:
3.1415926535897932384626433832795028841971693993751058209749
In??[23]:
R(golden_ratio)
Out[23]:
1.6180339887498948482045868343656381177203091798057628621354
In??[24]:
R(euler_gamma)
Out[24]:
0.57721566490153286060651209008240243104215933593992359880577

The Euler Gamma constant is defined as $$\lim \left[\sum_{k=1}^n \frac{1}{k}-\ln(n)\right]$$

In??[25]:
n = 1000
sum([1/k for k in range(1,n+1)]).n()-ln(n).n()
Out[25]:
0.577715581568207

Working with integers??

In??[26]:
n = 9577545463761

Explore 'n DOT and press TAB' to get list of methods that can be applied to the object $n$.

In??[27]:
n.binary()
Out[27]:
'10001011010111110010001100011111011111010001'
In??[30]:
L = n.bits()
print(L)

[1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1]
In??[31]:
k=8
k.coprime_integers(12) # Returns all integers less than 12 which are co-prime to 8.
Out[31]:
[1, 3, 5, 7, 9, 11]
In??[32]:
# Get help on coprime_integers
n.coprime_integers?
In??[33]:
n.is_prime()
Out[33]:
False
In??[34]:
n.prime_factors()
Out[34]:
[3, 11, 293, 990541469]
In??[35]:
n.next_prime()
Out[35]:
9577545463771
In??[36]:
n.previous_prime()
Out[36]:
9577545463757
In??[38]:
n = 653
factorial(n)
Out[38]:
224188141946934793623319769929397902043660534721713367369751882391342327695431674742296200915307205380517864097248718203552328417282189419765104063744228512187915360010203211872842272650182673182591817765376536855581242564642524150945506794078495473734942718457984827654888089453387450047980425829045485876292169462112693698619692358263988025416751046872335675043618161143795141615409052873537778007691159675386859078172836487400968268905941417603536483920337588055481538009301366587882867377617719746638567515775126561199996711110089145458637119127851080328625754407068391377179030453450522509984468149383363752584154882200843508033301368641831701400249984992319871716869052811896532706801586381206733310424362801140194421388281401635262322580002252956165809446403534418930145452923154413978375522536944082510734184569489449113214286508031252408106860022797696042286928222954815760954405358823430747706876017761829826801580009942033487559717638961947366588469437043378673229381591294943409754253556056981801095560312451747984420882191818519606434273546974275682967234480431761427428119124761687059488355276688107945421274402450896520342997248787098319826985445334410373608840834323424359263101240987638275510836837756707292886469556877970496634478480203998920058678182834846490440550099644811215994109789440469326807623302363416077687943660619810510528743886550320062090284269901731748238065664000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
In??[39]:
N = factorial(n)
N.ndigits()
Out[39]:
1557
In??[40]:
N.digits().count(0) # Counts the number of zeroes in n!
Out[40]:
301
In??[41]:
n = factorial(1000000)  
len(n.digits())
Out[41]:
5565709
In??[42]:
z = Partitions(10^8).cardinality()  
z.ndigits()
Out[42]:
11132
In??[43]:
a,b = 60336, 70356
In??[44]:
gcd(a,b)
Out[44]:
12
In??[45]:
xgcd(a,b)
Out[45]:
(12, 1629, -1397)
In??[46]:
# Verification of Bezout's identity.
d,x,y = xgcd(a,b)
print(d,x,y)
d == x*a+y*b

12 1629 -1397
Out[46]:
True
In??[2]:
F.<a> = GF(49)
print(type(F))

<class 'sage.rings.finite_rings.finite_field_givaro.FiniteField_givaro_with_category'>
In??[3]:
R.<x> = PolynomialRing(F)
print(R)
f = x^6 + x^3 - 2*x+27
f.factor()

Univariate Polynomial Ring in x over Finite Field in a of size 7^2
Out[3]:
(x + 3) * (x + 4) * (x^2 + (a + 3)*x + 4*a + 2) * (x^2 + (6*a + 4)*x + 3*a + 6)
In??[4]:
F.<a> = GF(7^4)
f = [4*a^3, 2*a^3 + a^2 + 3*a + 5,
     3*a^3 + 5*a^2 + 4*a + 2, 2*a^3 + 2*a^2 + 2]
fd = F.dual_basis(f, check=True);
show(fd)
\(\displaystyle \left[3 a^{3} + 4 a^{2} + 6 a + 2, a^{3} + 6 a + 5, 3 a^{3} + 6 a^{2} + 2 a + 5, 5 a^{2} + 4 a + 3\right]\)
In??[5]:
F.gens()
Out[5]:
(a,)
In??[6]:
F.cardinality()
Out[6]:
2401
In??[49]:
F.list()[5]
Out[49]:
2*a^3 + 3*a^2 + 4*a
In??[9]:
F.dual_basis(fd,check=True)
Out[9]:
[4*a^3, 2*a^3 + a^2 + 3*a + 5, 3*a^3 + 5*a^2 + 4*a + 2, 2*a^3 + 2*a^2 + 2]

Standard rings and fields

In??[51]:
ZZ,QQ,RR,CC,AA,SR,QQbar,RDF,CDF, GF(7)
Out[51]:
(Integer Ring,
 Rational Field,
 Real Field with 53 bits of precision,
 Complex Field with 53 bits of precision,
 Algebraic Real Field,
 Symbolic Ring,
 Algebraic Field,
 Real Double Field,
 Complex Double Field,
 Finite Field of size 7)

Exploring `solve' function??

Help on solve

In??[86]:
# help(solve)
In??[53]:
x, y = var('x, y')
solve([x+y==6, x-y==4], x, y)
Out[53]:
[[x == 5, y == 1]]
In??[54]:
solve([x^2+y^2==6, 2*x-y==4], x, y)
Out[54]:
[[x == -1/5*sqrt(14) + 8/5, y == -2/5*sqrt(14) - 4/5], [x == 1/5*sqrt(14) + 8/5, y == 2/5*sqrt(14) - 4/5]]
In??[17]:
var('x')
solve([x^4+2>0,x^2-9<0],x)
Out[17]:
[[-3 < x, x < 3]]
In??[19]:
var('y,z')
solve([x^2 * y * z == 18, x * y^3 * z == 24, x * y * z^4 == 6], [x, y, z],solution_dict=True)
Out[19]:
[{x: 3, y: 2, z: 1},
 {x: 1.337215067329613 - 2.685489874065195*I,
  y: -1.7004342714592282 + 1.052864325754712*I,
  z: 0.9324722294043555 - 0.3612416661871523*I},
 {x: 1.3372150673296135 + 2.685489874065194*I,
  y: -1.7004342714592282 - 1.0528643257547123*I,
  z: 0.9324722294043555 + 0.3612416661871523*I},
 {x: -2.550651407188846 - 1.5792964886320715*I,
  y: -0.5473259801441661 + 1.9236512863456376*I,
  z: -0.9829730996839015 - 0.18374951781657012*I},
 {x: -2.550651407188845 + 1.5792964886320704*I,
  y: -0.5473259801441662 - 1.9236512863456383*I,
  z: -0.9829730996839015 + 0.18374951781657012*I},
 {x: 0.2768050783899189 - 2.9872025288850637*I,
  y: 1.478017834441328 - 1.3473912872931377*I,
  z: -0.8502171357296144 - 0.526432162877356*I},
 {x: 0.27680507838990626 + 2.9872025288851*I,
  y: 1.4780178344413182 + 1.3473912872931144*I,
  z: -0.8502171357296144 + 0.526432162877356*I},
 {x: -0.8209889702162458 + 2.885476929518458*I,
  y: -1.205269272758512 - 1.59603445456048*I,
  z: 0.09226835946330206 - 0.9957341762950345*I},
 {x: -0.8209889702162483 - 2.885476929518458*I,
  y: -1.2052692727585128 + 1.5960344545604788*I,
  z: 0.09226835946330206 + 0.9957341762950345*I},
 {x: -1.8079039091377576 - 2.3940516818407116*I,
  y: 0.8914767115530776 - 1.7903265827101338*I,
  z: 0.7390089172206591 - 0.6736956436465571*I},
 {x: -1.8079039091377695 + 2.3940516818407187*I,
  y: 0.8914767115530766 + 1.7903265827101247*I,
  z: 0.7390089172206591 + 0.6736956436465571*I},
 {x: 2.2170267516619795 + 2.0210869309396733*I,
  y: 1.8649444588087118 + 0.722483332374306*I,
  z: -0.2736629900720828 - 0.9618256431728189*I},
 {x: 2.2170267516620012 - 2.021086930939692*I,
  y: 1.8649444588086936 - 0.7224833323742995*I,
  z: -0.2736629900720828 + 0.9618256431728189*I},
 {x: 2.797416688213066 - 1.0837249985614603*I,
  y: -1.9659461993678036 + 0.36749903563314074*I,
  z: -0.6026346363792563 - 0.7980172272802396*I},
 {x: 2.7974166882130658 + 1.08372499856146*I,
  y: -1.965946199367804 - 0.36749903563314085*I,
  z: -0.6026346363792563 + 0.7980172272802396*I},
 {x: -2.9489192990517044 + 0.5512485534497115*I,
  y: 0.184536718926604 + 1.9914683525900692*I,
  z: 0.4457383557765383 - 0.8951632913550623*I},
 {x: -2.9489192990517044 - 0.5512485534497117*I,
  y: 0.18453671892660403 - 1.9914683525900694*I,
  z: 0.4457383557765383 + 0.8951632913550623*I}]

Working with functions of single variable??

In??[21]:
f(x) = exp(1-x^2)+sin(x^2)*(2-x)

Use `f DOT +TAB' to list methods that can be applied to f.

In??[22]:
f.plot(figsize=4)
Out[22]:
No description has been provided for this image
In??[24]:
f.plot(-10,10,color='red',figsize=4,gridlines='minor')
Out[24]:
No description has been provided for this image
In??[25]:
f.limit(x=1)
Out[25]:
x |--> sin(1) + 1
In??[26]:
f.limit(x=1,dir='+')
Out[26]:
x |--> sin(1) + 1
In??[27]:
f.limit(x=1,dir='-')
Out[27]:
x |--> sin(1) + 1
In??[28]:
f.diff()
Out[28]:
x |--> -2*(x - 2)*x*cos(x^2) - 2*x*e^(-x^2 + 1) - sin(x^2)
In??[29]:
f5=f.diff(5)(x)
print(f5)

-32*(x - 2)*x^5*cos(x^2) - 32*x^5*e^(-x^2 + 1) - 160*(x - 2)*x^3*sin(x^2) - 80*x^4*sin(x^2) + 160*x^3*e^(-x^2 + 1) + 120*(x - 2)*x*cos(x^2) + 240*x^2*cos(x^2) - 120*x*e^(-x^2 + 1) + 60*sin(x^2)
In??[30]:
latex(f5)
Out[30]:
-32 \, {\left(x - 2\right)} x^{5} \cos\left(x^{2}\right) - 32 \, x^{5} e^{\left(-x^{2} + 1\right)} - 160 \, {\left(x - 2\right)} x^{3} \sin\left(x^{2}\right) - 80 \, x^{4} \sin\left(x^{2}\right) + 160 \, x^{3} e^{\left(-x^{2} + 1\right)} + 120 \, {\left(x - 2\right)} x \cos\left(x^{2}\right) + 240 \, x^{2} \cos\left(x^{2}\right) - 120 \, x e^{\left(-x^{2} + 1\right)} + 60 \, \sin\left(x^{2}\right)

The fifth derivative of $f$ is $$-32 \, {\left(x - 2\right)} x^{5} \cos\left(x^{2}\right) - 32 \, x^{5} e^{\left(-x^{2} + 1\right)} - 160 \, {\left(x - 2\right)} x^{3} \sin\left(x^{2}\right) - 80 \, x^{4} \sin\left(x^{2}\right) + 160 \, x^{3} e^{\left(-x^{2} + 1\right)} + 120 \, {\left(x - 2\right)} x \cos\left(x^{2}\right) + 240 \, x^{2} \cos\left(x^{2}\right) - 120 \, x e^{\left(-x^{2} + 1\right)} + 60 \, \sin\left(x^{2}\right)$$

In??[64]:
f.taylor(x,0,4)
Out[64]:
x |--> 1/2*x^4*e - x^3 - x^2*(e - 2) + e
In??[65]:
f.integral(x,0,4).n()
Out[65]:
2.92445245934350
In??[66]:
numerical_integral(f,0,4)
Out[66]:
(2.924452459343499, 4.900169303291803e-14)
In??[67]:
f.find_root(2,5)
Out[67]:
3.963327271283327
In??[68]:
f.find_local_maximum(2,3)
Out[68]:
(0.26505732334709253, 2.302028937091704)
In??[69]:
f.find_local_minimum(-1,0)
Out[69]:
(2.6908391241742633, -0.31522616679466264)
In??[70]:
## Symbolic integral
var("c b k")
assume(b>0)
assume(c>0)
integral(1/( x^2+2*k*x+(k^2+b)), x , -c, c).show()
\(\displaystyle \frac{\arctan\left(\frac{c + k}{\sqrt{b}}\right)}{\sqrt{b}} - \frac{\arctan\left(-\frac{c - k}{\sqrt{b}}\right)}{\sqrt{b}}\)
In??[33]:
# Improper integral
h(x)=exp(-x^2)
h.integral(x,0,infinity).show()
\(\displaystyle \frac{1}{2} \, \sqrt{\pi}\)

Implicit Derivative??

In??[72]:
g(x,y)=x^3+y^3+3*x*y+x-y
x0=-1.5
h(y)=g.substitute(x=x0)
y0=h.find_root(-2,2)
dyx=g.implicit_derivative(y,x)
m=dyx.subs(x=x0,y=y0)
l(x)=y0+m*(x-x0)
print('The derivative is give by')
show(dyx)
implicit_plot(g(x,y),(x,-3,2),(y,-3,2),color='red',figsize=5,gridlines='minor')+\
point((x0,y0),pointsize=30)+plot(l(x),(x,x0-0.2,x0+0.2))

The derivative is give by
\(\displaystyle -\frac{3 \, x^{2} + 3 \, y + 1}{3 \, y^{2} + 3 \, x - 1}\)
Out[72]:
No description has been provided for this image

Multivariable functions??

In??[34]:
var('x,y')
f(x,y) = (2*x^2-y^2)*exp(-x^2-2*y^2)
In??[??]:
plot3d(f(x,y),(x,-2,2),(y,-2,2))
In??[36]:
contour_plot(f(x,y),(x,-2,2),(y,-2,2),contours=20)
Out[36]:
No description has been provided for this image
In??[76]:
contour_plot(f(x,y),(x,-2,2),(y,-1.5,1.5),fill=False,contours=30,labels=True)
Out[76]:
No description has been provided for this image
In??[??]:
## Animating contours
animate( [contour_plot(f(x,y) , (x,-2,2),
(y,-1.5,1.5), contours=[k], fill=false ) for k in srange(-0.8,0.8,0.05)] )
In??[37]:
grad = f.gradient()
show(grad(x,y))
\(\displaystyle \left(-2 \, {\left(2 \, x^{2} - y^{2}\right)} x e^{\left(-x^{2} - 2 \, y^{2}\right)} + 4 \, x e^{\left(-x^{2} - 2 \, y^{2}\right)},\,-4 \, {\left(2 \, x^{2} - y^{2}\right)} y e^{\left(-x^{2} - 2 \, y^{2}\right)} - 2 \, y e^{\left(-x^{2} - 2 \, y^{2}\right)}\right)\)
In??[38]:
Hf = f.hessian()
show(Hf(x,y))
\(\displaystyle \left(\begin{array}{rr} 4 \, {\left(2 \, x^{2} - y^{2}\right)} x^{2} e^{\left(-x^{2} - 2 \, y^{2}\right)} - 16 \, x^{2} e^{\left(-x^{2} - 2 \, y^{2}\right)} - 2 \, {\left(2 \, x^{2} - y^{2}\right)} e^{\left(-x^{2} - 2 \, y^{2}\right)} + 4 \, e^{\left(-x^{2} - 2 \, y^{2}\right)} & 8 \, {\left(2 \, x^{2} - y^{2}\right)} x y e^{\left(-x^{2} - 2 \, y^{2}\right)} - 12 \, x y e^{\left(-x^{2} - 2 \, y^{2}\right)} \\ 8 \, {\left(2 \, x^{2} - y^{2}\right)} x y e^{\left(-x^{2} - 2 \, y^{2}\right)} - 12 \, x y e^{\left(-x^{2} - 2 \, y^{2}\right)} & 16 \, {\left(2 \, x^{2} - y^{2}\right)} y^{2} e^{\left(-x^{2} - 2 \, y^{2}\right)} + 16 \, y^{2} e^{\left(-x^{2} - 2 \, y^{2}\right)} - 4 \, {\left(2 \, x^{2} - y^{2}\right)} e^{\left(-x^{2} - 2 \, y^{2}\right)} - 2 \, e^{\left(-x^{2} - 2 \, y^{2}\right)} \end{array}\right)\)
In??[39]:
jacobian(f,[x,y])(x,y)
Out[39]:
[-2*(2*x^2 - y^2)*x*e^(-x^2 - 2*y^2) + 4*x*e^(-x^2 - 2*y^2) -4*(2*x^2 - y^2)*y*e^(-x^2 - 2*y^2) - 2*y*e^(-x^2 - 2*y^2)]
In??[40]:
contr = contour_plot(f(x,y),(x,-2,2),(y,-1.5,1.5),fill=False,contours=20)
grad_plot = plot_vector_field(grad,(x,-2,2),(y,-1.5,1.5),color='red')
show(contr+grad_plot,figsize=6)
No description has been provided for this image

Barth Surface $${\scriptsize f(x,y,z)={\left(x^{2} + y^{2} + z^{2} + t - 2\right)}^{2} {\left(x^{2} + y^{2} + z^{2} - 1\right)}^{2} {\left(5 \, t + 3\right)} - 8 \, {\left(t^{4} x^{2} - z^{2}\right)} {\left(t^{4} y^{2} - x^{2}\right)} {\left(t^{4} z^{2} - y^{2}\right)} {\left(x^{4} - 2 \, x^{2} y^{2} + y^{4} - 2 \, x^{2} z^{2} - 2 \, y^{2} z^{2} + z^{4}\right)}}$$ where $t$ is a parameter, typically, the golden ratio.

In??[??]:
## Barth surface
reset()
var('x,y,z,t')
t=(sqrt(5)-1)/2
f(x, y, z) = (3 + 5*t)*(-1 + x^2 + y^2 + z^2)^2*(-2 + t + x^2 + y^2 + z^2)^2+\
    8*(x^2-t^4*y^2)*(-(t^4*x^2) + z^2)*(y^2-t^4*z^2)*(x^4-2*x^2*y^2 + y^4-2*x^2*z^2-2*y^2*z^2 + z^4)
implicit_plot3d(f(x,y,z),(x,-2,2),(y,-2,2),(z,-2,2),plot_points=150)
In??[??]:
t=(sqrt(5)+1)/2
f(x, y, z) = (3 + 5*t)*(-1 + x^2 + y^2 + z^2)^2*(-2 + t + x^2 + y^2 + z^2)^2+\
    8*(x^2-t^4*y^2)*(-(t^4*x^2) + z^2)*(y^2-t^4*z^2)*(x^4-2*x^2*y^2 + y^4-2*x^2*z^2-2*y^2*z^2 + z^4)
implicit_plot3d(f(x,y,z),(x,-2,2),(y,-2,2),(z,-2,2),plot_points=150)
In??[??]:
t=1
f(x, y, z) = (3 + 5*t)*(-1 + x^2 + y^2 + z^2)^2*(-2 + t + x^2 + y^2 + z^2)^2+\
    8*(x^2-t^4*y^2)*(-(t^4*x^2) + z^2)*(y^2-t^4*z^2)*(x^4-2*x^2*y^2 + y^4-2*x^2*z^2-2*y^2*z^2 + z^4)
implicit_plot3d(f(x,y,z),(x,-2,2),(y,-2,2),(z,-2,2),plot_points=150)

Piriform Surface $$f(x,y,z)=x^4-x^3+y^2+z^2$$

In??[??]:
var('x,y,z')
f(x,y,z)=x^4-x^3+y^2+z^2
piriform = implicit_plot3d(f(x,y,z),(x,-0.5,1),(y,-0.5,1),(z,-0.5,0.5),plot_points=100,frame=True)
show(piriform)
In??[??]:
# solve(f(x,y,z)==0,[x,y,z])
In??[86]:
a1=0.7;b1=0.2;c1=-sqrt(-a1^4 + a1^3 - b1^2);c1
#a2=0.7;b2=-0.3;c2=-sqrt(-a2^4 + a2^3 - b2^2);c2
Out[86]:
-0.250798724079689
In??[??]:
##Tangennt Plane
T(x,y,z)=(x-a1)*f.diff(x)(x=a1,y=b1,z=c1)+(y-b1)*f.diff(y)(x=a1,y=b1,z=c1)+(z-c1)*f.diff(z)(x=a1,y=b1,z=c1);
pt = point3d((a1,b1,c1),color='red',size=5)
ep=0.2
Tgt=implicit_plot3d(T(x,y,z),(x,a1-ep,a1+ep),(y,b1-ep,b1+ep),(z,c1-ep,c1+ep),color='green',opacity=0.5)
show(piriform+pt+Tgt,frame=False)

Plotting in SageMath??

2d plotting??

https://doc.sagemath.org/html/en/reference/plotting/index.html

  • plot()
  • parametric_plot()
  • implicit_plot()
  • polar_plot()
  • region_plot()
  • list_plot()
  • scatter_plot()
  • bar_chart()
  • contour_plot()
  • density_plot()
  • plot_vector_field()
  • plot_slope_field()
  • matrix_plot()
  • complex_plot()
  • graphics_array()
  • multi_graphics()
  • plot_step_function()
In??[42]:
plot_step_function([(i, prime_pi(i)) for i in range(20)],figsize=4)
Out[42]:
No description has been provided for this image

Polar Plot??

In??[43]:
t = var('t'); n = 20/19
r1 = polar_plot(1 + 1/3*cos(n*t), (t, 0, n*36*pi),plot_points = 5000,figsize=4)
r2 = polar_plot(1 + 2*cos(n*t), (t, 0, n*36*pi), plot_points= 5000,figsize=4)
graphics_array([r1,r2])
Out[43]:
No description has been provided for this image
In??[90]:
C = graphs.CubeGraph(8)
P = C.graphplot(vertex_labels=False, vertex_size=0,
                graph_border=True)
P.show()
No description has been provided for this image
In??[91]:
A = random_matrix(ZZ, 100000, density=.00001, sparse=True)
matrix_plot(A, marker=',')
Out[91]:
No description has been provided for this image

3d Plotting??

https://doc.sagemath.org/html/en/reference/plot3d/index.html

Linear Algera with SageMath??

In??[44]:
v = vector([2,-3,1])
w = vector([3,1,-2])
In??[45]:
u = w+v
In??[??]:
v.plot()+w.plot(color='red')+u.plot(color='green')
In??[46]:
v.dot_product(u)
Out[46]:
15
In??[47]:
v.cross_product(w)
Out[47]:
(5, 7, 11)
In??[48]:
# Various norms
u.norm(),u.norm(1),v.norm(infinity)
Out[48]:
(sqrt(30), 8, 3)
In??[49]:
V = VectorSpace(QQ,6);V
Out[49]:
Vector space of dimension 6 over Rational Field
In??[50]:
v1 = vector(QQ,[-2,1,3,4,2,7])
v2 = vector(QQ,[1,2,3,-1,2,-1])
v3 = vector(QQ,[-5, 0, 3, 9, 2, 15])
v4 = vector(QQ,[2,1,0,5,-3,1])
v5 = vector(QQ,[11, 3, -3, 6, -11, -12])
v6 = vector(QQ,[2,3,5,7,11,13])
In??[51]:
B = [v1,v2,v3,v4,v5,v6]
In??[52]:
V.linear_dependence(B)
Out[52]:
[
(1, -1/2, 0, -3/2, 1/2, 0),
(0, 0, 1, -3, 1, 0)
]

We get list of scalares $\alpha_i$ so that $\sum \alpha_iv_i =0$.

In??[61]:
## Verification 
??= V.linear_dependence(B)
sum([??[0][i]*B[i] for i in range(len(B))])
Out[61]:
(0, 0, 0, 0, 0, 0)
In??[62]:
V.linear_dependence([v1,v2,v4])
Out[62]:
[

]
In??[63]:
V.span(B)
Out[63]:
Vector space of degree 6 and dimension 4 over Rational Field
Basis matrix:
[      1       0       0       0  293/61  171/61]
[      0       1       0       0 -824/61 -641/61]
[      0       0       1       0  496/61  374/61]
[      0       0       0       1   11/61   72/61]
In??[64]:
column_matrix(B).T.rref()
Out[64]:
[      1       0       0       0  293/61  171/61]
[      0       1       0       0 -824/61 -641/61]
[      0       0       1       0  496/61  374/61]
[      0       0       0       1   11/61   72/61]
[      0       0       0       0       0       0]
[      0       0       0       0       0       0]
In??[65]:
A = column_matrix(B).T;A
Out[65]:
[ -2   1   3   4   2   7]
[  1   2   3  -1   2  -1]
[ -5   0   3   9   2  15]
[  2   1   0   5  -3   1]
[ 11   3  -3   6 -11 -12]
[  2   3   5   7  11  13]
In??[66]:
A.row_space()
Out[66]:
Vector space of degree 6 and dimension 4 over Rational Field
Basis matrix:
[      1       0       0       0  293/61  171/61]
[      0       1       0       0 -824/61 -641/61]
[      0       0       1       0  496/61  374/61]
[      0       0       0       1   11/61   72/61]
In??[67]:
A.column_space()
Out[67]:
Vector space of degree 6 and dimension 4 over Rational Field
Basis matrix:
[ 1  0  2  0 -2  0]
[ 0  1 -1  0  1  0]
[ 0  0  0  1  3  0]
[ 0  0  0  0  0  1]
In??[68]:
A.left_kernel()
Out[68]:
Vector space of degree 6 and dimension 2 over Rational Field
Basis matrix:
[   1 -1/2    0 -3/2  1/2    0]
[   0    0    1   -3    1    0]
In??[69]:
A.nullity()
Out[69]:
2
In??[70]:
M = column_matrix([V.random_element() for i in range(6)]).T
M
Out[70]:
[    0     0   3/2   1/3   1/4     0]
[ -1/8     0  -1/7     0     1     2]
[   -1     6     9  -1/4 -1/10    -2]
[  586    -2  1/11   1/5   1/2     0]
[    1     1     1     7     4     0]
[   -1     1 -1/12     0     0    -1]
In??[75]:
G1,M1 = M.gram_schmidt()
In??[76]:
M1*G1
Out[76]:
[    0     0   3/2   1/3   1/4     0]
[ -1/8     0  -1/7     0     1     2]
[   -1     6     9  -1/4 -1/10    -2]
[  586    -2  1/11   1/5   1/2     0]
[    1     1     1     7     4     0]
[   -1     1 -1/12     0     0    -1]
In??[73]:
# Use help on M.gram_schmidt to know more about it.
In??[77]:
m = random_matrix(RDF, 1000) # 1000 X 1000 random matrix over RDF
In??[78]:
e = m.eigenvalues()  
w = [(i, abs(e[i])) for i in range(len(e))]
show(points(w,size=3,color='red'))
No description has been provided for this image

Jordan Canonical Form??

Example: Find a Jordan Canonical form of $A=\left(\begin{array}{rrrrrrrrrrr} 3 & 0 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 3 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 3 & 2 & 1 & 1 & 0 & -3 & 1 \\ 0 & 0 & 0 & 0 & -1 & 2 & 0 & 0 & 1 & -1 & 1 \\ 0 & 0 & 0 & 0 & 1 & 1 & 3 & 0 & -6 & 5 & -1 \\ 0 & 0 & 0 & 0 & 0 & 2 & 1 & 4 & 1 & -3 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 3 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -2 & 5 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -2 & 0 & 5 \end{array}\right)$

In??[79]:
A = matrix([[3 , 0 , 1 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0],
            [1 , 3 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0],
            [0 , 0 , 3 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0],
            [0 , 0 , 0 , 3 , 0 , 0 , 0 , 0 , 0 , 0 , 0],
            [0 , 0 , 0 , 0 , 3 , 2 , 1 , 1 , 0 , -3 , 1],
            [0 , 0 , 0 , 0 , -1 , 2 , 0 , 0 , 1 , -1 , 1],
            [0 , 0 , 0 , 0 , 1 , 1 , 3 , 0 , -6 , 5 , -1],
            [0 , 0 , 0 , 0 , 0 , 2 , 1 , 4 , 1 , -3 , 1],
            [0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 3 , 0 , 0],
            [0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , -2 , 5 , 0],
            [0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , -2 , 0 , 5 ]])
In??[80]:
p(x) = A.charpoly();p(x)
Out[80]:
x^11 - 37*x^10 + 619*x^9 - 6183*x^8 + 40986*x^7 - 189378*x^6 + 622566*x^5 - 1456542*x^4 + 2377269*x^3 - 2578473*x^2 + 1673055*x - 492075
In??[81]:
p.factor()
Out[81]:
x |--> (x - 3)^9*(x - 5)^2
In??[82]:
A.eigenvalues()
Out[82]:
[5, 5, 3, 3, 3, 3, 3, 3, 3, 3, 3]
In??[83]:
A.minpoly().factor()
Out[83]:
(x - 5) * (x - 3)^4
In??[84]:
A.jordan_form()
Out[84]:
[5|0|0 0 0 0|0 0 0|0 0]
[-+-+-------+-----+---]
[0|5|0 0 0 0|0 0 0|0 0]
[-+-+-------+-----+---]
[0|0|3 1 0 0|0 0 0|0 0]
[0|0|0 3 1 0|0 0 0|0 0]
[0|0|0 0 3 1|0 0 0|0 0]
[0|0|0 0 0 3|0 0 0|0 0]
[-+-+-------+-----+---]
[0|0|0 0 0 0|3 1 0|0 0]
[0|0|0 0 0 0|0 3 1|0 0]
[0|0|0 0 0 0|0 0 3|0 0]
[-+-+-------+-----+---]
[0|0|0 0 0 0|0 0 0|3 1]
[0|0|0 0 0 0|0 0 0|0 3]
In??[85]:
J,P = A.jordan_form(transformation=True)
In??[121]:
P.inverse()*A*P
Out[121]:
[5 0 0 0 0 0 0 0 0 0 0]
[0 5 0 0 0 0 0 0 0 0 0]
[0 0 3 1 0 0 0 0 0 0 0]
[0 0 0 3 1 0 0 0 0 0 0]
[0 0 0 0 3 1 0 0 0 0 0]
[0 0 0 0 0 3 0 0 0 0 0]
[0 0 0 0 0 0 3 1 0 0 0]
[0 0 0 0 0 0 0 3 1 0 0]
[0 0 0 0 0 0 0 0 3 0 0]
[0 0 0 0 0 0 0 0 0 3 1]
[0 0 0 0 0 0 0 0 0 0 3]

Google PageRank??

In??[122]:
n = int(input('Enter the number of website:'))
N = int(input('Enter the number of links:'))

Enter the number of website:10
Enter the number of links:75
In??[123]:
D = digraphs.RandomDirectedGNM(n,N)
D.show(figsize=6)
No description has been provided for this image
In??[124]:
H=D.adjacency_matrix()
H=H.T
print( "Ajacency matrix")
print(H)

Ajacency matrix
[0 1 1 1 1 1 1 1 1 1]
[0 0 1 1 1 1 1 1 1 1]
[1 1 0 1 1 1 1 1 1 1]
[1 1 1 0 1 1 1 0 1 1]
[0 0 1 1 0 0 1 1 1 1]
[0 0 0 1 1 0 1 1 1 1]
[1 1 1 1 1 1 0 1 1 1]
[1 1 1 1 1 0 1 0 1 1]
[0 1 1 1 0 1 1 1 0 1]
[1 1 1 1 0 0 0 1 0 0]
In??[125]:
## Probability matrix P=[p_{ij}] where p_{ij} is the probability of visiting i th website from jth website.
print ("Probability matrix or transition matrix")
P=column_matrix([H.columns()[i]/max(sum(H.columns()[i]),1) for i in range(n)])
print(P)

Probability matrix or transition matrix
[  0 1/7 1/8 1/9 1/7 1/6 1/8 1/8 1/8 1/9]
[  0   0 1/8 1/9 1/7 1/6 1/8 1/8 1/8 1/9]
[1/5 1/7   0 1/9 1/7 1/6 1/8 1/8 1/8 1/9]
[1/5 1/7 1/8   0 1/7 1/6 1/8   0 1/8 1/9]
[  0   0 1/8 1/9   0   0 1/8 1/8 1/8 1/9]
[  0   0   0 1/9 1/7   0 1/8 1/8 1/8 1/9]
[1/5 1/7 1/8 1/9 1/7 1/6   0 1/8 1/8 1/9]
[1/5 1/7 1/8 1/9 1/7   0 1/8   0 1/8 1/9]
[  0 1/7 1/8 1/9   0 1/6 1/8 1/8   0 1/9]
[1/5 1/7 1/8 1/9   0   0   0 1/8   0   0]
In??[126]:
i = int(input('Enter the starting website index:'))
L = [0]*n
L[i-1]=1
v = vector(L);v

Enter the starting website index:5
Out[126]:
(0, 0, 0, 0, 1, 0, 0, 0, 0, 0)
In??[127]:
S=1/n*matrix([[1]*n]*n)
#print(S)
print("The Google Matrix of the associated network is given by\n")
alpha=0.85 ## damping factor
G = alpha*P+(1-alpha)*S; 
print(G.n(digits=4))

The Google Matrix of the associated network is given by

[0.01500  0.1364  0.1213  0.1094  0.1364  0.1567  0.1213  0.1213  0.1213  0.1094]
[0.01500 0.01500  0.1213  0.1094  0.1364  0.1567  0.1213  0.1213  0.1213  0.1094]
[ 0.1850  0.1364 0.01500  0.1094  0.1364  0.1567  0.1213  0.1213  0.1213  0.1094]
[ 0.1850  0.1364  0.1213 0.01500  0.1364  0.1567  0.1213 0.01500  0.1213  0.1094]
[0.01500 0.01500  0.1213  0.1094 0.01500 0.01500  0.1213  0.1213  0.1213  0.1094]
[0.01500 0.01500 0.01500  0.1094  0.1364 0.01500  0.1213  0.1213  0.1213  0.1094]
[ 0.1850  0.1364  0.1213  0.1094  0.1364  0.1567 0.01500  0.1213  0.1213  0.1094]
[ 0.1850  0.1364  0.1213  0.1094  0.1364 0.01500  0.1213 0.01500  0.1213  0.1094]
[0.01500  0.1364  0.1213  0.1094 0.01500  0.1567  0.1213  0.1213 0.01500  0.1094]
[ 0.1850  0.1364  0.1213  0.1094 0.01500 0.01500 0.01500  0.1213 0.01500 0.01500]
In??[128]:
k = int(input('Enter the number of clicks:'))
pn = G^k*v.n(digits=6)
print('After ' +str(k)+' iteration the probability vector is: \n');
pn

Enter the number of clicks:200
After 200 iteration the probability vector is:
Out[128]:
(0.112500, 0.100319, 0.118983, 0.109670, 0.0797546, 0.0767971, 0.118983, 0.109149, 0.0929407, 0.0809033)
In??[129]:
## Sorting the website as per rank
PN = flatten(pn)
sorted_nodes = []
for i in range(n):
    sorted_nodes.append([PN[i],i])
sorted_nodes.sort(reverse=True)    
print(sorted_nodes)

[[0.118983, 6], [0.118983, 2], [0.112500, 0], [0.109670, 3], [0.109149, 7], [0.100319, 1], [0.0929407, 8], [0.0809033, 9], [0.0797546, 4], [0.0767971, 5]]

Optimization in SageMath??

Example: Find a minimizer of $$ -12 \, x^{4} - 24 \, x^{2} y^{2} - 12 \, y^{4} + 25 \, x^{2} - 6 \, x y + 25 \, y^{2}$$ starting with $(1,0.5)$ using the random/walk search method.

In??[130]:
var = ('x,y,z')
f(x,y)=25*x^2-12*x^4-6*x*y+25*y^2-24*x^2*y^2-12*y^4
In??[??]:
implicit_plot3d(z-f(x,y),(x,-0.6,0.6),(y,-0.6,0.6),(z,-0.5,1.5),
       color="orange",plot_points=200,opacity=0.6)
In??[136]:
zmin=-2
zmax=10
h=0.4
nc=round((zmax-zmin)/h)
cont=contour_plot(f(x,y),(x,-0.8,0.8),(y,-0.6,0.6),
                  fill=False,linewidths=1, contours=[zmin+ i*h for i in range(nc+1)],cmap='cool')
show(cont,figsize=4)
No description has been provided for this image

User defined function??

In??[137]:
import numpy as np
def random_search(g,alpha_choice,max_iter,w,num_samples):
    xpoints = [] # container for approximate points
    fvals = []   # Container for function values at approximate points.
    alpha = 0
    n = np.size(w)
    for k in range(1,max_iter+1):
        if alpha_choice == 'diminishing':
            alpha = 1/float(k)
        else:
            alpha = alpha_choice

        xpoints.append(w)
        fvals.append(g(w))
        dirs = np.random.randn(num_samples,n) 
        norms = sqrt(np.sum(dirs*dirs,axis = 1))[:,np.newaxis]
        dirs = dirs/norms
        w_new = w+alpha*dirs
        evals = np.array([g(p) for p in w_new])
        min_indx = np.argmin(evals)
        if g(w_new[min_indx]) < g(w):
            d = dirs[min_indx,:]
            w = w + alpha*d
    xpoints.append(w)
    fvals.append(g(w))
    return xpoints, fvals
In??[138]:
g = lambda w: 25*w[0]**2-12*w[0]**4-6*w[0]*w[1]+25*w[1]**2-24*w[0]**2*w[1]**2-12*w[1]**4
import numpy as np
# run random search algorithm 
# alpha_choice = 'diminishing'; 
alpha_choice = 0.1; 
w = np.array([-0.6,-0.5]); num_samples = 100; max_iter = 20;
xpoints,fvals = random_search(g,alpha_choice,max_iter,w,num_samples)
xpoints[-1],fvals[-1]
Out[138]:
(array([0.0136398, 0.0118859]), 0.00720895511436233)
In??[139]:
pts = point2d(xpoints,color='red')
directions = []
for i in range(max_iter):
    L = line([(xpoints[i][0],xpoints[i][1]),(xpoints[i+1][0],xpoints[i+1][1])])
    directions.append(L)
show(cont+pts+sum(directions),figsize=5)
No description has been provided for this image

Optimization using Nelder-Mead Method??

In??[140]:
## module downhill
''' x = downhill(F,xStart,side,tol=1.0e-6)
    Downhill simplex method for minimizing the user-supplied
    scalar function F(x) with respect to the vector x.
    xStart = starting vector x.
    side   = side length of the starting simplex (default is 0.1)
'''
from numpy import zeros,dot,argmax,argmin,sum
from math import sqrt

def downhill(F,xStart,side=0.1,tol=1.0e-6, maxit=500):
    n = len(xStart)                 # Number of variables
    x = zeros((n+1,n)) 
    f = zeros(n+1)
    
  # Generate starting simplex by adding a small number to each coordinate
    x[0] = xStart
    for i in range(1,n+1):
        x[i] = xStart
        x[i,i-1] = xStart[i-1] + side # increament i the cordinate by side       
  # Compute values of F at the vertices of the simplex     
    for i in range(n+1): f[i] = F(x[i])
    
  # Main loop
    for k in range(maxit):
      # Find highest and lowest vertices
        iLo = argmin(f)
        iHi = argmax(f)       
      # Compute the move vector d
        d = (-(n+1)*x[iHi] + sum(x,axis=0))/n
      # Check for convergence
        if sqrt(dot(d,d)/n) < tol: return x[iLo]
        
      # Try reflection
        xNew = x[iHi] + 2.0*d              
        fNew = F(xNew)        
        if fNew <= f[iLo]:        # Accept reflection 
            x[iHi] = xNew
            f[iHi] = fNew
          # Try expanding the reflection
            xNew = x[iHi] + d               
            fNew = F(xNew)
            if fNew <= f[iLo]:    # Accept expansion
                x[iHi] = xNew
                f[iHi] = fNew
        else:
          # Try reflection again
            if fNew <= f[iHi]:    # Accept reflection
                x[iHi] = xNew
                f[iHi] = fNew
            else:
              # Try contraction
                xNew = x[iHi] + 0.5*d
                fNew = F(xNew)
                if fNew <= f[iHi]: # Accept contraction
                    x[iHi] = xNew
                    f[iHi] = fNew
                else:
                  # Use shrinkage
                    for i in range(len(x)):
                        if i != iLo:
                            x[i] = (x[i] - x[iLo])*0.5
                            f[i] = F(x[i])
    print( "Too many iterations in downhill")
    print ("Last values of x were")
    return( x[iLo])
In??[141]:
g = lambda w: 25*w[0]**2-12*w[0]**4-6*w[0]*w[1]+25*w[1]**2-24*w[0]**2*w[1]**2-12*w[1]**4
xStart=np.array([-0.6,-0.6])
downhill(g,xStart,side=0.1,tol=1.0e-6, maxit=100)
Out[141]:
array([8.84449886e-08, 5.04749149e-07])
In??[142]:
## Aimation
f(x,y)=25*x^2-12*x^4-6*x*y+25*y^2-24*x^2*y^2-12*y^4
zmin=-2
zmax=10
h=0.4
nc=round((zmax-zmin)/h)
cont=contour_plot(f(x,y),(x,-0.8,0.8),(y,-0.6,0.6),fill=False,linewidths=1, contours=[zmin+ i*h for i in range(nc+1)],cmap='cool')
v1=vector([-0.6,-0.6])
v2=vector([-0.4,-0.6])
v3=vector([-0.6,-0.4])
S=[v1,v2,v3]
n=len(S[0])
m=len(S)
simplex=point(S)+polygon(S,fill=False,color='black')
p=cont+simplex
P = [p]
for k in range(20):
    vf=[f(x=S[i][0],y=S[i][1]) for i in range(3)]
    from numpy import argmax,argmin
    iLo = argmin(vf)
    iHi = argmax(vf)   
    sv=vector([sum([S[i][j] for i in range(m)]) for j in range(n)])
    d = (-(n+1)*S[iHi] + sv)/n
    xNew = S[iHi] + 2.0*d              
    fNew = f(x=xNew[0],y=xNew[1])
    if fNew <= vf[iLo]:
        S[iHi] = xNew
        vf[iHi] = fNew
        xNew = S[iHi] + d           
        fNew = f(x=xNew[0],y=xNew[1])
        if fNew <= vf[iLo]: # Accept expansion
            S[iHi] = xNew
            vf[iHi] = fNew
    else: # Try reflection again
        if fNew <= vf[iHi]:     # Accept reflection
            S[iHi] = xNew
            vf[iHi] = fNew
        else:
            xNew = S[iHi] + 0.5*d
            fNew = f(x=xNew[0],y=xNew[1])
            if fNew <= vf[iHi]: # Accept contraction
                 S[iHi] = xNew
                 vf[iHi] = fNew
            else: # Use shrinkage
                for i in range(3):
                      if i != iLo:
                          S[i] = (S[i] - S[iLo])*0.5
                          vf[i] = f(x=S[i][0],y=S[i][1])
    p=p+point(S)+polygon(S,fill=False,color='black')
    P.append(p)
A = animate(P)
A.show(delay=100)
No description has been provided for this image

Solving Differential Equations??

Solve the folloing differential equation using Rk4-method. \begin{align} x' & = -x -2y\tag{1}\\ y' & = x + \frac{y}{5}\tag{2}\\ x(0) & = 1\tag{3}\\ y(0) & = -1\tag{4} \end{align}

In??[143]:
reset()
x, y, t = var('x y t')
F = [-x -2*y, x - y/5]
P = desolve_system_rk4(F, [x, y], ics = [0,-1,-1], ivar = t, end_points = 10, step = 0.01)
Q = [[j,k] for i,j,k in P]
p = line(Q, axes_labels=['$x(t)$','$y(t)$'], fontsize=12, thickness=2)
p += arrow(Q[int(len(Q)/5)], Q[int(len(Q)/5) + 1])
n = sqrt(F[0]^2 + F[1]^2)
F_unit = [F[0]/n, F[1]/n]
p += plot_vector_field(F_unit, (x,-1.5,1.5), (y,-1.5,1.5), color='red',
                       axes_labels=['$x(t)$','$y(t)$'])
p
Out[143]:
No description has been provided for this image

Modeling a Pandemic??

The SIR model is described by the nonlinear system of below, where is the transmission rate and is the recovery rate.

\begin{align*} s(t) & = \text{Susceptible individuals}\\ i(t) & = \text{Infected individuals}\\ r(t) & = \text{Removed individuals} \end{align*}

\begin{equation*} s(t) + i(t) + r(t) = N. \end{equation*}

\begin{align*} S(t) & = s(t)/N = \text{the fraction of susceptible individuals}\\ I(t) & = i(t)/N = \text{the fraction of infected individuals}\\ R(t) & = r(t)/N = \text{the fraction of removed individuals}\\ 1 & = S(t) + I(t) + R(t)\text{.} \end{align*}

\begin{align*} \frac{dS}{dt} & = - \alpha SI\\ \frac{dI}{dt} & = \alpha SI - \beta I\\ \frac{dR}{dt} & = \beta I. \end{align*}

In??[10]:
Susceptible, Infected, Recovered, t = var('Susceptible, Infected, Recovered, t')
alpha = 1.0
beta = 0.2
F = [-alpha*Susceptible*Infected, alpha*Susceptible*Infected - beta*Infected, beta*Infected]
P = desolve_system_rk4(F,[Susceptible, Infected, Recovered],ics=[0,0.99,0.01,0],ivar=t,end_points=30,step=0.1)
Q = [ [i, j] for i,j,k,l in P]
R = [ [i, k] for i,j,k,l in P]
S = [ [i, l] for i,j,k,l in P]
LQ = line(Q, color='red', axes_labels=['$t$','$S(t), I(t)$'], legend_label='$S(t)$',
          legend_color='red', fontsize=12, thickness=1)
LR = line(R, color='blue', legend_label='$I(t)$', legend_color='blue', thickness=1)
LS = line(S, color='green', legend_label='$R(t)$', legend_color='green', thickness=1)
p = LQ + LR + LS
p.show(figsize=5)
No description has been provided for this image

Using Sage Interact??

To be done in sagecell

In??[??]:
x,y = var('x,y')
from sage.ext.fast_eval import fast_float
@interact
def _(f = input_box(default=y^2), g=input_box(default=-x*y+x^3-x),
      xmin=input_box(default=-1), xmax=input_box(default=1),
      ymin=input_box(default=-1), ymax=input_box(default=1),
      start_x=input_box(default=-0.5), start_y=input_box(default=0.0),
      step_size=(0.01,(0.01, 0.2)), steps=(600,(0, 1400)) ):
    ff = fast_float(f, 'x', 'y')
    gg = fast_float(g, 'x', 'y')
    steps = int(steps)

    points = [ (start_x, start_y) ]
    for i in range(steps):
        xx, yy = points[-1]
        try:
            points.append( (xx + step_size * ff(xx,yy), yy + step_size * gg(xx,yy)) )
        except (ValueError, ArithmeticError, TypeError):
            break

    starting_point = point(points[0], pointsize=50)
    solution = line(points)
    F=[f,g]
    n = sqrt(F[0]^2 + F[1]^2)
    F_unit = [F[0]/n, F[1]/n]
    vector_field = plot_vector_field(F_unit, (x,xmin,xmax), (y,ymin,ymax) )
    Arr = arrow(points[int(len(points)/5)], points[int(len(points)/5) + 2])

    result = vector_field + starting_point + solution+Arr
    pretty_print(html(r"$\displaystyle\frac{dx}{dt} = %s$  $ \displaystyle\frac{dy}{dt} = %s$" % (latex(f),latex(g))))
    result.show(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax)

Ask Sage??

https://ask.sagemath.org/questions/

References??

  1. Sagemath documentations: https://doc.sagemath.org/
  2. Computational Mathematics with Sagemath by Paul Zimmerman
  3. Online book on Linear Algebra by Robert Beezer
  4. Abstract Algebra: Theory and Applications by by Tom Judson, Sage material by Rob Beezer

    https://judsonbooks.org/aata/

  5. Series of YouTube video lectures (Twelve) on SageMath by Ajit Kumar

    https://ajitmathsoft.wordpress.com/sagemath/

  6. NPTEL course on Computational Matheatics with Sagemath by Ajit Kumar
  7. Group Threory: A expedition with SageMath by Ajit Kumar and Vikas Bist (2021)
  8. Learning and Experiencing Cryptography with CrypTool and SageMath by Bernhard Esslinger (2024)
  9. The Ordinary Differential Equations Project by Thomas W. Judson (2024)

    https://judsonbooks.org/odeproject/odeproject-html/frontmatter.html