Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download
Project: MAT 3770
Views: 3752
Kernel: SageMath 8.3

Notes and Tutorials

This page represents a collection of notes and tutorials for common functions, methods, definitions, etc. used in our class. This is not comprehensive and is intended as an aid for the course. For additional help, use the SAGE cheat sheet, SAGE tutorial, or stackoverflow.


Defining functions using SAGE

There are many ways to define functions in SAGE. The two main methods are:

  • defining a function explictly in terms of variables or

  • defining a variable and then using it in a function formula.

# first method f(x) = x^2+3
show(f)
x  x2+3\renewcommand{\Bold}[1]{\mathbf{#1}}x \ {\mapsto}\ x^{2} + 3
# second method y = var('y') g = sin(y)
show(g)
sin(y)\renewcommand{\Bold}[1]{\mathbf{#1}}\sin\left(y\right)

Note that the show command simply produces a nicely formated output using LaTeX.\LaTeX.

Differentiating functions in SAGE

There are several ways to differentiate a function using SAGE. If there are multiple variables, you will have to explicitly state the variable you are differentiating with respect to.

# as a method f.diff()
x |--> 2*x
# as a function derivative(f)
x |--> 2*x
# for multiple variables h(x,y) = x^2+y h.diff(y)
(x, y) |--> 1

Solving an equation in SAGE

The function solve has the format: solve(list of desired equations, variables to solve for {comma separated, not a list}).

# single variable example solve(3*x+1==0,x)
[x == (-1/3)]
# multiple variable example solve([x==y,x^2+y^2==1],x,y)
[[x == -1/2*sqrt(2), y == -1/2*sqrt(2)], [x == 1/2*sqrt(2), y == 1/2*sqrt(2)]]

The output of the solve command is a list where each element is a solution to the set of equations. Any expression of the form variable == function/formula/number is a symbolic expression object in SAGE. The right hand side of the equation given by this object can be extracted using the rhs method. (See the example below.)

# the sol varible below is the list of solutions sol = solve([x==y,x^2+y^2==1],x,y)
show(sol)
[[x=122,y=122],[x=122,y=122]]\renewcommand{\Bold}[1]{\mathbf{#1}}\left[\left[x = -\frac{1}{2} \, \sqrt{2}, y = -\frac{1}{2} \, \sqrt{2}\right], \left[x = \frac{1}{2} \, \sqrt{2}, y = \frac{1}{2} \, \sqrt{2}\right]\right]

Our system of equations has two solutions. We will extract the yy values of the first solution. (Note that you do not need to use all of the code below to accomplish this task. Here we show every intermediate step so as not to confuse.)

# remember that python/SAGE uses zero as the first index firstSol = sol[0]
show(firstSol)
[x=122,y=122]\renewcommand{\Bold}[1]{\mathbf{#1}}\left[x = -\frac{1}{2} \, \sqrt{2}, y = -\frac{1}{2} \, \sqrt{2}\right]
yEquation = firstSol[1]
show(yEquation)
y=122\renewcommand{\Bold}[1]{\mathbf{#1}}y = -\frac{1}{2} \, \sqrt{2}
yEquation.rhs()
-1/2*sqrt(2)
# single line version of the above work show(sol[0][1].rhs())
122\renewcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{2} \, \sqrt{2}

How to compute the sensitivity of yy with respect to xx

The sensitivity of yy with respect to xx is a measure for the percentage change in yy given a percentage change in xx. It is a unitless quantity defined by S(y,x)=dydxxy.S(y,x)=\dfrac{dy}{dx}\cdot\dfrac{x}{y}. If S(y,x)S(y,x) is close to 0, yy is not sensitive with respect to xx. A negative S(y,x)S(y,x) means that yy and xx have an inverse relationship (as xx goes up, yy goes down and vice versa).

# Note that we use X and Y instead of x and y to avoid problems with the definition of x and y in previous lines Y(X) = X^2+sin(X)
S(X) = Y.diff()*X/Y
# The method n is for numerical evaluations S(2).n()
1.46002690481524

The above shows that, when X=2X=2, if XX increases by 1%, then YY increases by 1.46%.

Graphing in SAGE

A plot in CoCalc is an object with various attributes. Graphics objects can be added together to produce multiple function plots on a single graph. Plots have many optional arguments that can add color, axes labels, and legends. Note that labels can be written in LaTeX.\LaTeX. Later on we will combine other types of plots together (like list_plot) and use the aspect_ratio optional argument to change how plots are displayed. These graphical objects are built on matplotlib, a module commonly used for plotting by data scientists that use Python.

p1 = plot(x^2,(x,-2,2),axes_labels=['$x$','$y$'],legend_label='$x^2$') p2 = plot(x,(x,-2,2),color='red',legend_label='$x$') p1+p2
Image in a Jupyter notebook

Finding Roots Numerically

The solve command cannot always find roots. The find_root method can find roots numerically, but an interval must be supplied. In the example below, we find a root of f(x)=x22f(x)=x^2-2.

eq = x^2-2 == 0 eq.find_root(0,2)
1.4142135623731364

Solving a System of Equations

The solve command can find solutions to systems of equations. We will find the solution to the following system of equations:

x+y=3x+y=3xy=5x-y=5
y = var('y') sol = solve([x+y==3,x-y==5],x,y)
sol
[[x == 4, y == -1]]
sol[0][0].rhs()
4
sol[0][1].rhs()
-1
contour_plot(x+sin(y),(x,-2,2),(y,-2,2))
Image in a Jupyter notebook
implicit_plot(x^2+y^2==1,(x,-2,2),(y,-2,2))
Image in a Jupyter notebook

3D Plotting

Plotting 3D surfaces in a Jupyter notebook can cause issues. There are two alternatives that you can use:

  1. Use a Sage worksheet to plot in 3D.

  2. Use the viewer='tachyon' option when plotting in 3D in a Jupyter notebook. Unfortunately, this option renders a static 3D image.

p = plot3d(x^2+sin(y),(x,-2,2),(y,-2,2)) show(p,viewer='tachyon')
Image in a Jupyter notebook

Finding Max/Mins of a Function in a Constrained Optimization Problem

Suppose you want to maximize/minimize f(x1,x2,,xn)f(x_1,x_2,\dots,x_n) subject to some constraints g1(x1,x2,,xn)c1,g2(x1,x2,,xn)c2,,gk(x1,x2,,xn)ckg_1(x_1,x_2,\dots,x_n)\leq c_1, g_2(x_1,x_2,\dots,x_n)\leq c_2, \dots, g_k(x_1,x_2,\dots,x_n)\leq c_k. Perform the following steps:

  1. Draw the feasible region (if possible).

  2. Look for internal critical points (f=0\nabla f=0).

  3. Use the method of Lagrange multipliers to find critical points on the boundary of the feasible region (f=λ1g1+λ2g2++λkgk\nabla f = \lambda_1\nabla g_1+\lambda_2\nabla g_2+\dots+\lambda_k\nabla g_k).

  4. Test all candidates from parts 2 and 3.

Monte Carlo Methods

The random() function in CoCalc produces uniformly distributed values on the interval [0,1)[0,1). By multiplying and adding the appropriate values, you can modify this function to produce numbers in any given interval. (Note: Other Python modules like numpy have more sophisticated function for generating random data.)

random()
0.17160677354490017

Newton's Method

Sometimes we are interested in finding the solution to the following system of equations:

f1(x1,x2,,xn)=0f2(x1,x2,,xn)=0fn(x1,x2,,xn)=0\begin{align*} f_1(x_1,x_2,\dots,x_n) &= 0\\ f_2(x_1,x_2,\dots,x_n) &= 0\\ \vdots \\ f_n(x_1,x_2,\dots,x_n) &= 0\\ \end{align*}

Let x=(x1,x2,,xn)\vec{x}=(x_1,x_2,\dots,x_n) and F(x)=(f1(x),f2(x),,fn(x))F(\vec{x})=(f_1(\vec{x}),f_2(\vec{x}),\dots,f_n(\vec{x})). We can use Newton's Method here to approximate the solution using an iterative procedure.

  1. Use a graphical method to determine a point x0\vec{x}_0 relatively close to the location of the true solution.

  2. Compute the Jacobian AA given by A:=A(x)=[df1dx1df1dxxdf1dxndf2dx1df2dxxdf2dxndfndx1dfndxxdfndxn]A:=A(\vec{x})=\begin{bmatrix}\frac{df_1}{dx_1} & \frac{df_1}{dx_x} & \cdots & \frac{df_1}{dx_n}\\ \frac{df_2}{dx_1} & \frac{df_2}{dx_x} & \cdots & \frac{df_2}{dx_n}\\ \vdots & \vdots & \ddots & \vdots \\\frac{df_n}{dx_1} & \frac{df_n}{dx_x} & \cdots & \frac{df_n}{dx_n}\end{bmatrix}

  3. Let xn+1=xnA1F(xn)\vec{x}_{n+1} = \vec{x}_n - A^{-1}F(\vec{x}_n) where A1A^{-1} is the inverse of AA. Note also that AA must be recomputed at each step.

Loading SAGE code from another file

Use the load() function to import SAGE code from a *.sage file into another file (like a Jupyter notebook).

load('/home/user/linProg.sage')

That is all for now.