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

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

Leave a comment here

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s