PDE Numerical Solver Using Finite Differences
9 views (last 30 days)
Show older comments
Hey everyone, I'm working on the following problem:
Solve Laplace's equation on the heating 3 by 3 heating block with the boundary conditions of 75, 100, 50, and 0. This is to be done by using the Liebmann method with an over-relaxation factor of 1.5 and and a stopping criteria (relative error) of 1%.
So far, my code looks like this: ________________________________________________________________________________________________
% This will attempt to solve the Laplace equation with the given Dirichlet
% conditions.
% We will use the Liebmann Method, which will employ the usage of the
% following equation:
% T(i,j) = ( T(i+1,j) + T(i-1,j) + T(i,j+1) + T(i,j-1) )/4
% We will also employ overrelaxation
% T(i,j,new) = l*T(i,j,new) + (1-l)*T(i,j,old)
% Where l is the overrelaxation factor.
% Along with the error condition
%
% |e(a,i,j)| = (T(i,j,new) - T(i,j,old))
% __________________________ * 100%
% T(i,j,new)
% Boundary condition Matrix
T = zeros(5,5);
T(:,1) = 75;
T(:,5) = 50;
T(5,2) = 100;
T(5,3) = T(5,2);
T(5,4) = T(5,3);
T(1,1) = 75/2;
T(1,5) = 150/2;
T(5,5) = 50/2;
T(5,1) = 75/2;
l = 1.5; % weighting factor
e = 0.01; % desired error
for i = 2:4
for j = 2:4
T(i,j) = (T((i+1),j) + T((i-1),j) + T(i,(j-1)) + T(i,(j+1)))./4;
T(i,j) = l.*T(i,j);
end
end
disp(['The new matrix should look like'])
disp(T)
______________________________________________________________________________________________
This gets me as far as the first iteration. But to get the final solution the program requires I need about nine iterations, but how do I employ the overrelaxation and the error condition in this to generate the final values?
All help is greatly appreciated.
Thanks!
1 Comment
Alberto
on 17 Dec 2014
You should include your double for loop inside a while that controls the iteration, so it can be runned several times.
You should create an absolute error variable that calculates the difference beetween 2 diferent iterations (you should save the old one ) and use norm to determine the difference, and a tolerance variable.
Answers (2)
pshymn
on 8 Mar 2017
Edited: pshymn
on 8 Mar 2017
% Lecture: Sayisal Isi Gecisi ve Akis 2
% Source: Numerical method for engineers, Chapra&Canale, 6th edition
% Example 29.1
% Method used: Liebmann Method (Gauss-Seidel Method)
T_ust=100; %top-side temperature of heated plate
T_sol=75; %left-side temperature of heated plate
T_sag=50; %right-side temperature of heated plate
T_alt=0; %bottom-side temperature of heated plate
n=100; %m and n are the dimension of the matrices, so it can be any number such as 10, 20...
m=100;
L=1.5; %overrelaxation value
T=zeros(m,n,'double');
T1=zeros(m,n,'double');
for i=1:m %here, the up&bottom-sides are defined.
T(i,1)=T_alt;
T(i,n)=T_ust;
end
for j=1:n %again, left and right sides are defined.
T(1,j)=T_sol;
T(m,j)=T_sag;
end
for i=2:m-1 %the nods whose temperatures will be found are defined. they are nods in the middle.
for j=2:n-1
T(i,j)=50;
T(i,j)=50;
end
end
for k=1:1000 %this is the iteration number,you can keep it around 20-30.but to be sure it's 1000.
for i=2:m-1
for j=2:n-1
T1(i,j)=0.25*(T((i+1),j) + T((i-1),j) + T(i,(j-1)) + T(i,(j+1))); %formula1
T(i,j)=L*T1(i,j)+(1-L)*T(i,j); %formula2
end
end
end
%Draw a graph
x=0:1/(m-1):1;
y=0:1/(n-1):1;
t2d=transpose(T);
contourf(x,y,t2d, 30)
xlabel({'X [metre]'});
ylabel({'Y [metre]'});
title('TEMPERATURE DISTRIBUTION IN THE PLATE')
% at the end, there will be a .fig file which demonstrates temperature distribution in colorful window.
1 Comment
Ashley Turner
on 3 Dec 2019
What is the top boundary was a function? I can you incorporate that into the code..
i.e - The upper boundary: T=f(x)=ax-x^2
See Also
Categories
Find more on Mathematics and Optimization in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!