Sage: Eigenstuff

This help sheet is very useful!

Suppose we want to find the eigenvalues and eigenvectors of the matrix
A = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}
We can find the eigenvalues using

A = matrix([[2,1],[1,2]])
A.eigenvalues()

A list of eigenvalues, with corresponding eigenvectors and multiplicities, is obtained from

A = matrix([[2,1],[1,2]])
A.eigenvectors_right()

The code

A = matrix([[2,1],[1,2]])
A.eigenmatrix_right()

produces two matrices: the first is a diagonal matrix with eigenvalues down the diagonal; the second is the transition matrix where the column vectors are eigenvectors. One way to use this function is

A = matrix([[2,1],[1,2]])
D,P = A.eigenmatrix_right()
show(P)

which designates the diagonal matrix to be D and the transition matrix to be P.

Posted in Sage, Uncategorized | Leave a comment

Sage: Numerically solving differential equations

Euler’s method by hand

We can easily make Sage do the computations in Euler’s method.
Here I explain how to build code to analyze the initial value problem
\frac{dy}{dt} = y^2
y(0) =1
First we define variable y and also the function f(y) = y^2 that determines the right hand side:

var('y')
f(y) = y^2

Then we define \Delta t and set it equal to 0.1.
We also define y_0, which we set equal to 1:

var('y')
f(y) = y^2

y0 = 1
deltaT = 0.1

Next we construct a data structure that will hold the various values of t_k and y_k.
To do this we make a variable called steps that tells us how many time steps we will take.
In this example, we set steps equal to 5.
We call that data structure eulerData.
For now, we define all the t_k values to be t_k = k(\Delta t) and all the y_k values to be the same as y_0.

var('y')
f(y) = y^2

y0 = 1
deltaT = 0.1
steps=5

eulerData = [[k*deltaT,y0] for k in range(0,steps+1)]

eulerData

The variable eulerData is organized as follows: the quantity eulerData[2] tells us the entries in the row corresponding to k=2.

var('y')
f(y) = y^2

y0 = 1
deltaT = 0.1
steps=5

eulerData = [[k*deltaT,y0] for k in range(0,steps+1)]

eulerData[2]

If we want only to see the time t_2, then we need to ask Sage to show us eulerData[2][0], while if we want Sage to show us y_2, we need to ask to be shown eulerData[2][1]. Try it:

var('y')
f(y) = y^2

y0 = 1
deltaT = 0.1
steps=5

eulerData = [[k*deltaT,y0] for k in range(0,steps+1)]

eulerData[2][0]

What we want to do now is systematically go through and update all the y_k values, starting at k=1 and ending at k=\texttt{steps}.
We can do this with a loop:

var('y')
f(y) = y^2

y0 = 1
deltaT = 0.1
steps=5

eulerData = [[k*deltaT,y0] for k in range(0,steps+1)]

for k in [1..steps]:
    eulerData[k][1]=eulerData[k-1][1]+deltaT*f(eulerData[k-1][1])
    
eulerData

Finally, we can plot the resulting approximate solution:

var('y')
f(y) = y^2

y0 = 1
deltaT = 0.1
steps=5

eulerData = [[k*deltaT,y0] for k in range(0,steps+1)]

for k in [1..steps]:
    eulerData[k][1]= eulerData[k-1][1] + deltaT*f(eulerData[k-1][1])
    
eulerPlot = list_plot(eulerData, color="red", plotjoined=true,marker='.',axes_labels=['$t$','$y$'])
eulerPlot.show()

Runge-Kutta

(to appear)

Runge-Kutta for systems

First we declare variables t, y_1, y_2.
We also define the function F that determines the right hand side, set the initial condition to be (t_0, y_1(t_0), y_2(t_0)) = (0,0.1, 0.1), and tell the computer that we we want to run the simulation until time t=20:

var('t,y1,y2')

F = [y1- y1^2 - y1*y2,-y2+2*y1*y2]
initCond = [0,.1, .1]
endTime = 20

Next we construct a variable numSoln which is the numerical solution.

var('t,y1,y2')

F = [y1- y1^2 - y1*y2,-y2+2*y1*y2]
initCond = [0,.1, .1]
endTime = 20

numSoln = desolve_system_rk4(F,[y1,y2],ics=initCond, ivar=t,end_points=endTime)

The object numSoln consists of triples [i,j,k] that represent (t_i, y_{1,i}, y_{2,i}).
In order to construct a parametric plot, we only need the y_1 and y_2 values. These we store in an object called parList.
We subsequently make a parametric plot from that list.

var('t,y1,y2')

