The Laplace transformation is a powerful tool to solve a vast class of ordinary differential equations. The principle is simple:
The mathematical theory of the Laplace transformation is not an
easy one, but for practical purposes it is often sufficient to simply
apply the transformation.
Murray R. Spiegel wrote a good tutorial about the topic, it was
published by McGraw-Hill under the title
Theory and Problems of Laplace Transforms.
The Laplace transformation maps a transcendental function onto a rational expression in a new variable:
laplace(exp(-b*t)*cos(w*t),t,s);
s + b -------------------- 2 2 2 w + s + 2 b s + b
assume(w > 0);
[w > 0]
Without that assumption, Maxima would ask later whether w is zero or nonzero.
The inverse Laplace transformation maps a rational expression onto a transcendental function:
ilt(%, s, t);
- b t %e cos(t w)
The knowledge of these functions is almost sufficient to solve an ordinary differential equation:
ode: 'diff(y(t), t, 2) + 5*'diff(y(t), t) + 4*y(t) = t;
It is important that you
2 d d --- (y(t)) + 5 (-- (y(t))) + 4 y(t) = t 2 dt dt
Next, we specify the initial conditions. We do that with the
atvalue(y(t), t=0, 0);
0
atvalue('diff(y(t), t), t= 0, 0);
0
Now we can compute the Laplace transformed of the equation:
lap_ode:laplace(ode,t,s);
2 1
s laplace(y(t), t, s) + 5 s laplace(y(t), t, s) + 4 laplace(y(t), t, s) = --
2
s
This is a linear equation in the unknown laplace(y(t), t, s). We solve it with solve:
sol: solve(%, 'laplace(y(t), t, s));
Note that you have to write the unknown with a quote. Without the quote, Maxima would try to evaluate the expression laplace(y(t), t, s).
1 [laplace(y(t), t, s) = ----------------] 4 3 2 s + 5 s + 4 s
The answer is a list with one equation, the solution of the linear equation.
Now, we can apply the inverse Laplace transformation. To apply apply a transformation to the elements of a list, we have to use map:
map( lambda( [eq], ilt(eq, s, t)), sol);
- t - 4 t
%e %e t 5
[y(t) = ----- - ------- + - - --]
3 48 4 16
This is again a list with one equation. Note that the inverse Laplace transformation was applied to both sides of the equation:
It is of course possible to solve this ordinary differential equation with the command ode2: Solution of the example with ode2
A system of ODE's
The following example is taken from Spiegel's book. It is a system of two ordinary differential equations
First we define the equations:
assume(s>0);
[s > 0]
diff_eq1: 'diff(f(x),x,2) + 'diff(g(x),x) + 3*f(x) = 15*exp(-x);
2
d d - x
-- (g(x)) + --- (f(x)) + 3 f(x) = 15 %e
dx 2
dx
diff_eq2: 'diff(g(x), x, 2) - 4*'diff(f(x), x) + 3*g(x) = 15*sin(2*x);
2
d d
--- (g(x)) - 4 (-- (f(x))) + 3 g(x) = 15 sin(2 x)
2 dx
dx
Next, we define the initial values. For this example, we have initial values for both functions and their derivatives at point x = 0.
atvalue (f(x), x=0, 35);
35
atvalue ('diff(f(x),x),x=0, -48);
- 48
atvalue (g(x), x=0, 27);
27
atvalue ('diff(g(x), x), x=0, -55);
- 55
Now we can compute the Laplace transforms:
lap_eq1:laplace(diff_eq1,x,s);
2 15
s laplace(g(x), x, s) + s laplace(f(x), x, s) + 3 laplace(f(x), x, s) - 35 s + 21 = -----
s + 1
lap_eq2:laplace(diff_eq2,x,s);
2 30
s laplace(g(x), x, s) + 3 laplace(g(x), x, s) - 4 (s laplace(f(x), x, s) - 35) - 27 s + 55 = ------
2
s + 4
These are two linear equations in the unknowns laplace(f(x),x, s) and laplace(g(x), x, s). We solve it with solve:
linsolve([lap_eq1,lap_eq2],['laplace(f(x),x,s),'laplace(g(x),x,s)]);
6 5 4 3 2
35 s - 13 s + 407 s + 185 s + 1020 s + 918 s - 72
[laplace(f(x), x, s) = ------------------------------------------------------,
7 6 5 4 3 2
s + s + 14 s + 14 s + 49 s + 49 s + 36 s + 36
6 5 4 3 2
27 s - 28 s + 50 s - 610 s - 787 s - 2022 s - 2250
laplace(g(x), x, s) = -------------------------------------------------------]
7 6 5 4 3 2
s + s + 14 s + 14 s + 49 s + 49 s + 36 s + 36
The denominators of the solutions are polynomials in s of degree 7. The computation of the inverse Laplace transform requires that the denominators are decomposed into linear and quadratic factors. The function ilt does not always factor denominators, so we do it:
factored: map(lambda( [eq], factor(eq)), %);
6 5 4 3 2
35 s - 13 s + 407 s + 185 s + 1020 s + 918 s - 72
[laplace(f(x), x, s) = ------------------------------------------------------,
2 2 2
(s + 1) (s + 1) (s + 4) (s + 9)
6 5 4 3 2
27 s - 28 s + 50 s - 610 s - 787 s - 2022 s - 2250
laplace(g(x), x, s) = -------------------------------------------------------]
2 2 2
(s + 1) (s + 1) (s + 4) (s + 9)
Now we can apply the inverse Laplace transform to both equations:
sol: map(lambda( [eq], ilt(eq, s, x)), factored);
We obtain a solution for the given system of ordinary differential equations:
- x
[f(x) = - 15 sin(3 x) + 2 cos(2 x) + 30 cos(x) + 3 %e ,
- x
g(x) = 30 cos(3 x) + sin(2 x) - 60 sin(x) - 3 %e ]
In principle we are done, but there is still one interesting technicality: How can we prove that the solution satisfies the given system of equations?
The proof requires four additional steps.
First, we write the two differential equations as a system:
ode_system: [diff_eq1, diff_eq2];
2
d d - x
[-- (g(x)) + --- (f(x)) + 3 f(x) = 15 %e ,
dx 2
dx
2
d d
--- (g(x)) - 4 (-- (f(x))) + 3 g(x) = 15 sin(2 x)]
2 dx
dx
Next we substitute the solution into that system:
ode_system, sol;
2
d - x
[--- (- 15 sin(3 x) + 2 cos(2 x) + 30 cos(x) + 3 %e )
2
dx
d - x
+ -- (30 cos(3 x) + sin(2 x) - 60 sin(x) - 3 %e )
dx
- x - x
+ 3 (- 15 sin(3 x) + 2 cos(2 x) + 30 cos(x) + 3 %e ) = 15 %e ,
d - x
- 4 (-- (- 15 sin(3 x) + 2 cos(2 x) + 30 cos(x) + 3 %e ))
dx
2
d - x
+ --- (30 cos(3 x) + sin(2 x) - 60 sin(x) - 3 %e )
2
dx
- x
+ 3 (30 cos(3 x) + sin(2 x) - 60 sin(x) - 3 %e ) = 15 sin(2 x)]
We obtain a list of two equations that contain symbolic derivatives. To enforce the computation of the derivatives, we have to evaluate the equations:
map (lambda( [eq], ev(eq, diff)), %);
- x [45 sin(3 x) + 3 (- 15 sin(3 x) + 2 cos(2 x) + 30 cos(x) + 3 %e ) - x - x - 6 cos(2 x) - 90 cos(x) + 6 %e = 15 %e , - x 3 (30 cos(3 x) + sin(2 x) - 60 sin(x) - 3 %e ) - 270 cos(3 x) - x - 4 (- 45 cos(3 x) - 4 sin(2 x) - 30 sin(x) - 3 %e ) - 4 sin(2 x) - x + 60 sin(x) - 3 %e = 15 sin(2 x)]
This looks better, but simplification is obviously needed:
trigsimp(%);
- x - x [15 %e = 15 %e , 15 sin(2 x) = 15 sin(2 x)]
These two identities complete the proof.
The next section shows that the equations of this example can also be solved with desolve.