# Why is my code not giving me a graph? (finite difference method) [closed]

For a finite difference method code, I after altered it a bit for my needs. I am no longer getting errors, however I am not getting my graph to show.

Mx = 100;
My = 100;
Lx = Ly = 1.;
dx = Lx/Mx;
dy = Ly/My;
V[0, j_] = 0;
V[Mx, j_] = 1;
V[i_, 0] = 0;
V[i_, My] = 0;

p[i_, j_] = 0;

var = Flatten[Table[V[i, j], {i, 1, Mx - 1}, {j, 1, My - 1}]];

eqs = Flatten[
Table[(V[i + 1, j] - 2 V[i, j] +
V[i - 1, j])/(dx)^2 + (V[i, j + 1] - 2 V[i, j] +
V[i, j - 1])/(dy)^2 ,
p[i, j], {i, 1, Mx - 1}, {k, 1, My - 1}]];

sol = Solve[eqs, var][[1]];

Vsol = Interpolation[
Flatten[Table[{i dx, j  dy, V[i, j]}, {i, 0, Mx, Mx/10}, {j, 0,
My, My/10}] /. sol, 1]];

ContourPlot[Vsol[x, y], {x, 0, 0.1}, {y, 0, 0.05},
FrameLabel -> {"x(m)", "y(m)"},
PlotLabel -> "Finite Difference Method Solution when p_v = 0"]


I am unsure for what I need to do here. Only thing I can think of is to maybe have all arguments be in terms of [i_ , j_] like in the declared variables part.

I know for this case (when p[i_, j_] = 0), I am to get a very similar solution as in this thread: Not sure how to set up the Laplacian/Poisson Equation

• What is the definition of V[i_,j_] ? As you note you need to define this. – Hugh Mar 24 '19 at 19:34

The problem was that Solve returned {} due to several typos. This should work:

Mx = 100;
My = 100;
Lx = Ly = 1.;
dx = Lx/Mx;
dy = Ly/My;
ClearAll[V];
V[0, j_] = 0;
V[Mx, j_] = 1;
V[i_, 0] = 0;
V[i_, My] = 0;
p[i_, j_] = 0;
var = Flatten[Table[V[i, j], {i, 1, Mx - 1}, {j, 1, My - 1}]];
eqs = Flatten[
Table[(V[i + 1, j] - 2 V[i, j] +
V[i - 1, j])/(dx)^2 + (V[i, j + 1] - 2 V[i, j] +
V[i, j - 1])/(dy)^2 == p[i, j], {i, 1, Mx - 1}, {j, 1,
My - 1}]];

sol = Solve[eqs, var][[1]];

Vsol = Interpolation[
Flatten[Table[{i dx, j dy, V[i, j]}, {i, 0, Mx}, {j, 0, My}] /.
sol, 1]];


This is however very inefficient. I would advise you to learn to use SparseArray in order to formulate this equations as a linear equation and to solve it with LinearSolve. That will be orders of magnitude faster.

• Haha, I definitely don't doubt that it was inefficient, this way made the most sense to me. However, I can't necessarily say I agree with the graph shown. I am doing the same problem as in the link to the thread I mentioned earlier. That solution was spot on in what I needed in that the top of the boundary was the value of 1, while in this case it seems to show that the corner has the value of 1 instead of the boundary. EDIT: I also should note that when I have the value of 1 set for V[Mx, j_], the values according to the legend are off by a very large factor. – LtGenSpartan Mar 24 '19 at 20:03
• The order in which you enter the conditions V[0, j_] = 0; V[Mx, j_] = 1; V[i_, 0] = 0; V[i_, My] = 0; is relevant; Mathematica tries to use the most special definition but in this case, there are several definitions that would have the same precedence. In this case, Mathematica uses the definition that was entered first. – Henrik Schumacher Mar 24 '19 at 20:56
• Yes, but when I have moved that value to the other boundaries, I'm still not getting a sensible graph. I have changed values of Mx, My, and p[i_, j_], but the graph is unchanging. However, increasing the values of Mx and My are giving me better values for the points, despite not changing how it looks. – LtGenSpartan Mar 24 '19 at 21:07
• That's another problem not related to the original question. (The sentence starting with "Yes, but..." is a good indicator for that.) As I said, this method is neither efficient nor robust, so I won't invest further time on fixing it. My suggestion to learn about SparseArray still holds. – Henrik Schumacher Mar 28 '19 at 9:40