unit33.mws

ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Unit 33 -- Using Laplace Transforms to Solve Systems

Prof. Douglas B. Meade

Industrial Mathematics Institute

Department of Mathematics

University of South Carolina

Columbia, SC 29208

URL: http://www.math.sc.edu/~meade/

E-mail: meade@math.sc.edu

Copyright © 2001 by Douglas B. Meade

All rights reserved

-------------------------------------------------------------------

>

Outline of Unit 33

>

Initialization

> restart;

> with( DEtools ):

> with( plots ):

> with( linalg ):

> with( inttrans ):

Warning, the name changecoords has been redefined

Warning, the name adjoint has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name hilbert has been redefined

> alias( X(s) = laplace( x(t), t, s ),

> Y(s) = laplace( y(t), t, s ) ):

>

33.A Solving Systems of ODEs via the Laplace Transform

The same algorithm is applied when using Laplace transforms to solve a system of linear ODEs as for a single linear ODE. The only difference is that the transform of the system of ODEs is a system of algebraic equations. Once this system is solved for the transform of each unknown function, the inverse Laplace transform is applied to obtain the solution of the system of ODEs.

To illustrate, consider the system of linear ODEs with sinusoidal forcing of an unspecified frequency.

> sys1 := diff( x(t), t ) = y(t),

> diff( y(t), t ) = -4*x(t) + sin(omega*t);

sys1 := diff(x(t),t) = y(t), diff(y(t),t) = -4*x(t)...

>

The Laplace transform of the system of ODEs is a system of linear algebraic equations

> L_sys1 := laplace( {sys1}, t, s );

L_sys1 := {s*X(s)-x(0) = Y(s), s*Y(s)-y(0) = -4*X(s...

>

which has solution

> XY_sol1 := solve( L_sys1, {X(s),Y(s)} );

XY_sol1 := {Y(s) = (y(0)*s^3+s*y(0)*omega^2+s*omega...

>

The inverse Laplace transform provides the general solution

> q1 := invlaplace( XY_sol1, s, t );

q1 := {y(t) = 8*x(0)*sin(2*t)/((omega-2)*(omega+2))...
q1 := {y(t) = 8*x(0)*sin(2*t)/((omega-2)*(omega+2))...

>

This solution agrees with the one obtained via dsolve with method=laplace :

> dsolve( { sys1 }, { x(t), y(t) }, method=laplace );

{y(t) = 8*x(0)*sin(2*t)/((omega-2)*(omega+2))-2*x(0...
{y(t) = 8*x(0)*sin(2*t)/((omega-2)*(omega+2))-2*x(0...

>

Regardless of the method used to obtain this solution, it can be cleaned up a little to facilitate additional analysis.

> xy_sol1 := map( collect, simplify( q1, symbolic ), {x(0),y(0)} );

xy_sol1 := {x(t) = 1/2*(-4*sin(2*t)+omega^2*sin(2*t...
xy_sol1 := {x(t) = 1/2*(-4*sin(2*t)+omega^2*sin(2*t...

>

For example, the general solution to the unforced ( omega = 0 ) problem is

> eval( xy_sol1, omega=0 );

{x(t) = 1/2*y(0)*sin(2*t)+x(0)*cos(2*t), y(t) = y(0...

>

For the general forcing function with initial conditions

> ic1 := x(0)=0, y(0)=0;

ic1 := x(0) = 0, y(0) = 0

>

the solution to the IVP is

> eval( xy_sol1, {ic1} );

{y(t) = (omega*cos(2*t)-omega*cos(omega*t))/((omega...

>

The user is responsible for identifying restrictions on the applicability of these results. In this case, note that the general solution is not defined when abs(omega) = 2 . In fact, when omega = 2 , the solution to the IVP with homogeneous initial conditions is

> dsolve( eval( {sys1,ic1}, omega=2 ), { x(t), y(t) }, method=laplace );

{x(t) = -1/4*t*cos(2*t)+1/8*sin(2*t), y(t) = 1/2*t*...

>

Note the linearly growing amplitude in each component of this solution.

>

[Back to ODE Powertool Table of Contents]

>