2.3 Ordinary differential equations

Symbolically solving ODEs can be done using Sage's interface with Maxima. Numerical solution of ODEs can be done using Sage's interface with Octave (an experimental package), or routines in the GSL (Gnu Scientific Library).

You can solve ODE's symbolically in Sage using the Maxima interface (do not type the ...):

sage: maxima.de_solve('diff(y,x,2) + 3*x = y', ['x','y'], [1,1,1])
y=3*x-2*%e^(x-1)
sage: maxima.de_solve('diff(y,x,2) + 3*x = y', ['x','y'])
y=%k1*%e^x+%k2*%e^-x+3*x
sage: maxima.de_solve('diff(y,x) + 3*x = y', ['x','y'])
y=(%c-3*(-x-1)*%e^-x)*%e^x
sage: maxima.de_solve('diff(y,x) + 3*x = y', ['x','y'],[1,1])
y=-%e^-1*(5*%e^x-3*%e*x-3*%e)

sage: maxima.clear('x'); maxima.clear('f')
sage: maxima.de_solve_laplace("diff(f(x),x,2) = 2*diff(f(x),x)-f(x)",\ 
...   ["x","f"], [0,1,2])
f(x)=x*%e^x+%e^x
            
sage: maxima.clear('x'); maxima.clear('f')            
sage: f = maxima.de_solve_laplace("diff(f(x),x,2) = 2*diff(f(x),x)-f(x)",\ 
...   ["x","f"])
sage: f
f(x)=x*%e^x*(?%at('diff(f(x),x,1),x=0))-f(0)*x*%e^x+f(0)*%e^x
sage: f.display2d()
                                               !
                                   x  d        !                  x          x
                        f(x) = x %e  (-- (f(x))!     ) - f(0) x %e  + f(0) %e
                                      dx       !
                                               !x = 0

If you have Octave and gnuplot installed,

sage: octave.de_system_plot(['x+y','x-y'], [1,-1], [0,2]) # optional octave required
yields the two plots $ (t,x(t)), (t,y(t))$ on the same graph (the $ t$ -axis is the horizonal axis) of the system of ODEs

$\displaystyle x' = x+y, x(0) = 1;\qquad y' = x-y, y(0) = -1,$   for$\displaystyle \quad 0 < t < 2.
$

Another way this system can be solved is to use the command desolve_system in calculus/examples.

sage: attach os.environ['SAGE_ROOT'] + '/examples/calculus/desolvers.sage'
sage: des = ["'diff(x(t),t)=x(t)+y(t)","'diff(y(t),t)=x(t)-y(t)"]
sage: vars = ["t","x","y"]
sage: ics = [0,1,-1]
sage: desolve_system(des,vars,ics)
[x(t)=cosh(2^(1/2)*t),y(t)=2*sinh(2^(1/2)*t)/sqrt(2)-cosh(2^(1/2)*t)]
The output of this command is not a pair of Sage functions.

Finally, Sage can solve linear DEs using power series:

sage: R.<t> = PowerSeriesRing(QQ, default_prec=10)
sage: a = 2 - 3*t + 4*t^2 + O(t^10)
sage: b = 3 - 4*t^2 + O(t^7)
sage: f = a.solve_linear_de(prec=5, b=b, f0=3/5)
sage: f
3/5 + 21/5*t + 33/10*t^2 - 38/15*t^3 + 11/24*t^4 + O(t^5)
sage: f.derivative() - a*f - b
O(t^4)

See About this document... for information on suggesting changes.