Successive over Relaxation (SOR) Method to Solve System of Equations

Problem: Develop a MATLAB code to solve the following system of algebraic equations using the Successive-over-Relaxation Method.

Solution:

Successive over Relaxation Method: This method is just the modification of the Gauss-Seidel method with an addition relaxation factor 𝛚. This method gives convergent solution as there is an option for under relaxation when 𝛚 is less than one. For different 𝛚, the following program can determine the solution. However, when 𝛚 is 0.05, the relatively better results are obtained.

MATLAB Program for Successive Over Relaxation (SOR) Method

function x = sor(a,b,x0,e,N,w)

% a - (nxn) matrix

% b - column vector of length n

% x0 - initial solution

% e - error definition

% N - no. of iterations

% w – relaxation factor

 

format long;

 

m=size(a,1); 

 

n=length(b);

 

p=length(x0);

 

q=length(e);

 

if (m ~= n)

    

    error('a and b do not have the same number of rows')

    

end

 

y0=x0;

 

for i=1:m

   

    [~,f]=max(abs(a(1:m,i:i))); 

                             

   % disp(f);  

 

    a([i,f],1:n)=a([f,i],1:n);                 

                                                                                       

end

 

for t=1:m

    

if a(t,t)==0

    

   [~,g]=max(abs(a(1:m,t:t)));            

        

   a([t,g],1:n)=a([g,t],1:n);

        

end

 

end

    

% disp(a);

 

c=1;

 

while(c<N)

 

for j=1:m  

    

    for k=1:p

        

        Z(1,k)=(a(j,1:n)*y0(1:p,1));

                                                                    

    end

    

    Y(j:j,1)=sum(Z(1,k)); 

       

   % disp(Y);

           

    R(j:j,1)=b(j:j,1)-Y(j:j,1);

    

   % disp (R);

                     

    x(j:j,1)=y0(j:j,1)+((R(j:j,1)*w)/a(j:j,j:j));

            

    y0(j:j,1)= x(j:j,1); 

    

   % disp (y0);    

   % disp (x);

               

end

 

if (abs(x(1:p,1)-x0(1:p,1)) < e(1:q,1))

    

    break;

    

end

 

disp (x);

disp(abs(x(1:p,1)-x0(1:p,1)));

    

x0=y0;

 

c=c+1;

 

end

 

x=x(1:p,1);

       

end

 

Program Outputs:

ans =

   2.322822494693882

   0.675693020694401

   2.160169880982512

  -0.350429027969758

  -0.428941223493063

  -0.394923172868219 

 

Error Percentage: 2.2%, 1.5%, 1.4%, 2.6%, 0.9% and 5.8%.

Gauss-Seidel Iteration Method to Solve System of Algebraic Equations

The following matrix represents a system of linear algebraic equations. In this problem, we are going to solve these equations by applying Gauss-Seidel iteration approach.

The following program solves system of linear algebraic equations iteratively with successive approximation by using most recent solution vectors. However, there is a condition to work for this program which is strictly diagonal dominance. And, the given matrix in question is not diagonally dominant although matrix switching operations are carried out. For this circumstance, convergence is not guaranteed and the program shows the problem while calculating error as it increases with no of iterations.

MATLAB Program for Gauss-Seidel Iterative Method

% Gauss-Seidel Iterative Algorithm

 

function x = gauss_seidal(a,b,x0,e,N) 

 

% a - (nxn) matrix

% b - column vector of length n

% x0 - initial solution

% e - error definition

% N - no. of iterations

 

format long;

 

m=size(a,1);   % get number of rows in matrix a

 

n=length(b);   % get length of b

 

p=length(x0);  % get length of x0

 

q=length(e);   % get length of e

 

if (m ~= n)

    

    error('a and b do not have the same number of rows')

    

end

 

y0=x0;

 

for i=1:m       % Interchanging row/column to make diagonally dominant matrix

   

    [~,f]=max(abs(a(1:m,i:i))); 

                             

   % disp(f);  

 

    a([i,f],1:n)=a([f,i],1:n);                 

                                                                                       

end

 

for

             % while switching rows. This loop just interchange between rows to avoid '0'

    

if a(t,t)==0

    

   [~,g]=max(abs(a(1:m,t:t)));            

        

   a([t,g],1:n)=a([g,t],1:n);

        

end

 

end

    

 disp(a);

 

c=1;

 

while(c<N)

 

