Sage Quickstart for Differential Equations

This Sage worksheet was developed for the MAA PREP Workshop "Sage: Using Open-Source Mathematics Software with Undergraduates" (funding provided by NSF DUE 0817071).

Solving differential equations is a combination of exact and numerical methods, and hence a great place to explore with the computer.  We have already seen one example of this in the calculus tutorial, which is worth reviewing.

{{{id=1| y = function('y',x) de = diff(y,x) + y -2 h = desolve(de, y) /// }}}

Forgetting about plotting for the moment, notice that there are three things one needs to solve a differential equation symbolically:

{{{id=9| show(expand(h)) ///
\newcommand{\Bold}[1]{\mathbf{#1}}c e^{\left(-x\right)} + 2
}}}

Since we did not specify any initial conditions, Sage (from Maxima) puts in a parameter.  If we want to put in an initial condition, we use 'ics' (Initial ConditionS).

{{{id=13| h = desolve(de, y, ics=[0,3]); h /// (2*e^x + 1)*e^(-x) }}}

And of course we have already noted that we can plot all this with a slope field.

{{{id=15| var('y') # Needed so we can plot Plot1=plot_slope_field(2-y,(x,0,3),(y,0,5)) Plot2=plot(h,x,0,3) Plot1+Plot2 /// }}}

If you wanted to make $y$ an abstract function again instead of a variable, you'd have to do that.  Another option is to let $z$ be the name of the vertical axis variable, but either way something will have to give, since in common speaking about these things we treat $y$ as both a variable and a function, which is much trickier to accomplish with a computer.

We can't cover all the variants of this in a Quickstart, but the documentation is good for symbolic solvers.

{{{id=12| desolvers? ///

File: /usr/local/sage-prep/local/lib/python2.6/site-packages/sage/calculus/desolvers.py

Type: <type ‘module’>

Definition: desolvers( [noargspec] )

Docstring:

Solving ordinary differential equations

This file contains functions useful for solving differential equations which occur commonly in a 1st semester differential equations course. For another numerical solver see ode_solver() function and optional package Octave.

Commands:

AUTHORS:

}}}

For instance, Maxima can do systems, as well as use Laplace transforms.  It is important to pay attention to the syntax; in this example, we have placed the differential equation in the body of the command, and had to specify that $f$ was the Dependent VARiable, as well as give initial conditions $f(0)=1$ and $f’(0)=2$, which gives the last list in the example.

{{{id=11| f=function('f',x) desolve_laplace(diff(f,x,2) == 2*diff(f,x)-f, dvar = f, ics = [0,1,2]) /// x*e^x + e^x }}} {{{id=23| g(x)=x*e^x+e^x derivative(g,x,2)-2*derivative(g,x)+g /// x |--> 0 }}}

You may notice that one of the options above was 'desolve_rk4'.  This is a fourth-order Runge-Kutta method, and returns appropriate (numerical) output.  Here, we must give the dependent variable and initial conditions.

{{{id=3| y = function('y',x) de = diff(y,x) + y -2 h = desolve_rk4(de, y, step=.05, ics=[0,3]) /// }}}

It can be fun to compare this with the original, symbolic solution.  Notice the use of the 'points' command from the advanced plotting tutorial.

{{{id=20| h1 = desolve(de, y, ics=[0,3]) plot(h1,(x,0,5),color='red')+points(h) /// }}}

The primary use of numerical routines from here is pedagogical in nature.  There are also more advanced numerical routines available via the GNU scientific library.  Using this is a little more sophisticated, but gives a wealth of options.

{{{id=22| ode_solver? ///

File: /usr/local/sage-prep/devel/sage/sage/gsl/ode.pyx

Type: <type ‘type’>

Definition: ode_solver( [noargspec] )

Docstring:

ode_solver() is a class that wraps the GSL libraries ode solver routines To use it instantiate a class,:

sage: T=ode_solver()

To solve a system of the form dy_i/dt=f_i(t,y), you must supply a vector or tuple/list valued function f representing f_i. The functions f and the jacobian should have the form foo(t,y) or foo(t,y,params). params which is optional allows for your function to depend on one or a tuple of parameters. Note if you use it, params must be a tuple even if it only has one component. For example if you wanted to solve y''+y=0. You need to write it as a first order system:

y_0' = y_1
y_1' = -y_0

In code:

sage: f = lambda t,y:[y[1],-y[0]]
sage: T.function=f

For some algorithms the jacobian must be supplied as well, the form of this should be a function return a list of lists of the form [ [df_1/dy_1,...,df_1/dy_n], ..., [df_n/dy_1,...,df_n,dy_n], [df_1/dt,...,df_n/dt] ].

There are examples below, if your jacobian was the function my_jacobian you would do:

sage: T.jacobian = my_jacobian     # not tested, since it doesn't make sense to test this

There are a variety of algorithms available for different types of systems. Possible algorithms are

The default algorithm is rkf45. If you instead wanted to use bsimp you would do:

sage: T.algorithm="bsimp"

The user should supply initial conditions in y_0. For example if your initial conditions are y_0=1,y_1=1, do:

sage: T.y_0=[1,1]

The actual solver is invoked by the method ode_solve(). It has arguments t_span, y_0, num_points, params. y_0 must be supplied either as an argument or above by assignment. Params which are optional and only necessary if your system uses params can be supplied to ode_solve or by assignment.

t_span is the time interval on which to solve the ode. If t_span is a tuple with just two time values then the user must specify num_points, and the system will be evaluated at num_points equally spaced points between t_span[0] and t_span[1]. If t_span is a tuple with more than two values than the values of y_i at points in time listed in t_span will be returned.

Error is estimated via the expression D_i = error_abs*s_i+error_rel*(a|y_i|+a_dydt*h*|y_i'|). The user can specify error_abs (1e-10 by default), error_rel (1e-10 by default) a (1 by default), a_(dydt) (0 by default) and s_i (as scaling_abs which should be a tuple and is 1 in all components by default). If you specify one of a or a_dydt you must specify the other. You may specify a and a_dydt without scaling_abs (which will be taken =1 be default). h is the initial step size which is (1e-2) by default.

ode_solve solves the solution as a list of tuples of the form, [ (t_0,[y_1,...,y_n]),(t_1,[y_1,...,y_n]),...,(t_n,[y_1,...,y_n])].

This data is stored in the variable solutions:

sage: T.solution               # not tested

EXAMPLES:

Consider solving the Van der Pol oscillator x''(t) + ux'(t)(x(t)^2-1)+x(t)=0 between t=0 and t= 100. As a first order system it is x'=y, y'=-x+uy(1-x^2). Let us take u=10 and use initial conditions (x,y)=(1,0) and use the runga-kutta prince-dormand algorithm.

sage: def f_1(t,y,params):
...      return[y[1],-y[0]-params[0]*y[1]*(y[0]**2-1.0)]

sage: def j_1(t,y,params):
...      return [ [0.0, 1.0],[-2.0*params[0]*y[0]*y[1]-1.0,-params[0]*(y[0]*y[0]-1.0)], [0.0,0.0] ]

sage: T=ode_solver()
sage: T.algorithm="rk8pd"
sage: T.function=f_1
sage: T.jacobian=j_1
sage: T.ode_solve(y_0=[1,0],t_span=[0,100],params=[10.0],num_points=1000)
sage: outfile = SAGE_TMP + 'sage.png'
sage: T.plot_solution(filename=outfile)

The solver line is equivalent to:

sage: T.ode_solve(y_0=[1,0],t_span=[x/10.0 for x in range(1000)],params = [10.0])

Let’s try a system:

y_0'=y_1*y_2
y_1'=-y_0*y_2
y_2'=-.51*y_0*y_1

We will not use the jacobian this time and will change the error tolerances.

sage: g_1= lambda t,y: [y[1]*y[2],-y[0]*y[2],-0.51*y[0]*y[1]]
sage: T.function=g_1
sage: T.y_0=[0,1,1]
sage: T.scale_abs=[1e-4,1e-4,1e-5]
sage: T.error_rel=1e-4
sage: T.ode_solve(t_span=[0,12],num_points=100)

By default T.plot_solution() plots the y_0, to plot general y_i use:

sage: T.plot_solution(i=0, filename=outfile)
sage: T.plot_solution(i=1, filename=outfile)
sage: T.plot_solution(i=2, filename=outfile)

The method interpolate_solution will return a spline interpolation through the points found by the solver. By default y_0 is interpolated. You can interpolate y_i through the keyword argument i.

sage: f = T.interpolate_solution()
sage: plot(f,0,12).show()
sage: f = T.interpolate_solution(i=1)
sage: plot(f,0,12).show()
sage: f = T.interpolate_solution(i=2)
sage: plot(f,0,12).show()
sage: f = T.interpolate_solution()
sage: f(pi)
0.5379...

The solver attributes may also be set up using arguments to ode_solver. The previous example can be rewritten as:

sage: T = ode_solver(g_1,y_0=[0,1,1],scale_abs=[1e-4,1e-4,1e-5],error_rel=1e-4, algorithm="rk8pd")
sage: T.ode_solve(t_span=[0,12],num_points=100)
sage: f = T.interpolate_solution()
sage: f(pi)
0.5379...

Unfortunately because python functions are used, this solver is slow on system that require many function evaluations. It is possible to pass a compiled function by deriving from the class ode_sysem and overloading c_f and c_j with C functions that specify the system. The following will work in notebook:

%cython
cimport sage.gsl.ode
import sage.gsl.ode
include 'gsl.pxi'

cdef class van_der_pol(sage.gsl.ode.ode_system):
    cdef int c_f(self,double t, double *y,double *dydt):
        dydt[0]=y[1]
        dydt[1]=-y[0]-1000*y[1]*(y[0]*y[0]-1)
        return GSL_SUCCESS
    cdef int c_j(self, double t,double *y,double *dfdy,double *dfdt):
        dfdy[0]=0
        dfdy[1]=1.0
        dfdy[2]=-2.0*1000*y[0]*y[1]-1.0
        dfdy[3]=-1000*(y[0]*y[0]-1.0)
        dfdt[0]=0
        dfdt[1]=0
        return GSL_SUCCESS

After executing the above block of code you can do the following (WARNING: the following is not automatically doctested):

sage: T = ode_solver()                     # not tested
sage: T.algorithm = "bsimp"                # not tested
sage: vander = van_der_pol()               # not tested
sage: T.function=vander                    # not tested
sage: T.ode_solve(y_0 = [1,0], t_span=[0,2000], num_points=1000)   # not tested
sage: T.plot_solution(i=0, filename=SAGE_TMP + '/test.png')        # not tested
}}}

We can even do power series solutions.  We do have to keep in mind that we must then define a special ring, and always define the precision.

{{{id=4| R. = PowerSeriesRing(QQ, default_prec=10) a = -1 + 0*t b = 2 + 0*t h=a.solve_linear_de(b=b,f0=3,prec=10) /// }}}

This power series solution is pretty good for a while!

{{{id=24| plot(h,(t,-2,5))+plot(2+e^-x,(x,-2,5),color='red',linestyle='--') /// }}}

This was just an introduction; there are a lot of resources for differential equations using Sage, particularly from David Joyner, who wrote much of the original code for using Maxima in this way.  

{{{id=6| /// }}}