Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

2D cylindrical heat equation

  1. Nov 9, 2014 #1
    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. jcsd
  3. Nov 14, 2014 #2
    Thanks for the post! This is an automated courtesy bump. Sorry you aren't generating responses at the moment. Do you have any further information, come to any new conclusions or is it possible to reword the post?
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook