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

Leapfrog Method with RK2 as a Start-up Scheme in MATLAB

Let's consider the following differential equation:

y'' + 10 y' + 500y = 0
Where, the initial conditions are,

y(0) = -0.025, and y'(0) = -1

We will solve this differential equation using a multi-step method, Leapfrog, where second order Runge-Kutta approach (RK2) is added as a start-up scheme in the algorithm. The following MATLAB program implements the Leapfrog method with initialization with RK2 method.

MATLAB Program:

close all;
clc;
h = 0.02;     % Step size
Tmax = 0.5;    % Maximum time
N = Tmax / h;  % Maximum number of steps
alpha=0.5;
t = linspace(0,0.5,N+1);  % Time range

% Analytical solution of the differential equation
y_real = -(((9.*sqrt(19))/760).*exp(-5.*t).*sin(5.*sqrt(19).*t)) - ((1/40).*exp(-5.*t).*cos(5.*sqrt(19).*t));
plot(t,y_real);
hold on

%Numerical solution
f=@(t,y) [y(2); -500*y(1)-10*y(2)]; % Governing system of equations

% Initial Conditions
Y = [-0.025; -1];
% Initialization with second order Runge-Kutta method
k1 = h.*f(t(1),Y(:,1));
    k2 = h.*f(t(1)+alpha.*h, Y(:,1)+alpha.*k1);
    Y(:,2) = Y(:,1) + (1-1/2/alpha).*k1 + k2/2/alpha;

% Leapfrog method steps
for i=2:N
    Y(:,i+1) = Y(:,i-1) + 2.*h.*f(t(i),Y(:,i));
end
plot(t,Y(1,:),'o:')
legend('Exact Solution','Leapfrog Solution','Location','NorthEast')
title('When h = 0.001')
xlabel('t')

ylabel('y')


Program Output:

The following plot shows the numerical and analytical solution of the differential equation with respect to time.

Numerical and Analytical Solution with Time

Implicit Euler Method by MATLAB to Solve an ODE

In this example, an implementation of the Implicit Euler approach by MATLAB program to solve an ordinary differential equation (ODE) is presented. Let's consider a differential equation, which is defined as,

dv/dt = p(t) v + q(t)

Where,

p(t) = 5(1+t)

and,

q(t) = (1+t)e-t

The initial value is, v(0) = 1; and the time period is 0 < t < 10.


Implicit Euler approach is unconditionally stable. The implementation of Implicit Euler scheme may be represented as,

v_n+1 = (v_n + hq_n+1) / (1 + hp_n+1)

The following MATLAB program implements this scheme:

MATLAB program:

close all;
clc;
y1(1) = 0; % Initial condition
h = 0.1; % Time step
t1(1) = 0;
i = 1;

% Implementation of the Implicit Euler approach
while (t1(i) <= 30)
    q=(1+t1(i))*exp(-t1(i));
    p= 5*(1+t1(i));
       
    y1(i+1) = (y1(i)+(h*q))/(1+h*p);
    t1(i+1) = t1(i) + h;
    i = i + 1;
end
plot(t1,y1,'-o')
hold on


Program Output:
The following plot shows the progression of the solution as the time increases.
Solution with respect to Time

Explicit Euler Method by MATLAB to Solve an ODE

In this example, an implementation of the Explicit Euler approach by MATLAB program to solve an ordinary differential equation (ODE) is presented. Let's consider a differential equation, which is defined as,

dv/dt = p(t) v + q(t)

Where,

p(t) = 5(1+t)

and,

q(t) = (1+t)e^-t

The initial value is, v(0) = 1;, and the time period is 0 < t < 10.


The implementation of Explicit Euler scheme may be represented as,

v_n+1 = v_n (1 - ph) + hq 


The following MATLAB program implements this scheme:


close all;
clc;
y_initial = 1; % Defines initial condition, v(0)=1
h = 1.1; % Defines time step
[ t, v ] = expliciteuler( @f_ode, [ 0.0,10 ], y_initial, h );
plot(t,v,'--*')
hold on


function f_vt = f_ode ( t, v )
% This function defines the differential equation
% t is the independent variable
% v is the dependent variable
% f_vt represents the dv/dt
p= 5*(1+t);
q=(1+t)*exp(-t);
f_vt = p*v+q;


