# Extension to Higher-Order Equations Matlab Help

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).

A mass and spring with viscous surface friction

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.

Displacement and velocity of the mass as a function of time.

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

An ~ature-cOn~\Ied d~ motor.

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(

Posted on July 29, 2015 in Numerical Calculus And Differential Equations