ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Unit 35 -- Application: Transfer Functions and Frequency Response Curves

Outline of Unit 35

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(t)=laplace(x(t),t,s) ):

>

The transfer function for a linear ODE is the Laplace transform of the fundamental solution to the ODE. While this definition applies to an ODE any order, only second-order linear ODEs will be discussed here. In particular, consider

> ode1 := m*diff( x(t), t$2 ) + b*diff( x(t), t ) + k*x(t) = g(t);

>

The fundamental solution can be obtained as the solution for the IVP with and homogeneous initial conditions.

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

>

The definition suggests finding the transfer function from the fundamental solution of the ODE :

> q1 := dsolve( { eval(ode1,g(t)=Dirac(t)), ic1 }, x(t), method=laplace );

> q2 := laplace( rhs(q1), t, s );

> H1 := 1/expand(simplify(1/q2));

>

However, it is generally more efficient to simply take the Laplace transform of the IVP and solve for the Laplace transform of the fundamental solution. That is,

> H2 := solve( eval( laplace( eval(ode1,g(t)=Dirac(t)), t, s ), {ic1} ), X(t) );

>

These results are obviously equivalent.

>

Note that since = = 1, an equilvalent definition of the transfer function is the ratio of the Laplace transform of the solution (with homogeneous initial conditions) and the Laplace transform of the forcing function. That is, = = .

>

The frequency response curves, or Bode plots, are plots of the magnitude and argument (in degrees) of the transfer function evaluated at ( ) as functions of the frequency . That is, if the transfer function is H, then and . The functions M and N correspond to the amplitude and phase shift, respectively, of the steady-state solution to the linear ODE with sinusoidal forcing function, , i.e., .

For example, when

> H := unapply( H1, s );

>

the amplitude and argument (in degrees) of the steady-state solution with frequency are

> M := unapply( evalc( abs( H(I*omega) ) ), (m,b,k) );

> N := unapply( evalc( argument( H(I*omega) ) ) * 180/Pi, (m,b,k) );

>

The quantity is the gain, in decibels (dB). Bode plots typically display the gain (in decibels) and phase angle (in degrees) versus the angular frequency on a logarithmic scale for .

>

For example, when , , and ,

> semilogplot( 20*log[10](M(1,1/4,2)), omega=0.1..10,

> axes=framed, labels=[omega,"gain (dB)"],

> labeldirections=[HORIZONTAL,VERTICAL] );

>

> semilogplot( N(1,1/4,2), omega=0.1..10, axes=framed,

> labels=[omega,"phase shift (degrees)"],

> labeldirections=[HORIZONTAL,VERTICAL] );

>

One important use of the Bode plot for the magnitude is to determine the resonant frequency. The resonant frequency for a system occurs when the magnitude of the steady-state response is maximized.

> dM := diff( M(m,b,k), omega );

> q3 := solve( dM=0, omega );

> res_freq := omega = op(select( w->evalb(eval(w,{m=1.,b=0,k=1.})>0), [q3] ));

>

Note that for the undamped system, the resonant frequency is found to be

> simplify( eval( res_freq, b=0 ) );

>

as expected.

>

Another use of Bode plots is to investigate the behavior of a system as one of the parameters changes. To illustrate, the gain and phase angle are investigated as the mass and spring constant are held fixed and the damping coefficient decreases to zero.

Consider the damped spring-mass system with = 1, = 4, and damping coefficients selected from

> damp_vals := [2^(-i) $ i=-3..8, 0];

>

Also, chose an interval of positive frequencies to use as the domain for the Bode plots.

> freq_range := 0.1..10;

>

The animated Bode plots of the magnitude for each damping coeffcient is

> display( [seq( semilogplot( 20*log[10](M(1,b,4)), omega=freq_range,

> axes=framed, labels=[omega,"gain (dB)"],

> labeldirections=[HORIZONTAL,VERTICAL] ),

> b=damp_vals )], insequence=true );

>

Note that there is no (interior) maximum amplitude in the first two frames of this animation ( = 8 and = 4). This is consistent with the fact that a damped spring-mass system has a resonant frequency only when < .

The animated Bode plots for the phase angle for each damping coefficient is

> display( [seq( semilogplot( N(1,b,4), omega=freq_range,

> axes=framed, labels=[omega,"phase shift (degrees)"],

> labeldirections=[HORIZONTAL,VERTICAL] ),

> b=damp_vals )], insequence=true );

>

Note the jump in the phase shift when . The argument jumps to 180 degrees instead of -180 degrees because Maple's arctangent function returns values in ( , ].

And, now, the two Bode plots with all 13 frames superimposed on the same set of axes.

> semilogplot( [seq( 20*log[10](M(1,b,4)), b=damp_vals )],

> omega=freq_range, axes=framed,

> labels=[omega,"gain (dB)"],

> labeldirections=[HORIZONTAL,VERTICAL] );

>

> semilogplot( [seq( N(1,b,4), b=damp_vals )],

> omega=freq_range, axes=framed,

> labels=[omega,"phase shift (degrees)"],

> labeldirections=[HORIZONTAL,VERTICAL] );

>

>