function [x, y] = expliciteuler( f_ode, xRange, y_initial, h )
% This function uses Euler’s explicit method to solve the ODE
% dv/dt=f_ode(t,v); x refers to t and y refers to v
% f_ode 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(1,k+1) = x(1,k) + h;
y(:,k+1) = y(:,k) + h * f_ode( x(k), y(:,k) );

end

Program Output:

The following plot shows the value of the solution with respect of time increment.

time versus solution plot

Learning Mathematica, Lesson 3: Integration

Integration becomes very handy if you choose to use Wolfram Mathematica. You can evaluate definite and indefinite integrals using two ways.

1. You may use "Integrate[f, x]" command for the indefinite integrals, and "Integrate[f, {x, upper limit, lower limit}]" for the definite integrals.

2. Use int to enter  and then use  to enter the lower limit, then  for the upper limit:

For example, if you would like to find out the following integral of a function using both approaches,


Integration in Mathematica

Now, let's look at another example, where you will find detailed Mathematica codes to determine the integral of the following:

integral example







Note that the following Mathematica program also implements a numerical integration scheme by using a built-in "NIntregate" command. Look carefully the systaxes used in the following codes to display the results, you may represent them in your own fashion.


Mathematica Codes:
(*Showing Exact Integration of First Problem*)
Print["Exact solution of \!\(\*SubsuperscriptBox[\(\[Integral]\), \
\(-1\), \(1\)]\)(2\!\(\*SuperscriptBox[\(x\), \(2\)]\)+3)\
\[DifferentialD]x"]
FullSimplify[\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(-1\), \(1\)]\(\((2
\*SuperscriptBox[\(x\), \(2\)] + 3)\) \[DifferentialD]x\)\)]

(*Showing Numerical Integration of First Problem*)
Print["Numerical Integration by Approximation: \
\!\(\*SubsuperscriptBox[\(\[Integral]\), \(-1\), \
\(1\)]\)(2\!\(\*SuperscriptBox[\(x\), \(2\)]\)+3)\[DifferentialD]x"]
NIntegrate[(2 x^2 + 3), {x, -1, 1}]
Print["Numerical Integration by Gaussian Rule Order-1: \
\!\(\*SubsuperscriptBox[\(\[Integral]\), \(-1\), \
\(1\)]\)(2\!\(\*SuperscriptBox[\(x\), \(2\)]\)+3)\[DifferentialD]x"]
a = 2 x^2 + 3 /. x -> -Sqrt[1/(3)];
b = 2 x^2 + 3 /. x -> Sqrt[1/(3)];

c = N[a + b]



Program Output:
Program Output



A Mathematica Program for the Newton Raphson Method

The following Mathematica program implements the widely used Newton Raphson approach. Any functions may be considered here, and the following program is ready to evaluate it.

(*Mathematica Codes for Newton-Raphson Approach*)
(*Write Your Own Function Below*)
Sigma = 100000*(Sqrt[epsilon] + 1/80*Sin[Pi*epsilon/0.002]) - 4000;
SigmaDerivative = D[Sigma, epsilon];
epsilontable = {0.0002};
Error = {1};
SetTolerance = 0.0005;
MaximumIteration = 100;
i = 1;

(*Newton-Raphson Algorithm Implementation*)
While[And[i <= MaximumIteration, Abs[Error[[i]]] > SetTolerance],
  epsilonnew = 
   epsilontable[[i]] - (Sigma /. 
       epsilon -> epsilontable[[i]])/(SigmaDerivative /. 
       epsilon -> epsilontable[[i]]); 
  epsilontable = Append[epsilontable, epsilonnew]; 
  Errornew = (epsilonnew - epsilontable[[i]])/epsilonnew; 
  Error = Append[Error, Errornew];
  i++];

L = Length[epsilontable];
SolutionTable = 
  Table[{i - 1, epsilontable[[i]], Error[[i]]}, {i, 1, L}];
SolutionTable1 = {"Iteration Number", "Epsilon", "Error"};
L = Prepend[SolutionTable, SolutionTable1];
Print["Showing Results for the Newton-Raphson Approach when Sigma is \
4,000 and Initial Guess is 0.0002"]
ScientificForm[L // MatrixForm, 4]


Program Output:

Program Output

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