2D cylindrical heat equation

Tags:
1. Nov 9, 2014

range.rover

Hello friends, I am new for numerical methods and programming. i have been trying to devolop a program in 2D poisson heat equation in cylinder (r,angle) by finite difference method
2u/∂r2 + 1/r * ∂u/∂r + 1/r2 * ∂2u/∂θ2 = Q(u,θ)

discritized equation :-

ui+1,j − 2uij + ui−1,j/(∆r)2 + 1/ri * ui+1,j − ui−1,j/2∆r +1/r2i * ui,j+1−2uij +ui,j−1/(∆θ)2 = Q(u,θ)

i have been trying to develop this program, can anybody help me in making it to run faster for n = 400;

clear all
close all
clc
tic
n = 5; % mesh size
s = n-2;
r = linspace(1,5,n)-1;

angle =((1:1:n)-1)/(n)*2*pi;

m = (s*n); %number of nodes

dr = r(3)-r(2); % distance between radius

dangle = angle(3)- angle(2); % change of angle between radius

Q = 20*dr^2; % heat sourse

% right hand side Matrix

rhs = zeros(m+1,1);
rhs( : ) = Q;
bou = ones(m+1,1);

% boundary values
bou( : ) = 0;
tic
for i = n-1:n-2:m+1
rhs(i,1) = rhs(i,1) + bou(i);
end
toc

% evaluating the co effitiants of tridiagonal sparse matrix

tic
for i = 2:n-1
a(i) = 2*(1 + ((dr^2)/(r(i)*dangle)^2));
b(i) = 1+(1/(2*r(i)));
c(i) = 1-(1/(2*r(i)));
d(i) = dr^2/(r(i)^2*dangle^2);
end
toc
% building up the matrix for A.x = B

f = c(2);
I = speye(n);
a = (a(2:end))';
b = [0;b(2:n-2)'];
p = [c(3:end)';0];
% d = (d(2:end))';
T = spdiags([p -1*a b],[-1 0 1],n-2,n-2);
oops = full(T);
FirstM = kron(I,T);
% z = full(FirstM);
d = (d(2:end));
e = repmat(d,1,n-1);
tic
for i = 1: (m-s)
FirstM(i,(i+s)) = e(i);

FirstM((i+s),i) = e(i);
end
toc
tic
for i = 1: s
FirstM((m-s+i),i) = e(i);
FirstM(i,(m-s+i)) = e(i);
end
toc
% z = full(FirstM);

FirstM(1,m+1) = 0;
FirstM(m+1,1) = 0;
tic
for i = 1:n-2:m+1
FirstM(i,m+1) = f;
FirstM(m+1,i) = 1;
end
toc
FirstM(m+1,m+1) = -4;
z = full(FirstM);

% solving A*x = B for solution
Tsol = abs(mldivide(FirstM,rhs));
toc
T1soln = Tsol(1:m,1);
Tsoln = reshape(Tsol,n-1,n-1);
% Tsoln = Tsoln';

% Tsoln = reshape(T1soln,n,s);
% angle = angle(1,2:n);
% r = r(1,2:n);
[x,y] = meshgrid(r,angle);
[xs,ys] = pol2cart(x,y);
surf(xs,ys,Tsoln,'EdgeColor','none')

.

Last edited: Nov 9, 2014
2. Nov 14, 2014