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:

van-der-pol