F = [y1- y1^2 - y1*y2,-y2+2*y1*y2]
initCond = [0,.1, .1]
endTime = 20

numSoln = desolve_system_rk4(F,[y1,y2],ics=initCond, ivar=t,end_points=endTime)

parList = [ [j,k] for i,j,k in numSoln]
parPlot = list_plot(parList, color='red', plotjoined=true,thickness=2)

Next we want to plot a vector field as before.
To do this we need to construct a vector version of F. We do this and then build the plot.

var('t,y1,y2')

F = [y1- y1^2 - y1*y2,-y2+2*y1*y2]
initCond = [0,.1, .1]
endTime = 20

numSoln = desolve_system_rk4(F,[y1,y2],ics=initCond, ivar=t,end_points=endTime)

parList = [ [j,k] for i,j,k in numSoln]
parPlot = list_plot(parList, color='red', plotjoined=true,thickness=2)

vectorF = vector(F)
normalF = vectorF/vectorF.norm()
vfPlot = plot_vector_field(normalF,(y1,-.2,1.2),(y2,-.2,1.2),axes_labels=['$y_1$','$y_2$'])

mainPlot = vfPlot + parPlot
mainPlot.show()

A fun example

For fun, let’s numerically solve the Van der Pol oscillator, which is described by
\frac{dx}{dt} = y
\frac{dy}{dt} = \mu (1-x^2)y - x
for some parameter \mu. For this example, we set \mu=1.

If we run this code

var('t,x,y')
mu=1

Field = vector([y,mu*(1-x^2)*y-x])
InitialCondition = [0,.1,.1]
EndTime=20

NumSoln = desolve_system_rk4(Field, [x,y], ics=InitialCondition, ivar=t, end_points=EndTime)
ParPlot = list_plot([[j,k] for i,j,k in NumSoln], plotjoined=true)

FieldPlot = plot_vector_field(Field/Field.norm(),(x,-3,3),(y,-4,4))

MainPlot = FieldPlot + ParPlot
MainPlot.show()

we see the following picture:

van-der-pol

Posted in Sage, Uncategorized | Leave a comment

Sage: Parametric curves and vector fields

This is one in a series of posts containing short snippets of Sage code for my students to use: see my page on Sage.

Vector fields can be plotted using the vector_field command:

var('y1,y2')

vfield = vector([y1+y2,y1-y2])

plot_vector_field(vfield,(y1,-5,5),(y2,-5,5))

Sometimes it is convenient to plot the normalized version of the vector field:

var('y1,y2')

vfield = vector([y1+y2,y1-y2])
nfield = vfield/vfield.norm()

plot_vector_field(vfield,(y1,-5,5),(y2,-5,5))

Parametric curves can be plotted using the parametric_plot function:

var('t')

y1(t) = 2*exp(-.1*t)*cos(t)
y2(t) = 2*exp(-.1*t)*sin(t)

parametric_plot([y1(t),y2(t)], (t,0,5*pi))

Here’s a fun example that has a (normalized) vector field as well as a parametric curve that follows the vector field.

var('t,y1,y2')

vfield = vector([y1+y2,y1-y2])
nfield = vfield/vfield.norm()

nplot = plot_vector_field(nfield,(y1,-2,5),(y2,-2,5))

y1(t) = (1+sqrt(2))*exp(sqrt(2)*t) + (1-sqrt(2))*exp(-sqrt(2)*t)
y2(t) = exp(sqrt(2)*t) + exp(-sqrt(2)*t)

pplot = parametric_plot([y1(t),y2(t)], (t,-1,.5))

nplot + pplot

Related to vector fields are the slope field plots used in the differential equations course. Suppose we have the differential equation

\frac{dy}{dt} = y(1-y).

The corresponding slope field is given by:

var('t','y')
plot_slope_field(y*(1-y), (t,0,2), (y,-1,2))

Here is a fun example:

var('t','y')
plot_slope_field(y*(1-y)*sin(t), (t,0,10), (y,-1,2))
Posted in Sage, Uncategorized | Leave a comment

Sage: plotting functions

This is one in a series of posts containing short snippets of Sage code for my students to use: see my page on Sage.

Plotting a single function

Simple plotting of functions can be done with the plot command:

var('t')
f(t) = t^2
plot(f,-1,2)

A slightly more complicated piece of code that yields the same result is

var('t')
f(t) = t^2
plot(f(t),(t,-1,2))

We can restrict the vertical range of the plot with the code

var('t')
f(t) = t^2
plot(f(t),(t,-1,2), ymin=-.5, ymax = 3)

