Welcome to the World of Modelling and Simulation

What is Modelling?

This blog is all about system dynamics modelling and simulation applied in the engineering field, especially mechanical, electrical, and ...

Finite Difference Approach by MATLAB for the First and Second Derivatives

The following MATLAB program determines the first and second derivatives of the data given in the problem applying the finite difference schemes and developing a custom user defined function firstsecondderivatives(x,y).

The finite difference schemes used are the following:

For first derivative:
  • At first point (three point forward)
  • At last point (three point backward)
  • At all other points (two point central)

For second derivative:
  • At first point (four point forward)
  • At last point (four point backward)
  • At all other points (three point central)

DATA SET
 x -1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
 f(x) -3.632 -0.3935 1 0.6487 -1.282 -4.518 -8.611 -12.82 -15.91 -15.88 -9.402 9.017


function [dy, ddy] = firstsecondderivatives(x,y)

% The function calculates the first & second derivative of a function that is given by a set
% of points. The first derivatives at the first and last points are calculated by
% the 3 point forward and 3 point backward finite difference scheme respectively.
% The first derivatives at all the other points are calculated by the 2 point
% central approach.

% The second derivatives at the first and last points are calculated by
% the 4 point forward and 4 point backward finite difference scheme respectively.
% The second derivatives at all the other points are calculated by the 3 point
% central approach.

n = length (x);
dy = zeros;
ddy = zeros;

% Input variables:
% x: vector with the x the data points.
% y: vector with the f(x) data points.

% Output variable:
% dy: Vector with first derivative at each point.
% ddy: Vector with second derivative at each point.

dy(1) = (-3*y(1) + 4*y(2) - y(3)) / 2*(x(2) - x(1)); % First derivative
ddy(1) = (2*y(1) - 5*y(2) + 4*y(3) - y(4)) / (x(2) - x(1))^2; % Second derivative

for i = 2:n-1
   
dy(i) = (y(i+1) - y(i-1)) / 2*(x(i+1) - x(i-1));
ddy(i) = (y(i-1) - 2*y(i) + y(i+1)) / (x(i-1) - x(i))^2;

end

dy(n) = (y(n-2) - 4*y(n-1) + 3*y(n)) / 2*(x(n) - x(n-1));
ddy(n) = (-y(n-3) + 4*y(n-2) - 5*y(n-1) + 2*y(n)) / (x(n) - x(n-1))^2;

figure(1)
subplot (3,1,1)
plot (x, y, 'linewidth', 2, 'color', 'k'), grid
title('x VS f(x)')
ylabel ('f(x)')

subplot (3,1,2)
plot (x, dy, 'linewidth', 2, 'color', 'b'), grid
title('x VS dy')
ylabel ('dy: First Derivative')

subplot (3,1,3)
plot (x, ddy, 'linewidth', 2, 'color', 'r'), grid
title('x VS dy')
ylabel ('ddy: Second Derivative'), xlabel('x')

figure(2)
plot (x, y, 'linewidth', 2, 'color', 'k');
hold on;
plot (x, dy, 'linewidth', 2, 'color', 'b');
hold on;
plot (x, ddy, 'linewidth', 2, 'color', 'r'), legend
ylabel ('y, dy, ddy'), xlabel('x')
end


The following operations are done in MATLAB command window to run the above function.

Program Outputs:

>> x = -1:0.5:4.5;
>> y=[-3.632 -0.3935 1 0.6487 -1.282 -4.518 -8.611 -12.82 -15.91 -15.88 -9.402 9.017];
>> firstsecondderivatives(x,y)

ans =

  Columns 1 through 11

    2.0805    2.3160    0.5211   -1.1410   -2.5833   -3.6645   -4.1510   -3.6495   -1.5300    3.2540   12.4485

  Column 12

   12.1947

We can also plot the derivatives, below is the figure for the function, first and second order derivatives respectively.

Plotting the first and second order numerical derivatives
Plotting the first and second order numerical derivatives

Learning Mathematica, Lesson 2: Solving Euler-Bernoulli Beam Equation

In this lesson, I would like to show the advantages of the Mathematica built-in solver to evaluate the analytical solution of a differential equation. For example, if we want to solve the well-known fourth order Euler-Bernoulli equation to solve a problem of a cantilever beam, the Mathematica code has very user-friendly features to do so. In following figure, we see that the beam has one fixed end and the other end has a concentrated load acting downward.

The following Mathematica program determines the analytical displacement function and plots the result. Most of the lines are commented as explanations for the line of codes to understand the respective operation.

Cantilever beam with a fixed end and concentrated load at free end
Cantilever beam with fixed end


