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: