## 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 Short Introduction to MATLAB for Beginners

In this post, you will find series of basic tutorial videos on MATLAB. As you may know that MATLAB is a very powerful computational tool, which has been widely used in academics as well as industries. You may find these videos useful if you haven't yet used this tool.

### 2D Laplace Equation Solution by 5 Point Finite Difference Approximation

The temperature distribution of a rectangular plate is described by the following two dimensional (2D) Laplace equation:

Txx + Tyy = 0

The width (w), height (h), and thickness (t) of the plate are 10, 15, 1 cm, respectively. There is a heat source at the top edge, which is described as, T = 100 sin (πx / w) Celsius, and all other three edges are kept at 00 C. There is no heat flow along the plate thickness since the faces are insulated. Also, assume that there is no internal heat generation so that the problem becomes 2D temperature distribution.

Solve the Laplace equation with a 5 point finite difference approximation by MATLAB where Δx = Δy = 2.5 cm. The temperature distribution is defined by T(x,y)

The following MATLAB program develops a 5-point approximation stencil to solve the 2D Laplace equation.

MATLAB Program:

close all;
clc;
% Width, height & thickness of the thin plate (dimension, cm)
W=10.0;
h=15.0;
t=1.0;
% Gird size definition along x and y axes
xx=2.5;
yy=2.5;

imax=(W/xx)-1; % Array dimension along the width
jmax=(h/yy)-1; % Array dimension along the height
N_array=imax*jmax; % Total number of grid points

% Generating the Matrix, A
for j=1:N_array
for i=1:N_array
if i==j
a(i,j)=-4.0;
end
if j==(i+1)
a(i,j)=1.0;
end
if j==(i-1)
a(i,j)=1.0;
end
if j==(i+5)
a(i,j)=1.0;
end
if j==(i-5)
a(i,j)=1.0;
end
if i==(imax+2)
a(i,i+1)=0.0;
end
if i==(imax+3)
a(i,i-1)=0.0;
end
if i==(imax+2+imax+2)
a(i,i+1)=0.0;
end
if i==(imax+2+imax+2+1)
a(i,i-1)=0.0;
end
end
end

for i=1:imax
x(i)=(i)*xx;
Tb(i)=100*sin(pi*x(i)/W); % Top edge temperature (in Celcius)
end

% Generating the Matrix, b
k=imax+2;
for i=1:imax
b(k)=-Tb(i);
k=k+imax+2;
end

%Displaying A and b matrices
B=b';
T=GaussElimination(a,B); % Calling the function defined elsewhere
TT(:,1) = T(1:5);
TT(:,2) = T(6:10);
TT(:,3) = T(11:15);
TT = [TT;Tb];

imagesc([0:xx:xx*imax],[0:yy:yy*jmax],TT)
h=colorbar;
colormap jet;
ylabel(h,'Temperature, C')
set(gca,'YDir','normal')
xlabel('Width, cm')
ylabel('Height, cm')

The following function is created in MATLAB to solve the system of algebraic equations using Gauss elimination method:

function x = GaussElimination(A,b)
% Gauss Elimination without Pivoting
% x = GaussElimination(A,b)
% input:
% A = Coefficient Matrix
% b = Right Hand Side Vector
% output:
% x = Solution Vector

[m,n] = size(A);
if m~=n
error('Matrix A must be square');
end
nb = n+1;
ATM = [A b];

% Forward Elimination
for k = 1:n-1
for i = k+1:n
factor = ATM(i,k)/ATM(k,k);
ATM(i,k:nb) = ATM(i,k:nb)-factor*ATM(k,k:nb);
end
end

% Backward Substitution
x = zeros(n,1);
x(n) = ATM(n,nb)/ATM(n,n);
for i = n-1:-1:1
x(i) = (ATM(i,nb)-ATM(i,i+1:n)*x(i+1:n))/ATM(i,i);
end

disp('Solution Vector')
disp(x)
disp('Coefficient Matrix')
a1 = ATM(:,1:n);
disp(a1)
disp('The b Vector')
b1 = ATM(:,nb);
disp(b1)
disp('Verification of the Solution')
C = a1*x-b1;
disp(C)

Program Output:

Solution Vector
1.3046
3.3734
7.4183
15.8087
33.4596
1.8450
4.7707
10.4910
22.3568
47.3190
1.3046
3.3734
7.4183
15.8087
33.4596

Coefficient Matrix
Columns 1 through 11

-4.0000    1.0000         0         0         0    1.0000         0         0         0         0         0
0   -3.7500    1.0000         0         0    0.2500    1.0000         0         0         0         0
0         0   -3.7333    1.0000         0    0.0667    0.2667    1.0000         0         0         0
0         0         0   -3.7321    1.0000    0.0179    0.0714    0.2679    1.0000         0         0
0         0         0         0   -3.7321    0.0048    0.0191    0.0718    0.2679    1.0000         0
0         0         0         0         0   -3.7321    1.0718    0.0192    0.0051    0.0013    1.0000
0         0         0         0         0         0   -3.4050    1.0824    0.0220    0.0055    0.2872
0         0         0         0         0         0         0   -3.3673    1.0839    0.0210    0.0964
0         0         0         0         0         0         0         0   -3.3638    1.0786    0.0343
0         0         0         0         0         0         0         0         0   -3.3861    0.0124
0         0         0         0         0         0         0         0         0         0   -3.7047
0         0         0         0         0         0         0         0         0         0         0
0         0         0         0         0         0         0         0         0         0         0
0         0         0         0         0         0         0         0         0         0         0
0         0         0         0         0         0         0         0         0         0         0

Columns 12 through 15

0         0         0         0
0         0         0         0
0         0         0         0
0         0         0         0
0         0         0         0
0         0         0         0
1.0000         0         0         0
0.3179    1.0000         0         0
0.1088    0.3219    1.0000         0
0.0385    0.1094    0.3206    1.0000
1.0947    0.0323    0.0114    0.0037
-3.3489    1.1156    0.0393    0.0124
0   -3.2968    1.1193    0.0365
0         0   -3.2919    1.1072
0         0         0   -3.3318

The b Vector
0
0
0
0
-70.7107
-0.0907
-0.3887
-1.4838
-5.5569
-120.7386
-0.5983
-1.9828
-5.5408
-14.9919
-111.4802

Verification of the Solution
1.0e-13 *

-0.0022
0
0
-0.0711
0.1421
0.0053
-0.0078
-0.0644
-0.1421
0
-0.0011
-0.0044
0.0089
0.0355
0

The following plot shows the desired temperature distribution within the plate.

### Hyperbolic Equation Solution by a Second Order Upwind Approximation by MATLAB

Let's assume that the temperature distribution inside a pipe is governed by the following one dimensional unsteady hyperbolic partial-differential-equation

Tt + u Tx = 0

Consider negligible heat diffusion and the initial velocity u is, 0.1 cm/s. The boundary conditions are described below:

T(x, 0) = 200 x;       0 ≤ x ≤ 0.5

T(x, 0) = 200 (1 – x);       0.5 ≤ x ≤ 1

Develop an algorithm in MATLAB to solve this problem using the second-order approximation in time and the second-order upwind approximation in space. Consider, Δx = 0.05 cm; Δt = 0.05.

The following MATLAB program develops the second order upwind method to solve the given one-dimensional unsteady hyperbolic convection equation.

% Hyperbolic equation - Convection equation: Upwind Method
close all;
clc;

% Input Properties from Problem
L=1.0; % (dimension in cm)
Lmax=5.0;
alpha=0.01;
u=0.1; % (dimension in cm/s)
Tt=10; % Simulation Time Period

% Mesh/Grid size
dx=0.05;
dt=0.05;

c=u*dt/dx; % Depends on dt values
d=alpha*((dt/dx)^2.0);

% Matrix Parameters

imax=(Lmax/dx)+1; % Array dimension along the width
tmax=(Tt/dt); % Array dimension along the height
ndim=imax-2;
tdim=tmax-2;

% Total no of iterations
iteration=10000;

% Convergence Criterion
tolerance=0.000001;

% Array and Variables
% a = Matrix of Coefficients A(i,j)
% b = Right Side Vector, b(i)
% x = Solution Vector, x(i)

% Initial Condition
f(:,1)=0.0;
for i=1:imax-2
x(i)=(i-1)*dx;

if x(i)<=1.0
f(i,1)=0.0;
end
if x(i)>1.0 && x(i)<=1.5
f(i,1)=200*(x(i)-1.0);
end
if x(i)>1.5 && x(i)<=2.0
f(i,1)=200*(L-x(i)+1.0);
end
if x(i)>2.0
f(i,1)=0;
end
end

fprintf('imax = %f\n',imax);
fprintf('c=%f\n',c);

% Upwind Method

for t=1:tmax
time(t+1)=(t)*dt;
for i=2:imax-2
% Second Order Approximation
f(i,t+1)=f(i,t)-c*(f(i,t)-f(i-1,t))-0.5*c*(1.0-c)*(f(i,t)-2.0*f(i-1,t)+f(i-2,t));
end
end

hold on
plot(x,f(:,1))
plot(x,f(:,find(time==1)),'o-')
plot(x,f(:,find(time==5)),'o-')
plot(x,f(:,find(time==10)),'o-')

xlabel('x, cm')
ylabel('Temperature, C')
legend('t=0s','t=1s','t=5s','t=10s')
title ('c=0.1')

## Program Output:

Second-order approximation in time and second-order upwind approximation in space: