PDE Numerical Solver Using Finite Differences

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!

[Alberto]

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.

Answer 1

[pshymn]

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

[Ashley Turner]

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

Answer 2

[pshymn]