PREP Quickstart, Differential Equations
system:sage


<h1 style="font-size: 2em;">Sage Quickstart for Differential Equations</h1>
<p>This&nbsp;<a href="http://www.sagemath.org" target="_blank">Sage</a>&nbsp;worksheet was developed for the MAA PREP Workshop "Sage: Using Open-Source Mathematics Software with Undergraduates" (funding provided by NSF DUE 0817071).</p>
<p>Solving differential equations is a combination of exact and numerical methods, and hence a great place to explore with the computer. &nbsp;We have already seen one example of this in the calculus tutorial, which is worth reviewing.</p>

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

<p>Forgetting about plotting for the moment, notice that there are three things one needs to solve a differential equation symbolically:</p>
<ul>
<li>an abstract function $y$;&nbsp;</li>
<li>a differential equation, which here we put in a separate line;</li>
<li>the actual Differential Equation SOLVE command.</li>
</ul>

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

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

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

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

{{{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
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

<p>If you wanted to make $y$ an abstract function again instead of a variable, you'd have to do that. &nbsp;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.</p>
<p>We can't cover all the variants of this in a Quickstart, but the documentation is good for symbolic solvers.</p>

{{{id=12|
desolvers?
///
<html><!--notruncate-->
<div class="docstring">
    
  <p><strong>File:</strong> /usr/local/sage-prep/local/lib/python2.6/site-packages/sage/calculus/desolvers.py</p>
<p><strong>Type:</strong> &lt;type &#8216;module&#8217;&gt;</p>
<p><strong>Definition:</strong> desolvers( [noargspec] )</p>
<p><strong>Docstring:</strong></p>
<p>Solving ordinary differential equations</p>
<p>This file contains functions useful for solving differential equations
which occur commonly in a 1st semester differential equations
course. For another numerical solver see <tt class="xref docutils literal"><span class="pre">ode_solver()</span></tt> function
and optional package Octave.</p>
<p>Commands:</p>
<ul class="simple">
<li><tt class="docutils literal"><span class="pre">desolve</span></tt> - Computes the &#8220;general solution&#8221; to a 1st or 2nd order
ODE via Maxima.</li>
<li><tt class="docutils literal"><span class="pre">desolve_laplace</span></tt> - Solves an ODE using laplace transforms via
Maxima. Initials conditions are optional.</li>
<li><tt class="docutils literal"><span class="pre">desolve_system</span></tt> - Solves any size system of 1st order odes using
Maxima. Initials conditions are optional.</li>
<li><tt class="docutils literal"><span class="pre">desolve_rk4</span></tt> - Solves numerically IVP for one first order
equation, returns list of points or plot</li>
<li><tt class="docutils literal"><span class="pre">desolve_system_rk4</span></tt> - Solves numerically IVP for system of first
order equations, returns list of points</li>
<li><tt class="docutils literal"><span class="pre">eulers_method</span></tt> - Approximate solution to a 1st order DE,
presented as a table.</li>
<li><tt class="docutils literal"><span class="pre">eulers_method_2x2</span></tt> - Approximate solution to a 1st order system
of DEs, presented as a table.</li>
<li><tt class="docutils literal"><span class="pre">eulers_method_2x2_plot</span></tt> - Plots the sequence of points obtained
from Euler&#8217;s method.</li>
</ul>
<p>AUTHORS:</p>
<ul class="simple">
<li>David Joyner (3-2006) - Initial version of functions</li>
<li>Marshall Hampton (7-2007) - Creation of Python module and testing</li>
<li>Robert Bradshaw (10-2008) - Some interface cleanup.</li>
<li>Robert Marik (10-2009) - Some bugfixes and enhancements</li>
</ul>


</div>
</html>
}}}

<p>For instance, Maxima can do systems, as well as use Laplace transforms. &nbsp;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&rsquo;(0)=2$, which gives the last list in the example.</p>

{{{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
}}}

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

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

<p>It can be fun to compare this with the original, symbolic solution. &nbsp;Notice the use of the 'points' command from the advanced plotting tutorial.</p>

{{{id=20|
h1 = desolve(de, y, ics=[0,3])
plot(h1,(x,0,5),color='red')+points(h)
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

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

{{{id=22|
ode_solver?
///
<html><!--notruncate-->
<div class="docstring">
    
  <p><strong>File:</strong> /usr/local/sage-prep/devel/sage/sage/gsl/ode.pyx</p>
<p><strong>Type:</strong> &lt;type &#8216;type&#8217;&gt;</p>
<p><strong>Definition:</strong> ode_solver( [noargspec] )</p>
<p><strong>Docstring:</strong></p>
<blockquote>
<p><tt class="xref docutils literal"><span class="pre">ode_solver()</span></tt> is a class that wraps the GSL libraries ode
solver routines To use it instantiate a class,:</p>
<div class="highlight-python"><div class="highlight"><pre class="literal-block"><span class="gp">sage: </span><span class="n">T</span><span class="o">=</span><span class="n">ode_solver</span><span class="p">()</span>
</pre></div>
</div>
<p>To solve a system of the form <tt class="docutils literal"><span class="pre">dy_i/dt=f_i(t,y)</span></tt>, you must
supply a vector or tuple/list valued function <tt class="docutils literal"><span class="pre">f</span></tt> representing
<tt class="docutils literal"><span class="pre">f_i</span></tt>.  The functions <tt class="docutils literal"><span class="pre">f</span></tt> and the jacobian should have the
form <tt class="docutils literal"><span class="pre">foo(t,y)</span></tt> or <tt class="docutils literal"><span class="pre">foo(t,y,params)</span></tt>.  <tt class="docutils literal"><span class="pre">params</span></tt> which is
optional allows for your function to depend on one or a tuple of
parameters.  Note if you use it, <tt class="docutils literal"><span class="pre">params</span></tt> must be a tuple even
if it only has one component.  For example if you wanted to solve
<span class="math">y''+y=0</span>. You need to write it as a first order system:</p>
<div class="highlight-python"><pre class="literal-block">y_0' = y_1
y_1' = -y_0</pre>
</div>
<p>In code:</p>
<div class="highlight-python"><div class="highlight"><pre class="literal-block"><span class="gp">sage: </span><span class="n">f</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">t</span><span class="p">,</span><span class="n">y</span><span class="p">:[</span><span class="n">y</span><span class="p">[</span><span class="mf">1</span><span class="p">],</span><span class="o">-</span><span class="n">y</span><span class="p">[</span><span class="mf">0</span><span class="p">]]</span>
<span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">function</span><span class="o">=</span><span class="n">f</span>
</pre></div>
</div>
<p>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 <tt class="docutils literal"><span class="pre">[</span> <span class="pre">[df_1/dy_1,...,df_1/dy_n],</span> <span class="pre">...,</span>
<span class="pre">[df_n/dy_1,...,df_n,dy_n],</span> <span class="pre">[df_1/dt,...,df_n/dt]</span> <span class="pre">]</span></tt>.</p>
<p>There are examples below, if your jacobian was the function
<tt class="docutils literal"><span class="pre">my_jacobian</span></tt> you would do:</p>
<div class="highlight-python"><div class="highlight"><pre class="literal-block"><span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">jacobian</span> <span class="o">=</span> <span class="n">my_jacobian</span>     <span class="c"># not tested, since it doesn&#39;t make sense to test this</span>
</pre></div>
</div>
<p>There are a variety of algorithms available for different types of systems. Possible algorithms are</p>
<ul class="simple">
<li><tt class="docutils literal"><span class="pre">rkf45</span></tt> - runga-kutta-felhberg (4,5)</li>
<li><tt class="docutils literal"><span class="pre">rk2</span></tt> - embedded runga-kutta (2,3)</li>
<li><tt class="docutils literal"><span class="pre">rk4</span></tt> - 4th order classical runga-kutta</li>
<li><tt class="docutils literal"><span class="pre">rk8pd</span></tt> - runga-kutta prince-dormand (8,9)</li>
<li><tt class="docutils literal"><span class="pre">rk2imp</span></tt> - implicit 2nd order runga-kutta at gaussian points</li>
<li><tt class="docutils literal"><span class="pre">rk4imp</span></tt> - implicit 4th order runga-kutta at gaussian points</li>
<li><tt class="docutils literal"><span class="pre">bsimp</span></tt> - implicit burlisch-stoer (requires jacobian)</li>
<li><tt class="docutils literal"><span class="pre">gear1</span></tt> - M=1 implicit gear</li>
<li><tt class="docutils literal"><span class="pre">gear2</span></tt> - M=2 implicit gear</li>
</ul>
<p>The default algorithm is <tt class="docutils literal"><span class="pre">rkf45</span></tt>. If you instead wanted to use
<tt class="docutils literal"><span class="pre">bsimp</span></tt> you would do:</p>
<div class="highlight-python"><div class="highlight"><pre class="literal-block"><span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">algorithm</span><span class="o">=</span><span class="s">&quot;bsimp&quot;</span>
</pre></div>
</div>
<p>The user should supply initial conditions in y_0. For example if
your initial conditions are y_0=1,y_1=1, do:</p>
<div class="highlight-python"><div class="highlight"><pre class="literal-block"><span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">y_0</span><span class="o">=</span><span class="p">[</span><span class="mf">1</span><span class="p">,</span><span class="mf">1</span><span class="p">]</span>
</pre></div>
</div>
<p>The actual solver is invoked by the method <tt class="xref docutils literal"><span class="pre">ode_solve()</span></tt>.  It
has arguments <tt class="docutils literal"><span class="pre">t_span</span></tt>, <tt class="docutils literal"><span class="pre">y_0</span></tt>, <tt class="docutils literal"><span class="pre">num_points</span></tt>, <tt class="docutils literal"><span class="pre">params</span></tt>.
<tt class="docutils literal"><span class="pre">y_0</span></tt> 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 <tt class="docutils literal"><span class="pre">ode_solve</span></tt> or by
assignment.</p>
<p><tt class="docutils literal"><span class="pre">t_span</span></tt> is the time interval on which to solve the ode.  If
<tt class="docutils literal"><span class="pre">t_span</span></tt> 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 <tt class="docutils literal"><span class="pre">t_span[0]</span></tt> and <tt class="docutils literal"><span class="pre">t_span[1]</span></tt>. If
<tt class="docutils literal"><span class="pre">t_span</span></tt> is a tuple with more than two values than the values of
y_i at points in time listed in <tt class="docutils literal"><span class="pre">t_span</span></tt> will be returned.</p>
<p>Error is estimated via the expression <tt class="docutils literal"><span class="pre">D_i</span> <span class="pre">=</span>
<span class="pre">error_abs*s_i+error_rel*(a|y_i|+a_dydt*h*|y_i'|)</span></tt>.  The user can
specify <tt class="docutils literal"><span class="pre">error_abs</span></tt> (1e-10 by default), <tt class="docutils literal"><span class="pre">error_rel</span></tt> (1e-10 by
default) <tt class="docutils literal"><span class="pre">a</span></tt> (1 by default), <tt class="docutils literal"><span class="pre">a_(dydt)</span></tt> (0 by default) and
<tt class="docutils literal"><span class="pre">s_i</span></tt> (as scaling_abs which should be a tuple and is 1 in all
components by default).  If you specify one of <tt class="docutils literal"><span class="pre">a</span></tt> or <tt class="docutils literal"><span class="pre">a_dydt</span></tt>
you must specify the other.  You may specify <tt class="docutils literal"><span class="pre">a</span></tt> and <tt class="docutils literal"><span class="pre">a_dydt</span></tt>
without <tt class="docutils literal"><span class="pre">scaling_abs</span></tt> (which will be taken =1 be default).
<tt class="docutils literal"><span class="pre">h</span></tt> is the initial step size which is (1e-2) by default.</p>
<p><tt class="docutils literal"><span class="pre">ode_solve</span></tt> solves the solution as a list of tuples of the form,
<tt class="docutils literal"><span class="pre">[</span> <span class="pre">(t_0,[y_1,...,y_n]),(t_1,[y_1,...,y_n]),...,(t_n,[y_1,...,y_n])]</span></tt>.</p>
<p>This data is stored in the variable solutions:</p>
<div class="highlight-python"><div class="highlight"><pre class="literal-block"><span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">solution</span>               <span class="c"># not tested</span>
</pre></div>
</div>
<p>EXAMPLES:</p>
<p>Consider solving the Van der Pol oscillator <span class="math">x''(t) +
ux'(t)(x(t)^2-1)+x(t)=0</span> between <span class="math">t=0</span> and <span class="math">t= 100</span>.  As a first
order system it is <span class="math">x'=y</span>, <span class="math">y'=-x+uy(1-x^2)</span>. Let us take <span class="math">u=10</span>
and use initial conditions <span class="math">(x,y)=(1,0)</span> and use the runga-kutta
prince-dormand algorithm.</p>
<div class="highlight-python"><div class="highlight"><pre class="literal-block"><span class="gp">sage: </span><span class="k">def</span> <span class="nf">f_1</span><span class="p">(</span><span class="n">t</span><span class="p">,</span><span class="n">y</span><span class="p">,</span><span class="n">params</span><span class="p">):</span>
<span class="gp">... </span>     <span class="k">return</span><span class="p">[</span><span class="n">y</span><span class="p">[</span><span class="mf">1</span><span class="p">],</span><span class="o">-</span><span class="n">y</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span><span class="o">-</span><span class="n">params</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span><span class="o">*</span><span class="n">y</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span><span class="o">*</span><span class="p">(</span><span class="n">y</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span><span class="o">**</span><span class="mf">2</span><span class="o">-</span><span class="mf">1.0</span><span class="p">)]</span>

<span class="gp">sage: </span><span class="k">def</span> <span class="nf">j_1</span><span class="p">(</span><span class="n">t</span><span class="p">,</span><span class="n">y</span><span class="p">,</span><span class="n">params</span><span class="p">):</span>
<span class="gp">... </span>     <span class="k">return</span> <span class="p">[</span> <span class="p">[</span><span class="mf">0.0</span><span class="p">,</span> <span class="mf">1.0</span><span class="p">],[</span><span class="o">-</span><span class="mf">2.0</span><span class="o">*</span><span class="n">params</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span><span class="o">*</span><span class="n">y</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span><span class="o">*</span><span class="n">y</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span><span class="o">-</span><span class="mf">1.0</span><span class="p">,</span><span class="o">-</span><span class="n">params</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span><span class="o">*</span><span class="p">(</span><span class="n">y</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span><span class="o">*</span><span class="n">y</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span><span class="o">-</span><span class="mf">1.0</span><span class="p">)],</span> <span class="p">[</span><span class="mf">0.0</span><span class="p">,</span><span class="mf">0.0</span><span class="p">]</span> <span class="p">]</span>

<span class="gp">sage: </span><span class="n">T</span><span class="o">=</span><span class="n">ode_solver</span><span class="p">()</span>
<span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">algorithm</span><span class="o">=</span><span class="s">&quot;rk8pd&quot;</span>
<span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">function</span><span class="o">=</span><span class="n">f_1</span>
<span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">jacobian</span><span class="o">=</span><span class="n">j_1</span>
<span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">ode_solve</span><span class="p">(</span><span class="n">y_0</span><span class="o">=</span><span class="p">[</span><span class="mf">1</span><span class="p">,</span><span class="mf">0</span><span class="p">],</span><span class="n">t_span</span><span class="o">=</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span><span class="mf">100</span><span class="p">],</span><span class="n">params</span><span class="o">=</span><span class="p">[</span><span class="mf">10.0</span><span class="p">],</span><span class="n">num_points</span><span class="o">=</span><span class="mf">1000</span><span class="p">)</span>
<span class="gp">sage: </span><span class="n">outfile</span> <span class="o">=</span> <span class="n">SAGE_TMP</span> <span class="o">+</span> <span class="s">&#39;sage.png&#39;</span>
<span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">plot_solution</span><span class="p">(</span><span class="n">filename</span><span class="o">=</span><span class="n">outfile</span><span class="p">)</span>
</pre></div>
</div>
<p>The solver line is equivalent to:</p>
<div class="highlight-python"><div class="highlight"><pre class="literal-block"><span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">ode_solve</span><span class="p">(</span><span class="n">y_0</span><span class="o">=</span><span class="p">[</span><span class="mf">1</span><span class="p">,</span><span class="mf">0</span><span class="p">],</span><span class="n">t_span</span><span class="o">=</span><span class="p">[</span><span class="n">x</span><span class="o">/</span><span class="mf">10.0</span> <span class="k">for</span> <span class="n">x</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mf">1000</span><span class="p">)],</span><span class="n">params</span> <span class="o">=</span> <span class="p">[</span><span class="mf">10.0</span><span class="p">])</span>
</pre></div>
</div>
<p>Let&#8217;s try a system:</p>
<div class="highlight-python"><pre class="literal-block">y_0'=y_1*y_2
y_1'=-y_0*y_2
y_2'=-.51*y_0*y_1</pre>
</div>
<p>We will not use the jacobian this time and will change the
error tolerances.</p>
<div class="highlight-python"><div class="highlight"><pre class="literal-block"><span class="gp">sage: </span><span class="n">g_1</span><span class="o">=</span> <span class="k">lambda</span> <span class="n">t</span><span class="p">,</span><span class="n">y</span><span class="p">:</span> <span class="p">[</span><span class="n">y</span><span class="p">[</span><span class="mf">1</span><span class="p">]</span><span class="o">*</span><span class="n">y</span><span class="p">[</span><span class="mf">2</span><span class="p">],</span><span class="o">-</span><span class="n">y</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span><span class="o">*</span><span class="n">y</span><span class="p">[</span><span class="mf">2</span><span class="p">],</span><span class="o">-</span><span class="mf">0.51</span><span class="o">*</span><span class="n">y</span><span class="p">[</span><span class="mf">0</span><span class="p">]</span><span class="o">*</span><span class="n">y</span><span class="p">[</span><span class="mf">1</span><span class="p">]]</span>
<span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">function</span><span class="o">=</span><span class="n">g_1</span>
<span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">y_0</span><span class="o">=</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span><span class="mf">1</span><span class="p">,</span><span class="mf">1</span><span class="p">]</span>
<span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">scale_abs</span><span class="o">=</span><span class="p">[</span><span class="mf">1e-4</span><span class="p">,</span><span class="mf">1e-4</span><span class="p">,</span><span class="mf">1e-5</span><span class="p">]</span>
<span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">error_rel</span><span class="o">=</span><span class="mf">1e-4</span>
<span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">ode_solve</span><span class="p">(</span><span class="n">t_span</span><span class="o">=</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span><span class="mf">12</span><span class="p">],</span><span class="n">num_points</span><span class="o">=</span><span class="mf">100</span><span class="p">)</span>
</pre></div>
</div>
<p>By default T.plot_solution() plots the y_0, to plot general y_i use:</p>
<div class="highlight-python"><div class="highlight"><pre class="literal-block"><span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">plot_solution</span><span class="p">(</span><span class="n">i</span><span class="o">=</span><span class="mf">0</span><span class="p">,</span> <span class="n">filename</span><span class="o">=</span><span class="n">outfile</span><span class="p">)</span>
<span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">plot_solution</span><span class="p">(</span><span class="n">i</span><span class="o">=</span><span class="mf">1</span><span class="p">,</span> <span class="n">filename</span><span class="o">=</span><span class="n">outfile</span><span class="p">)</span>
<span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">plot_solution</span><span class="p">(</span><span class="n">i</span><span class="o">=</span><span class="mf">2</span><span class="p">,</span> <span class="n">filename</span><span class="o">=</span><span class="n">outfile</span><span class="p">)</span>
</pre></div>
</div>
<p>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.</p>
<div class="highlight-python"><div class="highlight"><pre class="literal-block"><span class="gp">sage: </span><span class="n">f</span> <span class="o">=</span> <span class="n">T</span><span class="o">.</span><span class="n">interpolate_solution</span><span class="p">()</span>
<span class="gp">sage: </span><span class="n">plot</span><span class="p">(</span><span class="n">f</span><span class="p">,</span><span class="mf">0</span><span class="p">,</span><span class="mf">12</span><span class="p">)</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
<span class="gp">sage: </span><span class="n">f</span> <span class="o">=</span> <span class="n">T</span><span class="o">.</span><span class="n">interpolate_solution</span><span class="p">(</span><span class="n">i</span><span class="o">=</span><span class="mf">1</span><span class="p">)</span>
<span class="gp">sage: </span><span class="n">plot</span><span class="p">(</span><span class="n">f</span><span class="p">,</span><span class="mf">0</span><span class="p">,</span><span class="mf">12</span><span class="p">)</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
<span class="gp">sage: </span><span class="n">f</span> <span class="o">=</span> <span class="n">T</span><span class="o">.</span><span class="n">interpolate_solution</span><span class="p">(</span><span class="n">i</span><span class="o">=</span><span class="mf">2</span><span class="p">)</span>
<span class="gp">sage: </span><span class="n">plot</span><span class="p">(</span><span class="n">f</span><span class="p">,</span><span class="mf">0</span><span class="p">,</span><span class="mf">12</span><span class="p">)</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
<span class="gp">sage: </span><span class="n">f</span> <span class="o">=</span> <span class="n">T</span><span class="o">.</span><span class="n">interpolate_solution</span><span class="p">()</span>
<span class="gp">sage: </span><span class="n">f</span><span class="p">(</span><span class="n">pi</span><span class="p">)</span>
<span class="go">0.5379...</span>
</pre></div>
</div>
<p>The solver attributes may also be set up using arguments to
ode_solver.  The previous example can be rewritten as:</p>
<div class="highlight-python"><div class="highlight"><pre class="literal-block"><span class="gp">sage: </span><span class="n">T</span> <span class="o">=</span> <span class="n">ode_solver</span><span class="p">(</span><span class="n">g_1</span><span class="p">,</span><span class="n">y_0</span><span class="o">=</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span><span class="mf">1</span><span class="p">,</span><span class="mf">1</span><span class="p">],</span><span class="n">scale_abs</span><span class="o">=</span><span class="p">[</span><span class="mf">1e-4</span><span class="p">,</span><span class="mf">1e-4</span><span class="p">,</span><span class="mf">1e-5</span><span class="p">],</span><span class="n">error_rel</span><span class="o">=</span><span class="mf">1e-4</span><span class="p">,</span> <span class="n">algorithm</span><span class="o">=</span><span class="s">&quot;rk8pd&quot;</span><span class="p">)</span>
<span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">ode_solve</span><span class="p">(</span><span class="n">t_span</span><span class="o">=</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span><span class="mf">12</span><span class="p">],</span><span class="n">num_points</span><span class="o">=</span><span class="mf">100</span><span class="p">)</span>
<span class="gp">sage: </span><span class="n">f</span> <span class="o">=</span> <span class="n">T</span><span class="o">.</span><span class="n">interpolate_solution</span><span class="p">()</span>
<span class="gp">sage: </span><span class="n">f</span><span class="p">(</span><span class="n">pi</span><span class="p">)</span>
<span class="go">0.5379...</span>
</pre></div>
</div>
<p>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:</p>
<div class="highlight-python"><pre class="literal-block">%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</pre>
</div>
<p>After executing the above block of code you can do the
following (WARNING: the following is <em>not</em> automatically
doctested):</p>
<div class="highlight-python"><div class="highlight"><pre class="literal-block"><span class="gp">sage: </span><span class="n">T</span> <span class="o">=</span> <span class="n">ode_solver</span><span class="p">()</span>                     <span class="c"># not tested</span>
<span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">algorithm</span> <span class="o">=</span> <span class="s">&quot;bsimp&quot;</span>                <span class="c"># not tested</span>
<span class="gp">sage: </span><span class="n">vander</span> <span class="o">=</span> <span class="n">van_der_pol</span><span class="p">()</span>               <span class="c"># not tested</span>
<span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">function</span><span class="o">=</span><span class="n">vander</span>                    <span class="c"># not tested</span>
<span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">ode_solve</span><span class="p">(</span><span class="n">y_0</span> <span class="o">=</span> <span class="p">[</span><span class="mf">1</span><span class="p">,</span><span class="mf">0</span><span class="p">],</span> <span class="n">t_span</span><span class="o">=</span><span class="p">[</span><span class="mf">0</span><span class="p">,</span><span class="mf">2000</span><span class="p">],</span> <span class="n">num_points</span><span class="o">=</span><span class="mf">1000</span><span class="p">)</span>   <span class="c"># not tested</span>
<span class="gp">sage: </span><span class="n">T</span><span class="o">.</span><span class="n">plot_solution</span><span class="p">(</span><span class="n">i</span><span class="o">=</span><span class="mf">0</span><span class="p">,</span> <span class="n">filename</span><span class="o">=</span><span class="n">SAGE_TMP</span> <span class="o">+</span> <span class="s">&#39;/test.png&#39;</span><span class="p">)</span>        <span class="c"># not tested</span>
</pre></div>
</div>
</blockquote>


</div>
</html>
}}}

<p>We can even do power series solutions. &nbsp;We do have to keep in mind that we must then define a special ring, and always define the precision.</p>

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

<p>This power series solution is pretty good for a while!</p>

{{{id=24|
plot(h,(t,-2,5))+plot(2+e^-x,(x,-2,5),color='red',linestyle='--')
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

<p>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. &nbsp;</p>

{{{id=6|

///
}}}