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

### A Taylor Series Expansion to Approximate a Function near a Point

The Taylor series expansion up to seventh terms for approximation of near point x =1 is,

Letting starting point xₒ is at 0 and the step size (xi - xₒ) is 1. Now, the value of eˣ at x = 1 is to be calculated by Taylor series expansion on the basis of the starting point and its derivatives at x = 0 with step size 1. A MATLAB program is written to evaluate this value from zero order approximation to sixth order approximation and show the improvement of accuracy and minimization of truncation error as well.
The true value of is 2.718281828459046. The truncation error is the difference between true and approximate values. The following MATLAB program describes each steps of the calculation.

Taylor Series Implementation by a MATLAB Program

% Taylor Series expansion to approximate exp(x) near point x=1

function e = taylorseries(x0,x1,h)

e0=exp(x0);              % calculates exp(x) value when x=0

derivative1=1;           % all the six derivatives calculated at x=0
derivative2=1;
derivative3=1;
derivative4=1;
derivative5=1;
derivative6=1;

% zero order approximation

e=e0;
error=exp(x1)-e;
disp (e);        % display of exp(x) and error at command prompt
disp(error);

% first order approximation

e=e0+(derivative1*h);
error=exp(x1)-e;
disp (e);
disp(error);

% second order approximation

e=e0+(derivative1*h)+(derivative2*(h^2/factorial(2)));
error=exp(x1)-e;
disp (e);
disp(error);

% third order approximation
e=e0+(derivative1*h)+(derivative2*(h^2/factorial(2)))+(derivative3*(h^3/factorial(3)));
error=exp(x1)-e;
disp (e);
disp(error);

% fourth order approximation
e=e0+(derivative1*h)+(derivative2*(h^2/factorial(2)))+(derivative3*(h^3/factorial(3)))+(derivative4*(h^4/factorial(4)));
error=exp(x1)-e;
disp (e);
disp(error);

% fifth order approximation
e=e0+(derivative1*h)+(derivative2*(h^2/factorial(2)))+(derivative3*(h^3/factorial(3)))+(derivative4*(h^4/factorial(4)))+(derivative5*(h^5/factorial(5)));
error=exp(x1)-e;
disp (e);
disp(error);

% sixth order approximation
e=e0+(derivative1*h)+(derivative2*(h^2/factorial(2)))+(derivative3*(h^3/factorial(3)))+(derivative4*(h^4/factorial(4)))+(derivative5*(h^5/factorial(5)))+(derivative6*(h^6/factorial(6)));
error=exp(x1)-e;
disp (e);
disp(error);
end

Table 1: Outputs from the Above Program to Approximate at x=1

The above table tells us as the order of the Taylor series increases the truncation error also decreases. However, after sixth and higher order approximations, improvement of the accuracy is not much significant. So, we can stop the approximation to first seven terms of the Taylor series and the approximate value of ex near point x = 1 is 2.718.

### Finding Roots of a Equation by Newton-Raphson Method

In this tutorial, we are interested to the roots of the following function.

f(x) = sin x + 2 - 0.1 x

To find the value of x or roots of the equation, Newton-Raphson method is applied. For the proper initial guess for x, graphical method is applied to visualize the characteristics of this nonlinear function. This graphical technique also helps to find multiple roots for the unknown function. The following MATLAB commands generate the plot of the function, f(x). This function is plotted with respect to the values of x from -50 to +50.

>> X=-100:100;
>> plot(X,(sin(X)+2-(0.1*X)))

From the following figure 1, it is clear that the equation has seven roots and plotting the function helps to get the roots. The function changes its sign from, x = 10 - 11, 11 -12, 16 -17, 18 - 19, 21 - 22, 25 - 26, and 27 - 28. And, the roots of the equation must lie between these intervals. Now, to get good approximations of the roots, a Newton-Raphson algorithm is written in MATLAB.

Figure 1: Plot of the function f(x) to identify the intervals where roots of the function lie.

The program is developed by MATLAB user defined function called “sinusoidal”. The initial assumptions for approximating seven roots of x are 10.5, 11.5, 16.5, 18.5, 21.5, 25.5 and 27.5 respectively. The codes are given below. It is straightforward to run the program seven times as we have seven initial guesses for the corresponding roots.

## Newton-Raphson Algorithm by MATLAB

% Definition of a function "sinusoidal" to solve equation by Newton-Raphson
function [x,ea] = sinusoidal(X, etol)
format long;

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

% Program Initialization

solution = 21.5;

% Iterative Calculations

while (1)
solutionprevious = solution;

% Implementation of Newton-Raphson Method

solution = X - (sin(X)+2-(0.1*X))/(cos(X)-0.1);
X=solution;

if solution~=0                         % Approximate percent relative error

ea=abs((solution - solutionprevious)/solution)*100;
end

if ea<=etol,                  % Condition to meet specified error tolerance

break,
end

x = solution;

disp(x);

end

% Display of output parameters

disp(x);
disp(ea);

end

## Program Outputs:

>> sinusoidal (10.5,1e-10)
ans =
10.636782424281707
>> sinusoidal(11.5, 1e-10)
ans =
11.562055865585652
>> sinusoidal (16.5,1e-10)
ans =
16.107752970689450
>> sinusoidal (18.5,1e-10)
ans =
18.721338781142521
>> sinusoidal (21.5,1e-10)
ans =
21.809224297055675
>> sinusoidal (25.5,1e-10)
ans =
25.744697403517950
>> sinusoidal (27.5,1e-10)
ans =
27.435908953388580

So, the seven roots from the program are 10.637, 11.562, 16.108, 18.721, 21.809, 25.745 and 27.436. In all the cases, the error tolerance limit is kept extremely small (0.0000000001) to get more accurate and precise results.

#Blog #Blogger #NewtonRaphson #Matlab #Sinusoidal

### Explicit Euler Method to Solve System of ODEs in MATLAB

In this tutorial, I am going to show a simple way to solve system of first order ordinary differential equations (ODE) by using explicit Euler method. Let' say, we have following three first order ODEs.

Here, we have 3 ODEs, 3 dependent variables (x, y, z), and 1 independent variable, t. The initial conditions when time is zero are, x(0) = y(0) = z(0) =1. Our goal is to solve these differential equations with explicit Euler approach and plot the solutions afterwards.

Step 1: Define the Equations
The first step is to define all the differential equations in MATLAB. I did this by using MATLAB function handle, which is shown below.

Step 2: Choose a Numerical Approach
The next step is to select a numerical method to solve the differential equations. In this example, we will use explicit Euler method. I have created a function to implement the algorithm. The following image shows the application of the explicit Euler method.

Step 3: Call the Function
This is the final step where we need to call the function. In the previous step, we have created a function, which we are going to call in this step to solve the equations.

So, that's it. We are done with the process, and we can now visualize the solution. I have also attached the MATLAB codes for this problem at the end.

MATLAB Program:

close all;

clc;
format long;

% Defining the three differential equations of the problem
f = @(t,y) [-5*y(1)+5*y(2); 14*y(1)-2*y(2)-y(1)*y(3); -3*y(3)+y(1)*y(2)];
[x,y] = explicit_euler(f,[0,5],[1;1;1],0.001); % Calling Euler function
plot(x,y);
title('When Time Step is 0.001');
legend('x(t)', 'y(t)', 'z(t)', 'Location', 'NorthEast')
xlabel('t')
ylabel('Solutions')

function [x, y] = explicit_euler( f, xRange, y_initial, h )
% This function uses Euler’s explicit method to solve the ODE
% dv/dt=f(t,v); x refers to independent and y refers to dependent variables
% f defines the differential equation of the problem
% xRange = [x1, x2] where the solution is sought on
% y_initial = column vector of initial values for y at x1
% numSteps = number of equally-sized steps to take from x1 to x2
% x = row vector of values of x
% y = matrix whose k-th column is the approximate solution at x(k)
x(1) = xRange(1);
numSteps = ( xRange(2) - xRange(1) ) /h ;
y(:,1) = y_initial(:);
for k = 1 : numSteps
x(k + 1) = x(k) + h;
y(:,k+1) = y(:,k) + h * f( x(k), y(:,k) );
end

#ExplicitEulerMethod #Matlab #NumericalMethod #Blog #Blogger

### Discussion on Numerical Integration Approaches with MATLAB

In this tutorial, we are going to implement four numerical integration schemes in MATLAB. The methods are:

i)   Midpoint Rectangle Rule
ii)  Trapezoidal Rule
iii) Simpson’s 1/3 Rule
iv) Simpson’s 3/8 Rule

Let's say, we would like to integrate the following error function where the upper interval, a is 1.5. At first, we will evaluate the true integral of the error function, and then, we will find the integral numerically with the above methods (i) -(iv).

Finally, we will compare each method with the true value and make some comments based on the deviations if we have any. For each method, we will find out the following parameters:

Here, n refers to the numbers, 1,2,4,8,16, ...,128. I_num is the numerical value of integration, h is the step size, and E is the error, which is the absolute difference between true and numerical results. First, we are going to determine the true value of the error function while integrating from 0 to 1.5. Then, we will apply the corresponding methods and compare the absolute error.

MATLAB Code for Midpoint Rectangle Rule

close all;
clear;
clc;

% Defining the error function
f = @(x) (2/sqrt(pi))*exp(-x.^2);
a = 0;
b = 1.5; % Upper integration limit
n = [1 2 4 8 16 32 64 128];
L = length(n);

% Exact integral using MATLAB built-in function
I = integral(f,a,b);

fprintf('\n     n          h         I_NUM          E')

for j = 1:L
h(j)=(b-a)/n(j);

x_1 = a:h(j):b;
% Using MATLAB built-in command for repeated arrays
x_2 = repelem(x_1(2:end-1),2);
xh = [a x_2 b];

% Numerical integration by mid-point rectangular rule
I_NUM=0;
for i = 1:length(x_1)-1
I_NUM = I_NUM + f(a  +(i-1/2).*h).*h;
end

% Defining the error
E(j) = abs(I-I_NUM(j));

fprintf('\n %5.6g %11.6g %12.6g %10.6g', n(j), h(j), I_NUM(j), E(j))
fprintf('\n')

end

The following table shows the parameters determined by this approach:

MATLAB Code for Trapezoidal Rule

close all;
clear;
clc;

% Defining the error function
f = @(x) (2/sqrt(pi))*exp(-x.^2);
a = 0;
b = 1.5; % Upper integration limit
n = [1 2 4 8 16 32 64 128];
L = length(n);

% Exact integral using MATLAB built-in function
I = integral(f,a,b);

fprintf('\n     n          h         I_NUM          E')

for j = 1:L
h(j)=(b-a)/n(j);

% Numerical integration using Trapezoidal rule
I_NUM=0;
for i = 1:n(j)
I_NUM = I_NUM+(f(a+(i-1).*h) + f(a+i.*h)).*h/2;
end

% Defining the error
E(j) = I - I_NUM(j);

fprintf('\n %5.6g %11.6g %12.6g %10.6g', n(j), h(j), I_NUM(j), E(j))
fprintf('\n')
end

The following table shows the parameters determined by this approach:

MATLAB Code for Simpson’s 1/3 Rule

close all;
clear;
clc;

% Defining the error function
f = @(x) (2/sqrt(pi))*exp(-x.^2);
a = 0;
b = 1.5; % Upper integration limit
n = [1 2 4 8 16 32 64 128];
L = length(n);

% Exact integral using MATLAB built-in function
I = integral(f,a,b);

fprintf('\n     n          h         I_NUM          E')

for j = 1:L
h(j) = (b-a)/(2*n(j));

% Numerical integration using Simpson’s 1/3 rule

I_NUM = 0;

for i = 1:n(j)
I_NUM = I_NUM+(f(a+(2*i-2).*h) + 4*f(a+(2*i-1).*h) + f(a+(2*i).*h)).*h./3;
end
E(j) = abs(I - I_NUM(j));

fprintf('\n %5.6g %11.6g %12.6g %10.6g', n(j), h(j), I_NUM(j), E(j))
fprintf('\n')
end

The following table shows the parameters determined by this approach:

MATLAB Code for Simpson’s 3/8 Rule

close all;
clear;
clc;

% Defining the error function
f = @(x) (2/sqrt(pi))*exp(-x.^2);
a = 0;
b = 1.5; % Upper integration limit
n = [1 2 4 8 16 32 64 128];
L = length(n);

% Exact integral using MATLAB built-in function
I = integral(f,a,b);

fprintf('\n     n          h         I_NUM          E')

for j = 1:L
h(j) = (b-a)/(3*n(j));

% Numerical integration using Simpson’s 3/8 rule
I_NUM = 0;
for i=1:n(j)
I_NUM = I_NUM + (f(a+(3*i-3).*h) + 3*f(a+(3*i-2).*h) + 3*f(a+(3*i-1).*h) + f(a+(3*i).*h))*3.*h/8;
end

% Defining the error
E(j) = abs(I - I_NUM(j));

fprintf('\n %5.6g %11.6g %12.6g %10.6g', n(j), h(j), I_NUM(j), E(j))
fprintf('\n')

end

The following table shows the parameters determined by this approach:

We see from the above comparisons that both Simpson's rules give us a good approximation of the true value, especially when n is greater than 4. For the rectangular approach, the value reaches close to the true value when n is 128. The same is also evident for the trapezoidal rule. So, we can conclude that both the Simpson's rules are more efficient and accurate approaches compared to the rest.

#TrapezoidalRule #SimpsonsRule #Matlab #NumericalIntegration #Blog #Blogger

### Stress and Strain Analysis of Linear Finite Elements

In this tutorial, we will use principle of virtual work method to determine the stress and displacement of the nodes of linear one dimensional finite elements. Following figure 1 represents a cantilever beam which is stressed axially. Let's assume the parameters L, A, E, c have the same value of 1 unit for simplicity. Our goal is to compare the finite element results with the exact solutions of a typical cantilever problem.
For this problem, we need to consider two quadratic one dimensional element. So, the displacement is interpolated among the three nodes and has the following form:

 Fig.1: Showing one dimensional axially load bar given in the assignment problem.

 Fig.2: Showing one dimensional piecewise linear interpolation.

## Mathematica Program for Element 1

Clear[x, a0, a1, a2, u1, u2, u3, L]

(*Definition of quadratic form of displacement*)
u = a0 + a1*x + a2*x^2;

eq1 = (u /. x -> 0) - u1;
eq2 = (u /. x -> L/4) - u2;
eq3 = (u /. x -> L/2) - u3;

sol = Solve[{eq1 == 0, eq2 == 0, eq3 == 0}, {a0, a1, a2}];
u = u /. sol[[1]]

N1 = Coefficient[u, u1];
N2 = Coefficient[u, u2];
N3 = Coefficient[u, u3];

(*Shape Functions*)
Nn = {N1, N2, N3}

p = c*x;

(*Nodal Force Vectors*)
fe = Integrate[Nn*p, {x, 0, L/2}];
c = 1;
L = 1;
fe // MatrixForm

### Program Outputs

Displacement Vector for Element 1:

Shape Functions for Element 1:

Nodal Force Vector for Element 1:

## Mathematica Program for Element 2

Clear[x, a0, a1, a2, u1, u2, u3, L]

(*Definition of quadratic form of displacement*)
u = a0 + a1*x + a2*x^2;

eq1 = (u /. x -> 0) - u1;
eq2 = (u /. x -> L/4) - u2;
eq3 = (u /. x -> L/2) - u3;

sol = Solve[{eq1 == 0, eq2 == 0, eq3 == 0}, {a0, a1, a2}];
u = u /. sol[[1]]

N1 = Coefficient[u, u1];
N2 = Coefficient[u, u2];
N3 = Coefficient[u, u3];

(*Shape Functions*)
Nn = {N1, N2, N3}

p = c*(x + L/2);

(*Nodal Force Vectors*)
fe = Integrate[Nn*p, {x, 0, L/2}];
c = 1;
L = 1;
fe // MatrixForm

### Program Outputs

Displacement Vector for Element 2

Shape Functions for Element 2

## Local Stiffness Matrix for Element 1

Clear[x, a0, a1, a2, u1, u2, u3, L, EA]

(*Definition of quadratic form of displacement*)
u = a0 + a1*x + a2*x^2;

eq1 = (u /. x -> 0) - u1;
eq2 = (u /. x -> L/4) - u2;
eq3 = (u /. x -> L/2) - u3;

sol = Solve[{eq1 == 0, eq2 == 0, eq3 == 0}, {a0, a1, a2}];
u = u /. sol[[1]];

N1 = Coefficient[u, u1];
N2 = Coefficient[u, u2];
N3 = Coefficient[u, u3];

B = {{D[N1, x], D[N2, x], D[N3, x]}};

(*Local Stiffness Matrix*)
K = Integrate[EA*Transpose[B].B, {x, 0, L/2}];

L = 1;
EA = 1;
K // MatrixForm

## Comparison between the Exact and Solution using Four Elements

We know that the exact solution for the axially loaded bar shown in Fig. 1 may be written as,

Where, for the consecutive four elements, we may write the displacement for each element as,

From the assignment problem, the parameter values are,

From the virtual work method, we may write,

