## Sage code for trigonometric forcing

var('t')

omegaf = 3
limit=130

A(t) = 1-(1/(4-omegaf^2))*cos((omegaf-2)*t)
h(t) = cos(2*t)

AmpPlot = plot(A(t), (t,0,limit),linestyle='dashed', color='red',thickness=2)
hPlot = plot(h(t)*A(t), (t,0,limit))

mainPlot = AmpPlot+ hPlot
mainPlot.show(ymin = -5,ymax = 5, axes_labels=['$t$',' '],ticks = [[],[]])<span id="mce_SELREST_start" style="overflow:hidden;line-height:0;">﻿</span>


## The teachable student

A recent picture of Yeovil manager Darren Way shows him with a poster listing the characteristics of “The Coachable Player.”

The Coachable Player

• Humble
• Respectful
• Loves the game
• Stays in control
• Takes responsibility
• Thinks long term
• Keen to learn
• Excited by change
• Willing to try new things
• Unafraid of mistakes
• Trusts coaches

It strikes me that if you replace “coachable player” with “teachable student” then the list remains essentially unchanged.

Other things that should be on the list of characteristics of “the teachable student”?

## 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$.

## 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

We can use the built-in Runge-Kutta algorithm to generate numerical solutions.
Consider the following code, which provides a numerical solution to
$\frac{dy}{dt} = y^2$
on the time interval $0\leq t\leq 1$ with time steps of size $\Delta t = 0.25$.

var('y','t')

f(y) = y^2   # rhs of ODE

t0, y0 = 0, 1 # initial conditions

numsoln = desolve_rk4( f(y), y, ics=[t0,y0], ivar=t, end_points=1, step=0.25)

numsoln


To generate a list plot, we use the code

var('y','t')

f(y) = y^2   # rhs of ODE

t0, y0 = 0, 1 # initial conditions

numsoln = desolve_rk4( f(y), y, ics=[t0,y0], ivar=t, end_points=1, step=0.25)

numplot = list_plot(numsoln, marker=".",plotjoined=true)
numplot.show(figsize=[4,3])


### 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:

## 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 | 1 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.

### 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])


## 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]


## 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

## 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.

## 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.