6.5. Linear System Applications

Systems of linear equations are routinely solved using matrices in science and engineering. Here, examples of two areas of application are reviewed.

6.5.1. DC Electric Circuit

A direct current (DC) electrical circuit with only a battery and resistors provides a system of linear equations. We want to determine the current and voltage drop across each resistor in the circuit. We only need Ohm’s Law and either Kirchhoff’s Voltage Law or Kirchhoff’s Current Law to find a solution.

Note

An AC circuit with capacitors and inductors in addition to resistors is described by a system of differential equations that can also be solved with linear algebra. MATLAB and linear algebra make determining such solutions straight forward; however, additional mathematical abstractions are needed, so we’ll stick with the simpler DC circuit with only resistors.

Ohm's Law:

When a current of I amperes flows through a resistor of R ohms, the voltage drop, V, in volts across the resistor is

V = I\,R.

Kirchhoff's Voltage Law (KVL)

Each closed loop of a circuit is assigned a loop current. The current through resistors that are part of two loops is a combination of the loop currents. Then, the sum of the voltage drops around each loop is zero.

Kirchhoff's Current Law (KCL)

The sum of the currents entering and leaving each node of the circuit is zero. The voltages at each node of the circuit may be quatified in terms of any known voltages, voltage drops accross resistors, and other unknown node voltages.

Both the KVL and KCL methods lead directly to matrix equations. You will notice, however, some significant differences. The KVL method leads to a fairly small matrix that is quick to set-up; but after the matrix is solved, a little additional work is needed to relate the loop currents to voltage drops across each resistors. On the other hand, the KCL method leads to a larger, but more sparse matrix. Defining the KCL matrix equation is a little more involved, but the results are quickly related to voltage drops across the resistors.

In practice, when circuits are solved by hand, with the help of programs like MATLAB. the KVL method seems a likely plan of attack for most people. But it is easier for automatic circuit solving programs, such as PSpice, to use the KCL method.

6.5.1.1. KVL Method

The circuit is shown below with the loop currents and their directions identified. The constant values used in the MATLAB program are: R1 = 1000; R2 = 1500; R3 = 2000; R4 = 1000; R5 = 1500; v = 12.

../_images/KVLcircuit.png

The three loop equations are given in terms of the loop currents.

\spalignsys{ i_1R_1 + (i_1 - i_2)R_3 = v;
i_2R_2 + (i_2 - i_3)R_4 + (i_2 - i_1)R_3 = 0;
i_3R_5 + (i_3 - i_2)R_4 = 0}

In matrix notation, this becomes:

\spalignmat{(R_1+R_3) -R3 0;-R3 (R2+R3+R4) -R4; 0 -R4 (R4+R5)}
\spalignvector{i_1 i_2 i_3} = \spalignvector{v 0 0}

6.5.1.2. KCL Method

The KCL method has unknown variables for each current entering or exiting a node and the node voltages, thus it has more equations than the KVL method.

../_images/KCLcircuit.png

\spalignsys{ v_1 + i_1R_1 \. \. = v;
v_1 - i_3R_3 \. \. = 0;
v_1 - v_2 - i_2R_2 = 0;
v_2 - i_4R_4 \. \. = 0;
v_2 - i_5R_5 \. \. = 0;
i_1 - i_2 - i_3 = 0;
i_2 - i_4 - i_5 = 0}

\spalignmat{1 0 R1 0 0 0 0;
        1 0 0 0 -R3 0 0;
        1 -1 0 -R2 0 0 0;
        0 1 0 0 0 -R4 0;
        0 1 0 0 0 0 -R5;
        0 0 1 -1 -1 0 0;
        0 0 0 1 0 -1 -1}
\spalignvector{v_1 v_2 i_1 i_2 i_3 i_4 i_5} =
\spalignvector{v 0 0 0 0 0 0}

6.5.1.3. MATLAB Code

After doing the KVL calculation, the code computes some node voltages to help verify the correctness of the voltages from the KVL and KCL solutions.

%% Constants
R1 = 1000;
R2 = 1500;
R3 = 2000;
R4 = 1000;
R5 = 1500;
v = 12;

%% KVL Solution

R = [(R1+R3) -R3 0;
    -R3 (R2+R3+R4) -R4;
    0 -R4 (R4+R5)];
V = [v; 0; 0];
I = R \ V;      % V = RI
i1 = I(1);
i2 = I(2);
i3 = I(3);
V1 = R1*i1;
V2 = R2*i2;
V3 = R3*(i1-i2);
V4 = R4*(i2-i3);
V5 = R5*i3;
fprintf('\n\tKVL Solution:\n')
fprintf('Current:\n i1 = %f\n i2 = %f\n i3 = %f\n',i1,i2,i3);
fprintf('Voltages:\n V1 = %f\n V2 = %f\n V3 = %f\n V4 = %f\n V5 = %f\n', ...
            V1,V2,V3,V4,V5);
disp('Verify: ')
fprintf(' %f == %f ?\n',V4,V5);
fprintf(' %f == %f ?\n',V3,(V2+V4));
fprintf(' %f == %f ?\n',v,(V3+V1));

%% KCL Solution