The internal virtual work (IVW) may be expressed as,

And, the external virtual work may be presented as,

## MATHEMATICA PROGRAM

The following program is written to determine the nodal displacement using the finite element method considering the principle of virtual work.

Clear[x, a0, a1, a2, u1, u2, u3, L]

(*Definition of quadratic form of displacement*)
u = a0 + a1*x + a2*x^2;

eq1 = (u /. x -> 0) - u1;
eq2 = (u /. x -> L/4) - u2;
eq3 = (u /. x -> L/2) - u3;

sol = Solve[{eq1 == 0, eq2 == 0, eq3 == 0}, {a0, a1, a2}];
u = u /. sol[[1]];

(*Shape Functions*)
N1 = Coefficient[u, u1];
N2 = Coefficient[u, u2];
N3 = Coefficient[u, u3];
Nn = {N1, N2, N3};

B = {{D[N1, x], D[N2, x], D[N3, x]}};

c = 1;
L = 1;
EA = 1;

(*Axial Force on Element 1*)
p1 = c*x;
fe1 = Integrate[Nn*p1, {x, 0, L/2}];

(*Axial Force on Element 2*)
p2 = c*(x + L/2);
fe2 = Integrate[Nn*p2, {x, 0, L/2}];

K1 = Integrate[EA*Transpose[B].B, {x, 0, L/2}];

K2 = Integrate[EA*Transpose[B].B, {x, 0, L/2}];

Ndof = 5;
Ne = 2;
ID1 = {1, 2, 3};
ID2 = {3, 4, 5};

K = ConstantArray[0, {Ndof, Ndof}];
K[[ID1, ID1]] = K1 + K[[ID1, ID1]];
K[[ID2, ID2]] = K2 + K[[ID2, ID2]];

fe = ConstantArray[0, {Ndof, 1}];
fe[[ID1]] = fe1 + fe[[ID1]];
fe[[ID2]] = fe2 + fe[[ID2]];

IDE = {1};
IDN = {2, 3, 4, 5};

KR = Delete[K, IDE];
KR = Delete[Transpose[KR], IDE];
KR = Transpose[KR];

feR = Delete[fe, IDE];

(*Calculation of Nodal Displacements*)
DisR = LinearSolve[KR, feR];

Dis = ConstantArray[0, {Ndof, 1}];
Dis[[IDN]] = DisR + Dis[[IDN]];
Dis // MatrixForm

### Program Outputs

So, we see that the exact and finite element solutions are the same.

From the following two plots, we notice that although the nodal displacement is continuous from the FEA result, it is not true for the stress. The stress distribution is indeed discontinuous.

## MATHEMATICA PROGRAM FOR COMPARISON

(*Exact Solution*)
uexact = c*L^2/2/EA*x - c/6/EA*x^3;

(*Finite Element Approximation*)
N1 = Piecewise[{{4 x/L, 0 <= x <= L/4}, {4/L (L/2 - x),
L/4 <= x <= L/2}}, 0];
N2 = Piecewise[{{4/L (x - L/4), L/4 <= x <= L/2}, {4/L (3 L/4 - x),
L/2 <= x <= 3 L/4}}, 0];
N3 = Piecewise[{{4/L (x - L/2), L/2 <= x <= 3 L/4}, {4/L (L - x),
3 L/4 <= x <= L}}, 0];
N4 = Piecewise[{{4/L (x - 3 L/4), 3 L/4 <= x <= L}}, 0];

u1 = c*L^3/EA*(47/384);
u2 = c*L^3/EA*(11/48);
u3 = c*L^3/EA*(39/128);
u4 = c*L^3/EA*(1/3);

(*Calculating Nodal Displacement*)
u = u1*N1 + u2*N2 + u3*N3 + u4*N4;
uexactp = uexact /. {c -> 1, EA -> 1, L -> 1}

sigmaexact = D[uexactp, x]
up = u /. {c -> 1, EA -> 1, L -> 1}
sigmap = D[up, x]

(*Plotting to Show the Comparison between Exact and FEA Solution*)
Plot[{uexactp, up}, {x, 0, 1},
AxesLabel -> {"x", "Nodal Displacement"},
PlotLegends -> {"uexact", "u_Finite Element"}]
Plot[{sigmaexact, sigmap}, {x, 0, 1}, AxesLabel -> {"x", "sigma_11"},
PlotLegends -> {"sigma_exact", "sigma_Finite Element"}]