It is possible to add all sort of features to the plot: labels, legends, etc. Note that the labels can be done in LaTeX!

var('t')
f(t) = t^2
plot(f(t),(t,-1,2),ymin=-.5,ymax = 3, axes_labels=['$t$','gallons'], legend_label='$f(t)=t^2$',show_legend=true)

We can also play around with the style of the plot itself:

var('t')
f(t) = t^2
plot(f(t),(t,-1,2),ymin=-.5,ymax = 3, thickness=2, color='purple', linestyle='dashed')

For a whole list of options for plotting, see the Sage plotting page.

Finally, Sage can handle asymptotes… if we tell it to. First try out this code:

var('t')
f(t) = 1/((t-1)*(t+2))
plot(f(t),(t,-3,2),ymin=-10,ymax = 10)

The following, using detect_poles is much better:

var('t')
f(t) = 1/((t-1)*(t+2))
plot(f(t),(t,-3,2),ymin=-10,ymax = 10, detect_poles='show')

Plotting multiple functions

In order to plot multiple functions at the same time, we first take advantage of the feature that we can assign names to plots. When we do this, we need to use the .show() command to actually show the plot. We can incorporate many of the display options in to the .show() command:

var('t')
f(t) = t^2
fplot = plot(f(t), (t,-10,10))
fplot.show(xmin = -2, xmax=3, ymin=-1, ymax=5, axes_labels=['$t$',' '])

With these tools in mind, we can easily add the plot of a second function. Take a close look at each line of code here; what does each accomplish?

var('t')
f(t) = t^2
g(t) = 4-t^2

fplot = plot(f(t), (t,-10,10), color='purple')
gplot = plot(g(t), (t,-10,10), color='blue')

mainplot = fplot + gplot

mainplot.show(xmin = -2, xmax=3, ymin=-1, ymax=5, axes_labels=['$t$',' '])

Here is an example of something a bit fancier (and kind of fun)

var('t')
f(t) = t^2 - t^3
g = derivative(f,t)

fplot = plot(f(t), (t,-10,10), color='purple', legend_label='$f(t)$')
gplot = plot(g(t), (t,-10,10), color='blue', legend_label='$f^\prime(t)$')

mainplot = fplot + gplot

mainplot.show(xmin = -1, xmax=2, ymin=-1, ymax=1, axes_labels=['$t$',' '], show_legend=true)

Saving plots

We can save a plot as a pdf file using the command .save(). Here is a simple example:

var('t')
f(t) = t^2 - t^3

fplot = plot(f(t), (t,-1,2), ymax=1, ymin=-1,color='purple')

fplot.save(filename='fun-plot.pdf')

If I am going to be using the plot in a LaTeX document, I like to adjust the size. For example:

var('t')
f(t) = t^2 - t^3

fplot = plot(f(t), (t,-1,2), ymax=1, ymin=-1,color='purple')

fplot.save(filename='fun-plot.pdf', figsize=[4,3])

Of course, you can also add in all the other fancy stuff while saving the plot to a pdf!

var('t')
f(t) = t^2 - t^3

fplot = plot(f(t), (t,-1,2),color='purple',legend_label='$f(t)$')

fplot.save(filename='fun-plot.pdf', axes_labels=['$t$', ' '], show_legend=true,ymax=1, ymin=-1,figsize=[4,3])
Posted in Sage, Uncategorized | Leave a comment

Sage: Basics – computation, functions, calculus 1

This is one in a series of posts containing short snippets of Sage code for my students to use: see my page on Sage.

Basic computations can be done by simply entering in the expression you want to evaluate. For example,

2+3/5 - 5/9

If you want to display the result more nicely, you might try

show(2+3/5 - 5/9)

If you want to express the result as a decimal approximation, use

n(2+3/5 - 5/9)

If you want to obtain the LaTeX code for the result, use

latex(2+3/5 - 5/9)

Sometimes you it helps to use the command

simplify( )

Defining functions

In order to define functions, we need to have independent variables. Sage recognizes “x” as an independent variable. For other variables, we need to tell Sage that we are using them. For example:

var('t')
f(t) = t^2
f(4)

Here’s another example:

var('t','y')
f(t,y) = t*t - cos(t*y)
f(1,pi)

Notice three things about this previous example: (1) the cosine function is built in, (2) the constant pi is built in, and (3) we need to explicitly tell Sage to multiply — Sage does not interpret juxtaposition as multiplication.

It is also possible to define anonymous functions, though I rarely use this.

g = lambda t: t^2
g(7/3)

Basic calculus

Derivatives are rather straightforward:

var('t')
p(t) = t^2
q = derivative(p,t)
q(t)