A = [1 0 R1 0 0 0 0;
    1 0 0 0 -R3 0 0;
    1 -1 0 -R2 0 0 0;
    0 1 0 0 0 -R4 0;
    0 1 0 0 0 0 -R5;
    0 0 1 -1 -1 0 0;
    0 0 0 1 0 -1 -1];
b = [v 0 0 0 0 0 0]';
x = A \ b;
fprintf('\n\tKCL Solution:\n')
fprintf('Current:\n i1 = %f\n i2 = %f\n i3 = %f\n',x(3),x(4),x(5));
fprintf(' i4 = %f\n i5 = %f\n',x(6),x(7));
fprintf('Voltages:\n V1 = %f\n V2 = %f\n ', R1*x(3),R2*x(4));
fprintf('V3 = %f\n V4 = %f\n V5 = %f\n', R3*x(5),R4*x(6),R5*x(7));

This program displayed the following output. One can see that the KVL and KCL methods gave the same results. Remember that the currents reported for the KVL method are loop currents, not the currents through all of the resistors. The current through R_3, and R_4 are combinations of the loop currents.

    KVL Solution:
Current:
 i1 = 0.005928
 i2 = 0.002892
 i3 = 0.001157
Voltages:
 V1 = 5.927711
 V2 = 4.337349
 V3 = 6.072289
 V4 = 1.734940
 V5 = 1.734940
Verify:
 1.734940 == 1.734940 ?
 6.072289 == 6.072289 ?
 12.000000 == 12.000000 ?

    KCL Solution:
Current:
 i1 = 0.005928
 i2 = 0.002892
 i3 = 0.003036
 i4 = 0.001735
 i5 = 0.001157
Voltages:
 V1 = 5.927711
 V2 = 4.337349
 V3 = 6.072289
 V4 = 1.734940
 V5 = 1.734940

6.5.2. The Statics of Trusses

The field of statics is concerned with determining the forces acting on civil and mechanical systems. Knowing the magnitude of forces on parts of a system are needed to ensure that the systems are designed to be strong enough.

Here we consider a structure made of seven trusses. This example was taken from a text on statics and dynamics. [MERIAM78] The page from Meriam’s book is here. truss_problem.pdf

Statics problems such as this one are solved at two levels. First, we consider the external forces on the overall structure. This will tell us the forces on the supports of the structure. Three equations are used for this step.

  1. Sum of the moments is zero. Here we define a focal point and then calculate the moment of each external force as the product of the force value perpendicular to the structure times distance from the focal point. This ensures that the structure does not rotate.

    \sum_{\forall i} F_i\,d_i = 0

  2. Sum for forces in x direction is zero.

    \sum_{\forall i} F_{i,x} = 0

  3. Sum for forces in y direction is zero.

    \sum_{\forall i} F_{i,y} = 0

Next, use the method of joints to determine the force on each truss. We examine each joint and write equations for the forces in the x and the y directions. These forces must also sum to zero for each joint.

As you can see from the scanned page from Meriam’s text book, this is a fairly simple problem. When taken in the desired order, each equation reduces to having only one unknown variable.

When solved with a matrix, each truss force can be entered as a compression force. Trusses that are in tension will be found to have a negative force.

../_images/truss_prob.png

The matrix equation is overdetermined, meaning that we have more equations than unknowns. When this occurs, we multiply both sides of the equation by A^T.

A\,x = b \mbox{  becomes  } A^TA\,x = A^Tb

c60 = cosd(60);
c30 = cosd(30);
s30 = sind(30);
s60 = sind(60);
% Method of joints
%           AB AC   BC BD  CD  CE  DE    T Ey Ex
joints = [-s60  0    0  0   0   0   0    0 0  0; % Ay
          -c60 -1    0  0   0   0   0    0 0  0; % Ax
           s30  0 -c60 -1   0   0   0    0 0  0; % Bx
           c30  0  c30  0   0   0   0    0 0  0; % By
             0  0 -s60  0 -s60  0   0    0 0  0; % Cy
             0  1  c60  0 -c60 -1   0    0 0  0; % Cx
             0  0    0  0  c30  0  c30 s30 0  0; % Dy
             0  0    0  1  c60  0 -c60 c30 0  0; % Dx
             0  0    0  0   0   0 -c30   0 1  0; % Ey
             0  0    0  0   0   1  c60   0 0 -1];% Ex
M =  [0 0 0 0 0 0 0   5 0  0];  % sum of moments
fy = [0 0 0 0 0 0 0 s30 1  0];  % sum of external y forces
fx = [0 0 0 0 0 0 0 c30 0 -1];  % sum of external x forces
b = [30 0 0 0 20 0 0 0 0 0 400 50 0]';
A = [joints; M; fy; fx];
A1 = A'*A;
B1 = A'*b;
name = {'AB','AC','BC','BD','CD','CE','DE','T','Ey','Ex'};
x = A1\B1;
x(8) = -x(8);       % T is in tension
for ii = 1:length(name)
    showForces(x(ii),name{ii});
end

function showForces(x, name)
    if x < 0
        f = 'tension';
        x = -x;
    else
        f = 'compression';
    end
    fprintf('%s has %f kN of %s force\n', name, x, f);
end