(*Determination of the Exact Solution of an Euler-Bernoulli Beam*)
Clear[y, P, x, EI, L]
y[x]; (*The Final Solution, y*)
(*Defining the constant parameters, just assuming unit values*)
h = 1; (*Beam height*)
P = 1; (*Concentrated load at the end*)
L = 10; (*Beam length*)
b = 1; (*Beam width*)
Ie = 1/ 12; (*Moment of inertia*)
Ee = 10 000; (*Young's Modulus*)
EI = Ee* Ie; (*Flexural Modulus*)
th = D[y[x], x]; (*Slope*)
V = EI* D[y[x], {x, 3}]; (*Shear force*)
M = EI* D[y[x], {x, 2}]; (*Bending moment*)
(*Defining the boundaries of the problem*)
y1 = y[x] /. x → 0;
y2 = y[x] /. x → L;
th1 = th /. x → 0;
th2 = th /. x → L;
M1 = M /. x → 0;
M2 = M /. x → L;
V1 = V /. x → 0;
V2 = V /. x → L;
(*Using DSolve a built-in solver to solve beam equation*)
s = DSolve[{EI* y''''[x] ⩵ 0, y1 ⩵ 0, th1 ⩵ 0, V2 ⩵ P, M2 ⩵ 0}, y, x];
Print["The Exact Displacement Function"]
y = Simplify[y[x] /. s[[1]]]
Plot[y, {x, 0, 5}, PlotStyle → {Green, Thick}]
Print["Displacement at the tip of the beam in Y direction"]


Mathematica generates the analytical solution function which is seen below. This feature is very useful as you may solve any differential equations analytically and get the results symbolically as well like this problem.


analytical solution function

Linear Least Squares Regression Analysis by a MATLAB program

A MATLAB program is developed to determine the coefficients by linear least squares regression where the function is, y = mx + b. Here,


x
-6
-4
-1
0
3
5
8
y
18
13
6
4
-1
-8
-15


% A function defined for least square regression
function [c,R2] = linearregression(x,y)

% Least-squares fit of data to y = c(1)*x + c(2)
% Here, c(1) = m; c(2) = b;
% Inputs:
% x,y = Vectors of independent and dependent variables
% Outputs:
% c = Coefficients of the given equation

if length(y)~= length(x)
    error("x and y are not compatible");
end

% Least squares algorithm implementation
x = x(:); y = y(:); % x and y are column vectors
A = [x ones(size(x))]; % m-by-n matrix of the system
c = (A'*A)\(A'*y); % Solving normal equations

if nargout>1 % Checking number of function outputs
r = y - A*c;
R2 = 1 - (norm(r)/norm(y-mean(y)))^2;
end


Program Output:

>> x = [-6 -4 -1 0 3 5 8];
>> y = [18 13 6 4 -1 -8 -15];
>> linearregression(x,y)

ans =

   -2.3140

    4.0814

A MATLAB Program to Implement the Jacobi Iteration

A MATLAB Program to Implement Jacobi Iteration to Solve System of Linear Equations:

The following MATLAB codes uses Jacobi iteration formula to solve any system of linear equations where the coefficient matrix is diagonally dominant to achieve desired convergence.

%-----------------------------------------------------------------%
function [x, rel_error] = jacobimethod(A, b, x0, tol, iteration)

% Inputs: A - Coefficient matrix
%         b - Input matrix
%       tol - Defining tolerance for solution
% iteration - Number of iterations

% Outputs: x - Solutions
%  rel_error - Relative error

D = diag(diag(A));   % Making coefficient matrix diagonal
R = A - D;           % Construction of another matrix "R"
N = 1;               % iteration counter
x = x0;
rel_error = tol * 2; % norm(x - x0)/norm(x);

% Implementation of Jacobi method to solve Ax = b
while (rel_error>tol && N <= iteration)
   
    xprev = x;   
    x = inv(D)*(b - R*xprev);
    rel_error = norm(x - xprev)/norm(x);

    Z1=x(1);
    disp(Z1);
    Z2=x(10);
    disp(Z2);
    Z3=x(16);
    disp(Z3);
   
    fprintf('\n Iteration %i: Relative error =%d ',N, rel_error)   
    N = N + 1;       
end
%----------------------------------------------------------------%

Now, define your input matrices in the MATLAB command window and use the above developed function to compute Jacobi iteration.

MATLAB Command Prompt:

>> A = [-2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    1 -2 1 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 1 -2 1 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 1 -2 1 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 1 -2 1 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 1 -2 1 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 1 -2 1 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 1 -2 1 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 1 -2 1 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 1 -2 1 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 1 -2 1 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 1 -2 1 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 1 -2 1 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 1 -2 1 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 1 -2 1;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -2];

>> b = [-1; -2; -2; -2; -2; -2; -2; -2; -2; -2; -2; -2; -2; -2; -2; 1];

>> x0 = [0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0];

>> tol=10^-5;

>> iteration=350;

>> jacobimethod(A, b, x0, tol, iteration)


ans =

   14.8499
   28.7009
   40.5541
   50.4104
   58.2708
   64.1360
   68.0066
   69.8830
   69.7653
   67.6537
   63.5477
   57.4473
   49.3516
   39.2600
   27.1715
   13.0852

A MATLAB Program to Determine the Roots of Equation by Secant method

Problem: Solve the following problem by Secant Method.

f(x) = x - 2e^-x

Determine the roots of the above function where x1 = 0 and x2 = 1.


Solution: The following MATLAB program solves this problem.

% Definition of a function "secantmethod" to solve equation by Secant Method

function [x,ea] = secantmethod(X,X0,etol)

format long;

% Input: X=0, etol=Tolerance definition for error
% Output x=Solution of equation, ea=Calculated error in loop

% Program Initialization

% Iterative Calculations

while (1)
   
    % Implementation of Secant Method
         
    solution = X-((X-(2*exp(-X)))*(X0-X))/((X0-(2*exp(-X0)))-(X-(2*exp(-X))));    
   
    solutionprevious=X;
    X0=X;
    X=solution;
                 
if (solution-solutionprevious)~=0          % Approximate percent relative error
                                
   ea=abs((solution - solutionprevious)/solution)*100;
      
end

if ea<=etol                 
  
    break,               
end

%x0=solutionprevious;
x = solution;

%disp(x0);
disp(x);

end

% Display of output parameters

disp(x);                     

disp(ea);

%disp(solution-(0.1*abs(ea)*solution));

end


Program Outputs:

>> secantmethod(1,0,1e-10)

   0.883298154248460

   0.851584208519236

   0.852612692886131

   0.852605503703827

   0.852605502013723

   0.852605502013723

     3.125167797776517e-13


ans =

   0.852605502013723


#SecantMethod