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:

This entry was posted in Sage, Uncategorized. Bookmark the permalink.