Anti-derivatives do not include the arbitrary constant:

var('t')
p(t) = t^2
q = integral(p,t)
q(t)

In order to compute a definite integral, simply put in the desired interval:

var('t')
p(t) = t^2
integral(p(t),(t,2,pi))

If you want a numerical approximation, use the numerical_integral command:

var('t')
p(t) = exp(t^2)
numerical_integral(p,2,pi)

Notice that the output of the numerical integral has two pieces: the first is the approximate value; the second is an estimate of the error. If you only want the estimate, use this:

var('t')
p(t) = exp(t^2)
numerical_integral(p,2,pi)[0]
Posted in Sage | Leave a comment

First steps with Sagemath

I am spending some time this summer playing around with Sagemath. The hope is that I can learn it well enough to use it in my differential equations course in the near future. In this post I list some resources that I am finding helpful. (The list is to be updated as I continue.)

Getting started

  • I have decided to work primarily in the SageMathCloud (SMC) environment. One advantage of this is that I don’t have to worry about whether certain files are on my home or office machines. The other is that there seem to be some interesting options for assigning/grading/returning projects to students in that setting. I’ve opted for the “basic” $7 per month plan.
  • Students will have the option of working in SMC, installing Sage on their personal machines, or using the SageMathCell. The SageMathCell seems very handy, as you can get right to the computing without a lot of fuss. (Of course, work there is not saved, etc.) Having students working on their own machines is fine… but there does not seem to be a graceful way to move work between SMC Worksheet and the Sage Notebook formats. Thus far, this is my biggest complaint with the Sage framework.
  • I also downloaded Sagemath to my laptop and to my office machine.
  • Actually getting started with SMC was more difficult than I had hoped. There is a lot of powerful functionality, but it is not so obvious where a beginner like me should begin. After thrashing around for a day or so, I ended up doing most of my work in Sage Worksheets, where I am able to write a lot of notes to myself about what I am doing. It is nice to be able to populate the notebook with markdown cells and to be able to use LaTeX when commenting on the code.

General resources

Resources for teaching with Sage

Posted in Sage, Uncategorized | Leave a comment

New Paper: Boundary Value Problems and Finite Differences

My article Boundary Value Problems and Finite Differences has just been appeared in the January 2016 issue of the College Mathematics Journal. Here is the abstract:

The solvability of boundary value problems differs greatly from that of initial value problems and can be somewhat difficult to make sense of in the context of a sophomore-level differential equations course. We present an approach that uses finite difference approximations to motivate and understand the theory governing the existence and uniqueness of solutions to boundary value problems at an elementary level.

Posted in Fun & Interesting, Uncategorized | Leave a comment

Using Overleaf for student papers

In my introductory differential equations course, I assign students to write several papers. I require the students to typeset their papers with LaTeX, and to use graphics imported from Mathematica.

This semester, the paper assignments are:

  • Equilibrium solutions for the logistic model. This first assignment is mathematically extremely simple. The purpose is to give the students experience with LaTeX typesetting, with generating plots using Mathematica, and with paper writing in a technical setting.
  • Comparing harvesting models. In this assignment student compare/contrast two different ways to modify the logistic model in order to account for harvesting. The purpose is to for students to think more deeply about the assumptions that go in to the models and about the predictions that the models make.
  • Predator-prey models. In this assignment students are explore a simple predator-prey system with variable coupling. The purpose of this assignment is for students to discover the bifurcation points of the system, and to interpret the bifurcation in the context of population modeling.
  • Gravitation. In this assignment, students study a small body orbiting a larger body according to Newtonian mechanics. The purpose is for the students to learn how to apply, and interpret, conserved quantities.

In all of these assignments, I want students to become proficient with LaTeX and to develop technical writing skills. I also want the students to become familiar with the process of using multiple pieces of software in order to create a single product.

One of the challenges associated to these assignments is getting LaTeX typesetting software up and running on whatever device students are using. The Mac labs on our campus all have the excellent (and free) TeXShop program installed, but it was a bit of a hassle getting  appropriate software installed on the students’ personal machines.

Fortunately, I discovered Overleaf, which allows students to write and compile LaTeX via a web browser, and to store their papers on a server in the cloud. Thus far, the advantages of using Overleaf are

  • All students, regardless of which machine or operating system they use, are using the same software.
  • Students can access their projects from any machine.
  • I can use Overleaf to share a template (requires this graphic file) for the students to use.
  • Students can sign up for Overleaf using their LC Google ID.

The main disadvantage of Overleaf seems to be that processing speed is sometimes slow. This is due to the fact that the compiling is happening remotely. Thus far, it has not been a large issue.

