I am trying to develop a code that will allow me to solve the following problem. Lets say that I have a square block, inside that square block I have a smaller square box thats at a certain temperature. The outside of the outer square is also at a known temperature. What I am trying to find is the temperature distribution in the larger square due to heat transfer from the small square.

The following code allows me to find the temperature distribution in the larger square.

    diffFormula = (T[i + 1, k] - 2 T[i, k] + T[i - 1, k])/x^2 +(T[i, k + 1] - 2 T[i, k] + T[i, k - 1])/y^2 == 0

The boundary conditions are

NP=5; (*number of grid points*) (*so a 6x6 matrix*)
Table[T[0, k] = 325, {k, 0, NP}];
Table[T[i, NP] = 325, {i, 0, NP}];
Table[T[i, 0] = 325, {i, 0, NP}];
Table[T[NP, k] = 325, {k, 0, NP}];

Now I setup the list of equations

eq = Flatten[Table[diffFormula, {i, 1, NP - 1}, {k, 1, NP - 1}]];

Followed by a list of variables

var = Flatten[Table[T[i, k], {i, 1, NP-1}, {k, 1, NP-1}]];

Followed by the solution

solution = N[Flatten[Solve[eq, var]]]
s = N[Table[T[i, k], {k, 0, NP}, {i, 0, NP}]] /. solution;

My problem is, lets say that I would like positions T[2,2],T[2,3],T[3,2],T[3,3] to represent the smaller square that is generating heat at 400. Simply setting these coordinates to 400 does not work as it passes this value to the variables list which causes interference with Solve. I've tried many ways of making the adjustment, however I can not find the right way to go about this.

I would really appreciate any help, Thank you.


1 Answer 1


Just remove the undesirable equations and variables:

diffFormula = (T[i + 1, k] - 2 T[i, k] + T[i - 1, k])/x^2 + 
              (T[i, k + 1] - 2 T[i, k] + T[i, k - 1])/y^2 == 0;

NP = 5;(*number of grid points*)
x = y = 0.05;
Table[T[0, k] = 325., {k, 0, NP}];
Table[T[i, NP] = 325., {i, 0, NP}];
Table[T[i, 0] = 325., {i, 0, NP}];
Table[T[NP, k] = 325., {k, 0, NP}];

inner = {T[2, 2], T[2, 3], T[3, 2], T[3, 3]};

eq = Flatten[Table[diffFormula, {i, 1, NP - 1}, {k, 1, NP - 1}]];
selectedeq = Select[eq, Count[#, Alternatives @@ inner, Infinity] < 2 &];

var = Flatten[Table[T[i, k], {i, 1, NP - 1}, {k, 1, NP - 1}]];
selectedvar = Complement[var, inner];

(# = 400) & /@ inner;

solution = Flatten@Solve[selectedeq, selectedvar];
s = Table[T[i, k], {k, 0, NP}, {i, 0, NP}] /. solution;


enter image description here

BTW, though the region you're dealing with is luckily simple, handling irregular region with FDM can be really frustrating. This post is an example. If FDM isn't necessary for you, have a look at the FEM capabilities of NDSolve new added in v10.

  • $\begingroup$ Thank you very much for the help. Ill have to go through the code that you added because at first glance I dont exactly understand what the adjustments do. I completely agree with your initial comment that I needed to remove the certain variables and associated equations however I did not know how to go about it. Im looking forward to going through your modifications. I was also planning to modify the square into a circle so thank you again for the link. Ill try my hand it. $\endgroup$
    – Valentin
    Jan 6, 2016 at 21:04
  • $\begingroup$ @Valentin Feel free to continue to ask in the comment if you still have difficulty in understanding the solution :) $\endgroup$
    – xzczd
    Jan 7, 2016 at 2:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.