# Poisson PDE in a rectangular domain

I found a solution for this problem, but this is in Scilab and I never use Scilab. Can anyone can help me to translate it in Mathematica? Here is the Link:http://imechanica.org/files/TorqueR.pdf

I solved this problem with NDSolve but now I want to implement it with FDM, like this.

Clear[a, b, n, m, Δx, Δy, G, θ,
sol, vars, ϕsol, eqns, tauyz, tauxz, ϕ, gammaxz, gammayz]
a = 0.04;
b = 0.04;
n = 5;
m = 5;
Δx = a/n;
Δy = b/m;
G = 13.125 10^6;
θ = 0.002;
ϕ[0, j_] = 0;
ϕ[n, j_] = 0;
ϕ[i_, 0] = 0;
ϕ[i_, m] = 0;
Ψ[i_, j_] = -2 G θ;
vars = Flatten[Table[ϕ[i, j], {i, 1, n - 1}, {j, 1, m - 1}]];
eqns =
Flatten[
Table[
(ϕ[i + 1, j] - 2 ϕ[i, j] + ϕ[i - 1, j])/Δx^2 +
(ϕ[i, j + 1] - 2 ϕ[i, j] + ϕ[i, j - 1])/Δy^2 == Ψ[i, j],
{i, 1, n - 1}, {j, 1, m - 1}]];
sol = Solve[eqns, vars][[1]];
ϕsol =
Interpolation[
Flatten[Table[{i Δx, j Δy, ϕ[i, j]}, {i, 0, n}, {j, 0, m}] /. sol, 1]];
tauyz = -D[ϕsol[x, y], x];
tauxz = D[ϕsol[x, y], y];
gammayz = tauyz/G;
gammaxz = tauxz/G;


but this only works for a rectangle, but I need to solve for a rectangle with hole in the center, like this:

• According to the link, the program is written with SciLab
– user9660
Commented Mar 31, 2015 at 16:52
• if you want Poisson solver using FDM on rectangular cross section, there is a demo that implements this with many other solvers as well. demonstrations.wolfram.com/… Commented Mar 31, 2015 at 17:05
• When the hole is also rectangular shape, then makes it easier, since the grid can be made to align there. with the rest of the region with no problem. So all what you have to do, is in the loop as you scan the region (updating the solution at each point), is simply check if you are inside the hole or not. Commented Mar 31, 2015 at 17:20
• @jens Ah, I see. But from reading that question it seems he has 10.0.0 which can be freely updated to 10.0.2 Commented Mar 31, 2015 at 17:52
• "I would like it for any geometry" - 如果这就是你的终级目标的话，那么，我在你的第一个问题里给出的代码本来就是适用于任何二维直角坐标系下的泊松方程的，你只需要把开头用于指定区域的部分稍微改改就行——‌​‌​你要是看不懂我的答案那你完全可以在下面追问。Translation: If this is your ultimate goal, then the code in my answer for your first question is completely suited for the task, you just need to modify the part defining the region i.e. rulei and ruleo. If you have any difficulty in understanding, feel free to continue to ask in the comment under my answer. Commented Apr 1, 2015 at 3:40

This is a very simple way to do it. I modified your code very little to add the condition, that is inside a hole, then set the solution to be the boundary condition at the edge of the hole inside, which I set to be zero.

  startRow = 4; endRow = 6; startCol = 4; endCol = 6;


Which tells the boundaries of the hole. And inside the Table added an If to check:

If[i >= startRow && i <= endRow && j >= startCol && j <= endCol,
ϕ[i, j] == 0,
.... same as before
],


With more mesh, the solution should become better.

Here is your code with the small modification

Clear[a, b, n, m, Δx, Δy, G, θ, sol, vars, ϕsol, eqns, tauyz, tauxz, ϕ, gammaxz, \
gammayz]
a = 0.04;
b = 0.04;
n = 100;
m = 100;
Δx = a/n;
Δy = b/m;
G = 13.125 10^6;
θ = 0.002;
ϕ[0, j_] = 0;
ϕ[n, j_] = 0;
ϕ[i_, 0] = 0;
ϕ[i_, m] = 0;
Ψ[i_, j_] = -2 G θ;
startRow = 40; endRow = 60; startCol = 40; endCol = 60;
vars = Flatten[Table[ϕ[i, j], {i, 1, n - 1}, {j, 1, m - 1}]];

eqns = Table[
If[i >= startRow && i <= endRow && j >= startCol && j <= endCol,
ϕ[i, j] == 0, (*assume B.C. in hole edge is zero *)
(ϕ[i + 1, j] - 2 ϕ[i, j] + ϕ[i - 1, j])/Δx^2 + (ϕ[i, j + 1] -
2 ϕ[i, j] + ϕ[i, j - 1])/Δy^2 == Ψ[i, j]
],
{i, 1, n - 1}, {j, 1, m - 1}
];

sol = Solve[Flatten@eqns, vars][[1]];

ϕsol = Interpolation[Flatten[Table[{i Δx, j Δy, ϕ[i, j]}, {i, 0, n}, {j, 0, m}] /. sol, 1]];
tauyz = -D[ϕsol[x, y], x];
tauxz = D[ϕsol[x, y], y];
gammayz = tauyz/G;
gammaxz = tauxz/G;