2.4.3 Solving Differential Equations

You can use Sage to investigate ordinary differential equations. To solve the equation x'+x-1=0:

sage: t = var('t')    # define a variable t
sage: x = function('x',t)   # define x to be a function of that variable
sage: DE = lambda y: diff(y,t) + y - 1
sage: desolve(DE(x(t)), [x,t])
'%e^-t*(%e^t+%c)'
This uses Sage's interface to Maxima [Max], and so its output may be a bit different from other Sage output. In this case, this says that the general solution to the differential equation is x(t) = et(et + c).

You can compute Laplace transforms also; the Laplace transform of t2et - sin(t) is computed as follows:

sage: s = var("s")
sage: t = var("t")
sage: f = t^2*exp(t) - sin(t)
sage: f.laplace(t,s)
2/(s - 1)^3 - 1/(s^2 + 1)

Here is a more involved example. The displacement from equilibrium (respectively) for a coupled spring attached to a wall on the left

|------\/\/\/\/\---|mass1|----\/\/\/\/\/----|mass2|
         spring1               spring2
is modeled by the system of 2nd order differential equations

$\displaystyle m_1 x_1'' + (k_1+k_2) x_1 - k_2 x_2 = 0, \quad
m_2 x_2''+ k_2 (x_2-x_1) = 0,
$

where mi is the mass of object i, xi is the displacement from equilibrium of mass i, and ki is the spring constant for spring i.

Example: Use Sage to solve the above problem with m1=2, m2=1, k1=4, k2=3, x1(0)=3, x'1(0)=0, x2(0)=3, x'2(0)=0.

Solution: Take the Laplace transform of the first equation (with the notation x=x1, y=x2):

sage: de1 = maxima("2*diff(x(t),t, 2) + 6*x(t) - 2*y(t)")
sage: lde1 = de1.laplace("t","s"); lde1
2*(-?%at('diff(x(t),t,1),t=0)+s^2*?%laplace(x(t),t,s)-x(0)*s)-
2*?%laplace(y(t),t,s)+6*?%laplace(x(t),t,s)
This is hard to read, but it says that

$\displaystyle -2x'(0) + 2s^2*X(s) - 2sx(0) - 2Y(s) + 6X(s) = 0
$

(where the Laplace transform of a lower case function like x(t) is the upper case function X(s)). Take the Laplace transform of the second equation:
sage: de2 = maxima("diff(y(t),t, 2) + 2*y(t) - 2*x(t)")
sage: lde2 = de2.laplace("t","s"); lde2
-?%at('diff(y(t),t,1),t=0)+s^2*?%laplace(y(t),t,s)+2*?%laplace(y(t),t,s)-
2*?%laplace(x(t),t,s)-y(0)*s
This says

$\displaystyle -Y'(0) + s^2Y(s) + 2Y(s) - 2X(s) - sy(0) = 0.
$

Plug in the initial conditions for x(0), x'(0), y(0), and y'(0), and solve the resulting two equations:
sage: var('s X Y')
(s, X, Y)
sage: eqns = [(2*s^2+6)*X-2*Y == 6*s, -2*X +(s^2+2)*Y == 3*s] 
sage: solve(eqns, X,Y)
[[X == (3*s^3 + 9*s)/(s^4 + 5*s^2 + 4), 
  Y == (3*s^3 + 15*s)/(s^4 + 5*s^2 + 4)]]
Now take inverse Laplace transforms to get the answer:
sage: var('s t')
(s, t)
sage: inverse_laplace((3*s^3 + 9*s)/(s^4 + 5*s^2 + 4),s,t)
cos(2*t) + 2*cos(t)
sage: inverse_laplace((3*s^3 + 15*s)/(s^4 + 5*s^2 + 4),s,t)
4*cos(t) - cos(2*t)
Therefore, the solution is

$\displaystyle x_1(t) = \cos(2t) + 2\cos(t), \quad x_2(t) = 4\cos(t) - \cos(2t).
$

This can be plotted parametrically using
sage: t = var('t')
sage: P = parametric_plot((cos(2*t) + 2*cos(t), 4*cos(t) - cos(2*t) ),\
...   0, 2*pi, rgbcolor=hue(0.9))
sage: show(P)
The individual components can be plotted using
sage: t = var('t')
sage: p1 = plot(cos(2*t) + 2*cos(t), 0, 2*pi, rgbcolor=hue(0.3))
sage: p2 = plot(4*cos(t) - cos(2*t), 0, 2*pi, rgbcolor=hue(0.6))
sage: show(p1 + p2)
(For more on plotting, see Section 2.5.)

REFERENCES: Nagle, Saff, Snider, Fundamentals of Differential Equations, 6th ed, Addison-Wesley, 2004. (see §5.5).

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