**Extension to Higher-Order Equations**

To use the ODE solvers to solve an equation higher than order 2, you must first write the equation as a set of first-order equations. This is easily done. Consider the second-order equation

5y + 7y +4y = f(/) (8.6-1)

Solve it for the highest derivative:

J f() 4 7 . Y = S t sY – sY (8.6-2

Define two new variables XI and X2 to be y and its derivative y. That is, define These definitions imply that X2 = ~f(t) – ~XI – ~X2

This form is sometimes called the Cauchy form or the state-variable form. CAUCHY FORM Now write a function file that computes the values of XI and X2 and stores them in a column vector. To do so, we must first have a function specified for STATE·VARIABLE

f(t). Suppose that f(t) = sin t. Then the required file is FORM ——-

function xdot = example1(t,x) % Computes derivatives of two equations

xdot(l) = x(2);

xdot(2) = (1/5)*(sin(t)-4*x(l)-7*x(2));

xdot = [xdot(l); xdot(2));

Note that xdot (1) represents XI, xdot (2) represents X2, x (1) represents XI, and x (2) represents X2. Once you become familiar with the notation for the state-variable form, you will see that the preceding code could be replaced with the following shorter form:

function xdot = example1(t,x) % Computes derivatives of two equations

xdot = [x(2); (liS) * (sin(t) -4*x(1) -7*x(2)));

Suppose we want to solve (8.6-1) for 0 :5 t :5 6 with the initial conditions y(O) = 3, y(O) = 9. Then the initial condition for the vector x is [3, 91. To use ode4 5,you type [t, x) = ode45(‘example1’, [0,6), [3,9));

Each row in the matrix x corresponds to a time returned in the column vector t. If you type p lot (t ,x) ,you will obtain a plot of both X I and X2 versus t. Note that x is a matrix with two columns; the first column contains the values of XI at the various times generated by the solver. The second column contains the values

of X2. Thus to plot only XJ, type plot (t,x (:,1) ).

Solution of Nonlinear Equations We mentioned earlier that when **solving nonlinear equations,**

sometimes it is possible to check the numerical results by using an approximation that reduces the equation to a linear one. The following example illustrates such an approach with a second-order equation.

**Matrix Methods**

We can use matrix operations to reduce the number of lines to be typed in the derivative function file. For example, the following equation describes the motion of a mass connected to a spring with viscous friction acting between the mass and the surface. Another force f(t) also acts on the mass (see Figure 8.6-3).

The following function file shows how to use matrix operations. In this example m = 1, c =2, k =5, and the applied force f is a constant equal to 10. function xdot = msd(t,x) % function file for mass with spring and damping. % position is first variable, velocity is second variable. global c f k m

A = [0, l;-k/m, -elm];

B=[O;l/m];

xdot = A*x + B*f;

Note that the output xdot will be a column vector because of the definition of matrix-vector multiplication. The characteristic roots are the roots of ms2 + cs + k = s2 + 2s + 5 = 0 and are s = -1 ± 2i. The time constant is 1, and the steady-state response will thus be reached after t = 4. The period of oscillation will be 7C • Thus if we choose a final time of 5, we will see the entire

response. Using the initial conditions Xl (0) = 0, X2(0) = 0, the solver is called as follows:

global c f·k m

m = l;c = 2;k = 5;

f = 10;

[t, x] = ode23 (‘msd’, [0, 5], [0, 0]);

plot(t,x),xlabel(‘Time (s)’), …

ylabel(‘Displaeement (m) and Velocity (m/s)’), …

gtext(‘Displacernent’), gtext(‘Velocity’)

The result is shown in Figure 8.6-4.

**Test Your Understanding**

18.6-1 Plot the position and velocity of a mass with a spring and damping having the parameter values m = 2, c = 3, and k = 7. The applied force is f = 35, the initial position is y(O) = 2, and the initial velocity is

y(O) =-3.

**Characteristic Roots from the eig Function**

We saw in Section 8.4 that you can obtain the characteristic polynomial and roots for a linear ODE by substituting yet) =Aesl for the dependent variable yet). However, when the ODE is in state-variable form, there is more than one dependent variable. In such cases you must substitute Xl(t) = AlesI, X2(/) = A2esl

for the dependent variables Xl (I), X2(t), …. For example, making these substitutions into the following

Xl = -3XI +X2

X2 = -XI -7X2

gives

sAlesl = -3Alesl +A2esl

sA2esI = -AleSI -7A2esl

Collect like terms and cancel the est tenns to obtain

(s+3)AI~A2=0

AI + (s + 7)A2 = 0

As we saw in Chapter 6, a nonzero solution will exist for A I and A2 if and only if the determinant is zero. This requirement leads to Ist 3 s ~\I =S2 + 10s + 22 =0 This is the characteristic equation, and its roots are s = -6.7321 ands = -3.2679. MATLAB provides the e ig function to compute the characteristic roots when the model is given in the state-variable form (8.6-7). Its syntax is eig (A), where

A is the matrix that appears in (8.6-7). (The function’s name is an abbreviation of eigenvalue, which is another name for characteristic root.) For example, the EIGENVALUE· matrix A for the equations (8.6-8) and (8.6-9) is [-3 1] -1 -7

To find the characteristic roots, type

A = [-3, 1;-1, -7];

r = eig(A)

The answer so obtained is r = [- 6 . 7321 , – 3 . 2679] , which agrees with the answer we found by hand. To find the time constants, which are the negative reciprocals of the real parts of the roots, you type tau = -l./real (r). The time constants are 0.1485 and 0.3060.

Programming Detailed Forcing Functions As a final example of higher-order equations, we now show how’ to program

detailed forcing functions within the derivative function file. We now use a de motor as the application. The equations for an armature-controlled dc motor (such as a permanent magnet motor) shown in Figure 8.6-5 are the following. They result from Kirchhoff’s voltage law and Newton’s law applied to a rotating

inertia. The motor’s current is i, and its rotational velocity is to.

di . L- = -Rl – Kew + v(t) (8.6-10

dt dw /- = Kri – cw (8.6-11)

dtl R, and I are the motor’s inductance, resistance, and inertia; K rand K t are the torque constant and back emf constant; c is a viscous damping constant; and v(t) is the applied voltage. These equations can be put into matrix form as follows.

where Xl = i and X2 = W. [;:] = [~; ~11′[;:] +mv(