This semester I devoted one full class day to showing students how to use Overleaf, and giving them a quick introduction to LaTeX. Students also worked through a self-paced module that introduced them to Mathematica. This took place in the computer lab, and seemed to work very well.

A second day was spent discussing elements of writing that are specific to math. This included a discussion of document structure, as well as conventions (not starting sentences with symbols, making all mathematics part of a sentence, etc.).

These two days, together with the detailed feedback I gave on the first assignment, seemed to be sufficient to get most students up and running with LaTeX and mathematical writing.

Posted in Differential equations, LaTeX, Teaching, Writing | Leave a comment

Mathematica — Euler’s method for systems, revisited

An earlier post discusses Mathematica code for Euler’s method. Here is some updated code.

A couple preliminary notes:

  • I like to clear all variables at the beginning.
  • For systems it is fun to superimpose the plot on top of the stream plot of the corresponding vector field

Example 1
For the IVP \frac{dy}{dt} = t-y, y(0) = 2 with \Delta t = 0.1 we can use the code

Clear["Global`*"]
f[t_, y_] := t - y
deltat = 0.1;
eulerValues = 
RecurrenceTable[{t[k + 1] == t[k] + deltat,
   y[k + 1] == y[k] + deltat*f[t[k],y[k]], 
   t[0] == 0., 
   y[0] == 2.}, 
   {t, y}, {k, 0, 10}];
Grid[eulerValues]
eulerPlot=ListLinePlot[eulerValues, PlotMarkers -> Automatic]
realPlot = Plot[3 Exp[-t] + t - 1, 
	{t, 0, 1}, PlotStyle -> Red];
Show[eulerPlot, realPlot]

Example 2
Consider \frac{dy}{dt} = e^y, y(0)=1, \Delta t = .25
We use

Clear["Global`*"]
deltat = .05;
f[y_] := Exp[y];
values = RecurrenceTable[{
    t[k + 1] == t[k] + deltat,
    y[k + 1] == y[k] + deltat*f[y[k]],
    t[0] == 0, y[0] == 1
    }, {t, y}, {k, 0, 10}];
Grid[values]
ListLinePlot[values, PlotMarkers -> Automatic]

Example 3
Here we consider the system
\frac{dx}{dt} = x(1-x) - xy
\frac{dy}{dt} = -y + 2 xy
with \Delta t = 0.1 and initial conditions x(0) = 1.5, y(0) = 1. The following code processes the first 100 time steps.

Clear["Global`*"]

f[x_, y_] := x (1 - x) - x*y
g[x_, y_] := -y + 2 x*y
streamPlot = StreamPlot[
	{f[x, y], g[x, y]}, 
	{x, -1, 2}, 
	{y, -1, 2}]
	
deltat = 0.1;
steps = 100;

eulerValues = RecurrenceTable[{
    t[0] == 0, x[0] == 1.5, y[0] == 1,
    t[k + 1] == t[k] + deltat, 
    x[k + 1] == x[k] + deltat*f[x[k], y[k]], 
    y[k + 1] == y[k] + deltat*g[x[k], y[k]]},
   {t, x, y},
   {k, 0, steps}];
eulerPlotValues = Table[{
    eulerValues[[k, 2]],
    eulerValues[[k, 3]]},
   {k, 1, steps + 1}];
eulerPlot = ListLinePlot[
  eulerPlotValues,
  PlotMarkers -> Automatic,
  PlotStyle -> Red]
  
Show[streamPlot, eulerPlot]
Posted in Differential equations, Mathematica | Leave a comment

Mathematica for linear algebra

Here is some mathematica code to be used by my linear algebra students. The post will be periodically updated as the course progresses.

Entering matrices

Mathematica views matrices as lists of columns. Thus to enter the matrix
\left(  \begin{array}{cccc}   1 & 1 & -1 & -1 \\   1 & 2 & 3 & 4 \\   1 & 3 & 5 & 7 \\  \end{array}  \right)
we type in

mat = {{1, 1, -1, -1}, {1, 2, 3, 4}, {1, 3, 5, 7}};

and hit [enter]. Here we have given the matrix the name mat and suppressed the output.

To view the matrix, type in the code

MatrixForm[mat]

and hit [enter].

Row reducing a matrix

To row reduce the matrix, type in the code

RowReduce[mat]

and hit [enter]. Notice that the output is just given in list format. To display the result as a matrix, you can use either

MatrixForm[RowReduce[mat]]

or

RowReduce[mat] // MatrixForm
Posted in Linear Algebra, Mathematica | Leave a comment