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.
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 ...
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).
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];
TT =
padarray(TT,[1 0],'pre');
TT =
padarray(TT, [0 1]);
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
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
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:
