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

Convergence of my code

  1. Aug 31, 2014 #1
    HI guys,this is my first programming experience , i have developed an matlab code for steady state heat conduction equation , on governing equation
    dt2 /dx2 + dt2/dy2 = -Q(x,y)

    i have solved this equation with finite difference method, As far as i know if we increase the mesh size it leads to decrease in the distance between the nodes and leads us to the true solution.i guess the final solution will be constant.

    but my code is not converging or moving towerds true solution can anybody help me out on this. Here the solution is just varying with the mesh size.
    what can be the reason behind this.

    this is my code.

    %%%%%% STRICTLY ONLY FOR ODD NUMBERS %%%%%%%%%

    clear all
    close all
    n = 101;% n grids and has n interior points per dimension
    m = n+2;
    x = linspace(0,1,m); % grid points x including boundaries
    y = linspace(0,1,m); % grid points y including boundaries
    dx = x(2)-x(1);
    h = dx;
    dy = dx;
    k=1; %thermal conductivity
    Q=8000; % heat source
    Iint = 2:n+1;
    Jint = 2:n+1;
    %Right hand side matrix
    if (mod(n,2)==0)

    %Main Matrix
    Told = zeros(m,m);

    % boundary values
    Told(1:m,1) = 0;
    Told(1:m,m) = 0;
    Told(1,1:m) = 0;
    Told(m,1:m) = 0;
    Tsoln = Told;
    B(:,1) = B(:,1) + Told(Iint,1)*k/h^2; % boundary valuues are added from the left hand side matrix including k and h
    B(:,n) = B(:,n) + Told(Iint,m)*k/h^2;
    B(1,:) = B(1,:) + Told(1,Jint)*k/h^2;
    B(n,:) = B(n,:) + Told(m,Jint)*k/h^2;

    % Matrix formations
    F = reshape(B,n*n,1);
    I = speye(n);
    e = ones(n,1);
    T = spdiags([e -4*e e],[-1 0 1],n,n);
    % v = full(T);
    S = spdiags([e e],[-1 1],n,n);
    % w = full(S);
    % v = kron(I,T);
    % w = kron(S,I);
    % z = (kron(I,T) + kron(S,I));
    % nny = full(z);
    A = (kron(I,T) + kron(S,I))* k / h^2; % A is a tridiagonal spurse matix multiplied and dvided by h^2
    % z = full(A);
    % Tvec = abs(A\F); Solving AX = B ;

    Tvec = abs(mldivide(A,F));
    Tsoln(Iint,Jint) = reshape(Tvec,n,n);
    sol = max(max(Tvec)) ;

    % Ploting
    figure, contour(x,y,Tsoln),
    title('Temperature (Steady State)'),xlabel('x'),ylabel('y'),colorbar
    figure,pcolor(x,y,Tsoln), shading interp
    title('Temperature (Steady State)'),xlabel('x'),ylabel('y'),colorbar
  2. jcsd
  3. Aug 31, 2014 #2

    Dr Transport

    User Avatar
    Science Advisor
    Gold Member

    What are your boundary conditions, more fully explain the problem and we can be more helpful
  4. Aug 31, 2014 #3
    thanks for your reply
    i am solving in a two dimensional staedy state heat conduction equtation with a heat source at the center in square plate with dirichlet boundary conditions, i mean to say values on boundries of four sides are same and constant.
    dt2 /dx2 + dt2/dy2 = -Q(x,y)
    i solved it by finite difference method
    B(i,j) is a right hand side matrix with heat source at the center.
    i am solving it by direct method, creating right hand side matrix and main spurse matrix. later i am dividing them for solution.
    My discritised equation is
    T(i,j) = ((T(i+1,j) +T(i-1,j) +T(i,j+1)+T(i,j-1) )*k)+(B(i,j)*dx^2)))/(4*k)
  5. Aug 31, 2014 #4
    Why are you multiplying T by k/h^2 ?

    If this is a 2D problem you are trying to solve, the cross sectional area shared by two adjacent cells is the length of the cell side times the depth...in 2D problems, one usually assumes 1 unit of length deep...

    ...and, so, I an thinking, you should be multiplying by k/(h*1)

    ...otherwise, you are changing the depth of your model and that is why you are getting different temperature every time.
  6. Aug 31, 2014 #5
    Or...if I am way off, never mind my point...I don't think I understand how you are setting up the problem, where your cells are, where your nodes are and whether the boundary nodes are numbered and part of your program (they shouldn't be). I certainly don't understand your choice of m (n = n+2), or why you presume your stiffness matrix A to be tridiagonal (it is only triadiagonal for a 1D mesh where every node had only 2 neighbors).
  7. Aug 31, 2014 #6
    Hey gsal, I appriciate your reply

    when i discritised this equation "d2t/dx2 + d2t/dx2 = -Q " i wrote the linear eauation
    ((T(i+1,j) +T(i-1,j) +T(i,j+1)+T(i,j-1) -4*T(i,j))*k)/ h2 = Q,
    above in my program h is the distance between the nodes
    this has lead me to the tridiagonal matrix or a sparse matrix ,framing a system of linear equations ax = b

    after seeing your message i tried it by using h^1, but it didnt work.i think we can not change that, as far as your second doubt, i am trying to solve n nodes of mtrix dimension of n+2,to set the boundary conditions to it.

    if you still feel that there is some thing wrongwith my code. feel free to edit my program and do it in your way. i want to see that too.
  8. Aug 31, 2014 #7
    Forgive me if I choose not to write the program for you. First, I don't use matlab; although it does not look like it would be much different from Python.

    Anyway, setting up this problem with a few loops is not that difficult, but you really need to know exactly what you need to do...if I were you (I already was and did it), I would really do my homework and deduce the equations by hand, first.

    They may sound like too many but I think you need at least a 3x3 grid (9 equations) to realized the pattern of the equations; otherwise, you should go for 4x4.

    If you welcome my advice, forget about your dt^2/dx^2...equation for now.

    I would take the 3x3 set of nodes, number them from 1 to 9 and setup Kirkchhoff's Equations of heat flowing away from each node toward its neighbors and BC "nodes" (these are not numbered, they are actual numeric values in the equations and the conductance connection to this temperature should half the one between actual nodes...you'll see).

    ...THEN, you will see a nice pattern on how each diagonal and off-diagonal entry looks like...I am thinking you really 4x4 cells, but you don't really need to set up all 16 equations, just the unique ones...10 ;-) ...4 corners, 4 sides, 1 diagonal (not corner), 1 off-diagonal (not side).

    To me, the temperature difference (dT) between two nodes equal the heat flow (Q) times the thermal resistance (R) in between, and so: dT = QR, where R=ρL/A = L/(κA) = dx/(κ*dy*1) = 1/κ, if dx=dy...so, I don't see how h^2 came to be in your equations.

    So, do Kirkchhoff's equations and organize them in the form Ax=b

    • b is the known 1D vector of heat injections into the system (from center and boundary conditions)
    • x is the 1D vector of temperatures to solve for
    • A is the 2D stiffness (conductance) matrix
  9. Aug 31, 2014 #8
    i appriciate for your conversation,
    i have done my homework and solved it numerically as you said taking 3*3 matrix and i happen to program my code in the same way i have calculated, now i get to know why you are facing the the problem,if you just run my program line to line you can check in matlab what am i doing, i have just done the same thing what you have written up there.
    i have checked my program many times by calculatng with hands numerically by taking 3*3 and they turn out to be same.so i had no clue why this program is not converging if i inrcrease the mesh size.
    let me clarify you again h^2 is distance between the nodes and it came in my equation while discritizing the main equation through Finite difference method.
    i really didnt understand what are you talking about (dT).
  10. Aug 31, 2014 #9
    I am a steady-state, finite-difference kind of guy and, so, I stay away from differential equations.

    In my equation, dT is nothing else than (T1 - T2), for example:
    T1 - T2 = Q12R12
    and so
    Qij = (Ti - Tj) / Rij
    and so, the summation of the various (Ti - Tj) / Rij for all j connected to the same i is equal to injection Qi at node i.

    Here, the denominator is thermal resistance connecting node i to node j, not just simply the distance between them, nor squared; like I illustrated before when I expanded/simplified dT=QR, in previous post, dx and dy cancel out leaving k/1, yielding the correct units on both sides of the equation...for a 2D problem, 1 unit deep, the units on both sides of Q = (T1-T2)/R are watts per meter squared (W/m2)...as it should be.

    If you do a dimensional analysis in your linearized equation, the left hand side has temperature, times thermal conductivity divided by linear squared: K*W/Km/m2 = W/m3...that's a volumetric heat density...how can that be the heat flow across the wall shared between two cells? The units need to be W/m2.

    So, I think I going back to what I said before...because the units of the left hand side are W/m3 and you have h2 in there, the meaning of Q on the other side is not an absolute, fix amount of watts every time, instead it is a density...so, every time you change the size of the one element with heat generation, you change the amount of total heat being injected into the system and when you solve whatever problem you are solving, you of course get a different answer...that's all I can speculate.
  11. Aug 31, 2014 #10
    Actually, let me take the part about W/m2 half-way back...sure, steady-state heat rate through a shared wall can only be in W/m2, but the heat from one node to the other (after multiplying by the area of the wall) should truly just be in plain W (watts), that's it.
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Similar Threads - Convergence code Date
Region of convergency Jul 11, 2017
Matlab ode solvers - adding a separate convergence criteria Dec 10, 2015
Mesh Convergence Issue In Ansys Feb 9, 2013
Convergence of an ODE May 17, 2012