for j=1:m   % Successive iteration, most recent solution used for next calculation

    

    for k=1:p

        

        Z(1,k)=(a(j,1:n)*y0(1:p,1));

                                                                    

    end

    

    Y(j:j,1)=sum(Z(1,k)); 

       

   % disp(Y);

           

    R(j:j,1)=b(j:j,1)-Y(j:j,1);

    

   % disp (R);

                     

    x(j:j,1)=y0(j:j,1)+(R(j:j,1)/a(j:j,j:j));

            

    y0(j:j,1)= x(j:j,1); 

    

   % disp (y0);    

   % disp (x);

end

 

if (abs(x(1:p,1)-x0(1:p,1)) < e(1:q,1)) % Error checking

    

    break;

    

end

 

disp (x);

disp(abs(x(1:p,1)-x0(1:p,1)));

    

x0=y0;

 

c=c+1;

 

end

 

x=x(1:p,1); % Returns solution

       

end

 

Program Outputs:

>> gauss_seidal([1 -1 4 0 2 9; 0 5 -2 7 8 4; 1 0 5 7 3 -2; 6 -1 2 3 0 8; -4 2 0 5 -5 3; 0 7 -1 5 4 -2], [19;2;13;-7;-9;2], [0; 0;0; 0; 0; 0],[0.05;0.05;0.05;0.05;0.05;0.05], 3)  % function calling

     6    -1     2     3     0     8    % rearranged matrix

     0     7    -1     5     4    -2

     1     0     5     7     3    -2

    -4     2     0     5    -5     3

     0     5    -2     7     8     4

     1    -1     4     0     2     9

   3.166666666666667   % first improved solutions

   0.285714285714286

   1.966666666666667

   1.019047619047619

  -1.703571428571429

  -0.593386243386244

   3.166666666666667   % error

   0.285714285714286

   1.966666666666667

   1.019047619047619

   1.703571428571429

   0.593386243386244

   2.840388007054674   % second improved solutions

   0.642705971277400

   1.390044091710759

  -0.732311665406903

  -0.241714380196523

  -0.586047738025251

   0.326278659611992   % error

   0.356991685563114

   0.576622574955908

   1.751359284454522

   1.461857048374905

   0.007338505360992

ans =                 % final solutions

   2.840388007054674

   0.642705971277400

   1.390044091710759

  -0.732311665406903

  -0.241714380196523

  -0.586047738025251

 

Percentage of error for the six solutions: 10.12%,  125%,  29%,  173%,  85%  and  1.1%. 

These errors may be the results of having a matrix, which is not diagonally dominant and the solution diverges as iteration increases.

Application of Gaussian Elimination to Solve System of Linear Algebraic Equations

Problem: In this example, we are going to solve the following system of linear algebraic equations by the widely-used Gaussian elimination method. Here, we are going to develop a MATLAB program that implements the Gaussian elimination to solve the following linear algebraic equations.

The following program implements Gaussian elimination method with partial pivoting and scaling to solve system of linear algebraic equations. The explanations of the codes are mentioned just right of each line.

MATLAB Program for Gauss Elimination Method

% Gaussian elimination with partial pivoting and scaling

 

function x = gausselimination(a,b)

 

% a - (nxn) matrix

% b - column vector of length n

 

format long;

 

m=size(a,1); % get number of rows in matrix a

 

n=length(b); % get length of b

 

if (m ~= n)

    

    error('a and b do not have the same number of rows')

    

end

 

a(:,n+1)=b;  % Forming (n,n+1) augmented matrix

 

for c=1:n

              

for r=c:m

     

% Scaling of the input matrix "a"

       

    Z(r,c) = a(r,c)/max(a(r,c:n));  % Scaling prior to pivoting 

 

    K(r:r,c:c)=Z(r,c); % Normalized coefficient parameters saved temporality 

                       % to another matrix "K"     

end

 

disp(K);

           

[~,f]=max(abs(K(1:r,c:c))); % Finding maximum absolute value from "K" matrix 

                            % for switching rows in matrix "a"      

disp(f);  

 

a([c,f],1:n+1)=a([f,c],1:n+1); % Switching between two rows 

                                                                                    

disp(a);

                                        

for i=c

% Making diagonal element of matrix "a" into 1.0

    

    a(i,i:n+1) = a(i,i:n+1)/a(i,i);   

 

for j=i+1:n

    

    % Making all elements below the diagonal element into zero

 

    a(j,i:n+1) = a(j,i:n+1) - a(j,i)*a(i,i:n+1);

    

end

 

end

 

disp(a);

 

end

 

% Process of back substitution

 

for j=n-1:-1:1

 

    a(j,n+1) = a(j,n+1) - a(j,j+1:n)*a(j+1:n,n+1);

 

end

    

% Returning solution

 

x=a(:,n+1);

 

end

Program Outputs:

ans =

 

  -1.761817043997860

   0.896228033874012

   4.051931404116157

  -1.617130802539541

   2.041913538501913

   0.151832487155935

Therefore, the solutions of the six linear algebraic equations are -1.76; 0.89; 4.05; -1.62; 2.04 and 0.15.

A PID Controller Design Method for DC Motor Speed Control

Control systems analysis and design focuses on three primary objectives:

•        producing the desired transient response
•        reducing steady state errors
•        achieving stability

Controller: Additional component or device that equalizes or compensates for the performance deficiency is called compensator or controller.

Plant: A system to be controlled.

Controller: Provides the excitation for the plant; Designed to control the overall system behavior.

There are several techniques available to the control systems engineer to design a suitable controller.

One of controller widely used is the proportional plus integral plus derivative (PID) controller, which has a transfer function:

Kp = Proportional gain

KI = Integral gain

Kd = Derivative gain

The signal (u) just past the controller is now equal to the proportional gain (Kp) times the magnitude of the error plus the integral gain (Ki) times the integral of the error plus the derivative gain (Kd) times the derivative of the error.

MATLAB Program:

Design requirements

Settling time less than 1 seconds

Overshoot less than 5%

Steady-state error less than 1%

Open-loop transfer function of the DC Motor:

Proportional Control:

j=3.2284e-6;

b=3.5077e-6;

k=0.0274;

r=4;

l=2.75e-6;

num=k;

kp =1.7;

den=[(j*l) ((j*r)+(l*b)) ((b*r)+k^2)  0];

sys=tf(num,den);

feedbk=feedback(kp*sys,1);

step(feedbk,0:.001:1)

Response Curve:

Proportional+Integral Control:

J=0.01;

b=0.1;

K=0.01;

R=1;

L=0.5;

num=K;

den=[(J*L) ((J*R)+(L*b)) ((b*R)+K^2)];

open_sys= tf(num,den);

subplot(2,1,1),step (open_sys,0:0.1:3);

step (open_sys,0:0.01:3);

title('Step Response for the Open Loop System');

kp = 100;

ki = 200;

kd = 0 ;

controller=tf([kd kp ki],[1 0]);

sys=feedback(controller*open_sys,1);

subplot(2,1,2),step (sys,0:0.01:3);

title('Closed loop Step:ki=100 kp=200');

Response Curve:

#Blog #Blogger #PID #PIDcontroller #DCmotor #ResponseCurve

A Short MATLAB Script for Spectrum Analysis or FFT

The following codes may be used for spectrum analysis for any vibration problems. The program is capable of isolating fundamental peaks from a noisy signal smoothly. It uses the MATLAB built-in Fast Fourier Transform (FFT) algorithm.

Codes for Spectrum Analysis (FFT):

NS=EFFORTS(:,1);                 % Data to Analyze
L=length(NS);
Fs=L;                            % Sample Rate
NFFT = 2^nextpow2(L);
F = Fs/2*linspace(0,1,NFFT/2+1);
win=hamming(L,'periodic');
Ywin=win.*NS;
X = fft(Ywin,NFFT)/L;            % MATLAB built-in FFT Algorithm
Xf=abs(X(1:NFFT/2+1));
Xf=Xf/max(Xf);                   % Normalizing Magnitude
plot(F,Xf)

Observation of a Typical Open Loop Dynamic Responses of Second Order Systems

Dynamic Responses

In vibration and control, we always deal with the following three types of the dynamic responses from a system due to the step input to that system. These responses are: over damped, critically damped, and under damped signals.

• Overdamped - where the damping ratio is higher than the expected or optimum. As the name suggests, 'damped' means to stabilize or settle down a signal, and over-damped is in the same perspective that the signal is sized or filtered overly to reach a desired value. In the following response curves, the yellow colored signal represents an overdamped response.
• Underdamped - where the damping ratio is below the desired limit. Again, from its title, 'under' defines that the signal is not sufficiently filtered or processed to reach the optimum or desired level. The light blue colored signal in the following response curve shows the underdamped signal response. We see that the signal has couple of overshoots and undershoots until it settles down to the desired value.
• Critically Damped - This is the response that we desire to achieve, which means that the system response to the step input should have a shape that is shown by the pink colored signal in the following response curve.

SIMULINK Representation

The following block diagram represents the generation of three types of responses using three SIMULINK transfer function blocks.

Response Curves

#OverDamped   #UnderDamped   #CriticallyDamped   #DynamicResponse

Modelling a Basic Second Order System in SIMULINK

Problem: Consider the mechanical system shown in the figure below. The system is at rest for t < 0. At t = 0 an impulse‐force of 8 N is applied to the mass. The impulse lasts 0.06 seconds. The displacement x is measured from the equilibrium position just before the mass is hit by the impulse force. Find the mathematical model of the system, draw a block diagram to represent it, and plot the displacement and velocity of the mass during the first 25 seconds.

Figure 1: Diagram of the System.

Solution: The mathematical model is given by the following equation:

The block diagram for the system is given below:

Figure 2: Block Diagram.

Figure 3: Simulink Block Diagram.

The output variables (i.e., displacement and velocity of the mass) displayed by the scope is shown in the figure below.

After examining the SIMULINK model, you can see that blocks from the following libraries have been used:

• Sources library: the Pulse Generator block (the Clock block was not used because it was not requested to send any variable to the MATLAB Workspace)
• Sinks library: The Scope block
• Math Operations library: The Sum and Gain blocks
• Continuous library: The Integrator block
• Signal Routing library: The Mux block

Figure 4: Response Curve.

You should have noticed that two of the gain blocks (damper and spring) have been rotated 180 degrees with respect to its default orientation. This has been done by right clicking on the corresponding block and then selecting flip block from the format menu. In this case the applied load was modeled using the Pulse Generator block. The dialog box containing the parameters for this block is shown below.

An Iterative Approach to Solve a System of Equations

In this tutorial, we are going to solve the following two algebraic equations iteratively.

X₁ = X₂ - 1

X₂ = 2.5 X₁ - 0.5

Using fixed point iterative method and taking initial assumptions according to the question as, X₁ = X₂ = 0, the solutions do not converge with this present arrangement of equations.

However, if we rearrange these equations in the following forms, the solutions indeed converge.

X2 = X1 + 1

X₁ = X₂ / 2.5 + 0.2

A MATLAB program has been developed to solve the modified equations iteratively by creating a user defined function called "iterative". This function may be used to solve other equations iteratively which is the essence of the user defined function. The previous forms of equations from the question are also fed into the program to see the solutions, but the program indicates divergence.

The following MATLAB program written in “m-file” subsequently enlightens essential commands in the program by additional notes.

MATLAB Codes:

% Definition of a function "iterative" to solve equations iteratively

function [x1, x2, ea1, ea2] = iterative(X1, X2, etol1, etol2)

% Input X1=0, X2=0, etol1, etol2=Tolerance definition for errors

% Output x1, x2=Solutions of equations, ea1, ea2= Calculated errors in loop

% Program Initialization

solution1 = 0;

solution2 = 0;

% Iterative Calculations

while (1)

    solutionprevious1 = solution1;

   

    solutionprevious2 = solution2;

   

    solution2 = X1+1;             % Problem equation 1

   

    solution1 = (X2/2.5)+0.2;     % Problem equation 2

    

    X1=solution1;

  

    X2=solution2;

if solution1~=0 && solution2~=0

% Approximate percent relative errors

    ea1=abs((solution1 - solutionprevious1)/solution1)*100;

   

    ea2=abs((solution2 - solutionprevious2)/solution2)*100;

end

if ea1<=etol1 && ea2<=etol2,  % Conditions to meet specified error tolerances

  

    break,

end

x1 = solution1;

x2 = solution2;

end

% Display of output parameters

disp(x1);                     

disp(x2);

disp(ea1);

disp(ea2);

end

Outputs:

>> iterative (0,0,1e-5,1e-5)      % Command input at the MATLAB command prompt

   1.0000

   2.0000

   3.4360e-06

   8.5899e-06

Therefore, the solutions obtained from the program are X₁ = 1 and X₂ = 2. And, the approximate percent relative errors for finding X₁, X₂ are 3.4360e-06 and 8.5899e-06 respectively which also satisfy the condition mentioned in the program.