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

First we define variable and also the function that determines the right hand side:

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

Then we define and set it equal to .

We also define , which we set equal to :

var('y')
f(y) = y^2
y0 = 1
deltaT = 0.1

Next we construct a data structure that will hold the various values of and .

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 .

We call that data structure **eulerData**.

For now, we define all the values to be and all the values to be the same as .

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 .

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 , then we need to ask Sage to show us **eulerData[2][0]**, while if we want Sage to show us , 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 values, starting at and ending at .

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

on the time interval with time steps of size .

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

We also define the function that determines the right hand side, set the initial condition to be , and tell the computer that we we want to run the simulation until time :

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 .

In order to construct a parametric plot, we only need the and 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

for some parameter . For this example, we set .

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: