# Laplace equation

I tried to run this code, but it doesn't work for NDSolve

Nx = 4; Ny = 4;

A = SparseArray[{{i_,
i_} -> -4., {i_, j_} /; ((i - j ) == 1 && Mod[j, Nx] != 0) ->
1., {i_, j_} /; (i - j == -1 && Mod[i, Nx] != 0) ->
1., {i_, j_} /; (i - j == Nx ) -> 1. , {i_, j_} /; j - i == Nx ->
1.}, {Nx Nx, Nx Nx}]

b = SparseArray[{{1, 1} -> 75, {4, 1} -> 50, {5, 1} -> 75, {8, 1} ->
50, {9, 1} -> 75, {12, 1} -> 50, {13, 1} -> 175, {14, 1} ->
100, {15, 1} -> 100, {16, 1} -> 150}];

tbl = (A.Table[T [k][t], {k, 1, Nx Nx}]) - b;

sol = NDSolve[
Flatten[{Table[1/(Nx + 1)^2 T[k]'[t] == tbl[[k]], {k, 1, Nx Nx}],
Table[T[k][0] == 0, {k, 1, Nx Nx}]}],
Table[T[k], {k, 1, Nx Nx}], {t, 0, 100}];

tbl2 = Join[{0},
Table[T[i][t] /. sol /. t -> 10, {i, 4 Nx + 1, 5 Nx}], {0}] //
Flatten

a = Table[{(i - 1)/(Nx + 1.0), tbl2[[i]]}, {i, 1, Nx + 2}]

plt1 = ListPlot[a, PlotStyle -> {Red, PointSize[0.02]}, Frame -> True,
ImageSize -> 400 {1, 1}, GridLines -> Automatic, AspectRatio -> 1]

• What are you trying to do? What do you expect to happen? – Gerli Apr 30 '15 at 14:15
• I tried to solve laplace equation by applying Finite difference method of solving PDEs by converting it to a system of ODEs @Gerli – Akram Ghanem Apr 30 '15 at 14:21
• In NDSolve replace tbl[[k]] with tbl[[k,1]]. As is, M returns T[i]==0 – LLlAMnYP Apr 30 '15 at 14:23
• Thanks @LLlAMnYP , NDSolve runs well, please see the rest of the code cuz it should give me numbers not symbols – Akram Ghanem Apr 30 '15 at 14:30
• The rest of the code gives T[17] thru T[20] which aren't defined or solved for anywhere. – LLlAMnYP Apr 30 '15 at 14:39

What follows is my latest edit. I have moved my previous attempts to solve this to the bottom of the answer.

EDIT 04.05.15

I have given this problem some more thought. It's peculiar, that Mathematica cannot solve it straight away, but it should in fact be quite easy. Using your notation, we are dealing with the following differential equation:

$$(N_x+1)^2 \hat A \cdot \vec T(t) - (N_x+1)^2 \vec b = \vec{T}{}'(t)$$

Straight away I would redefine $(N_x+1)^2 \hat A \rightarrow \hat A$ and $(N_x+1)^2 \vec b \rightarrow \vec b$. The solution of the associated homogeneous diffrential equation has a form of

$$\vec T(t) = \sum_k^{N_x^2} C_k e^{\lambda_k t} \vec V_k$$

where $\lambda_k$ and $\vec V_k$ are the eigenvalues and associated eigenvectors of the redefined matrix $\hat A$ and $C_k$ are arbitrary constants. The vector function $\vec T(t)$ must be augmented by a constant vector to account for the $-\vec b$ summand. I start from scratch with the following code:

Nx = 4;

A = (Nx +
1)^2 SparseArray[{{i_,
i_} -> -4., {i_, j_} /; ((i - j) == 1 && Mod[j, Nx] != 0) ->
1., {i_, j_} /; (i - j == -1 && Mod[i, Nx] != 0) ->
1., {i_, j_} /; (i - j == Nx) -> 1., {i_, j_} /; j - i == Nx ->
1.}, {Nx Nx, Nx Nx}]

b = Flatten@((Nx +
1)^2 SparseArray[{{1, 1} -> 75, {4, 1} -> 50, {5, 1} ->
75, {8, 1} -> 50, {9, 1} -> 75, {12, 1} -> 50, {13, 1} ->
175, {14, 1} -> 100, {15, 1} -> 100, {16, 1} -> 150}])


The constant vector augmenting the solution of the homogeneous equation is found as follows:

p = LinearSolve[A,b]
(* {-40.9091, -29.4508, -26.6098, -30.6818, -59.1856, -50.2841,
-46.3068, -46.1174, -70.5492, -66.1932, -62.2159, -57.4811,
-81.8182, -81.7235, -78.8826, -71.5909} *)


I'll call this vector $\vec p$. To satisfy the initial conditions we must find the constants $C_k$. At $t=0$ the function is simply

$$\vec T(0) = \vec p + \sum_k^{N_x^2} C_k \vec V_k = 0$$

Let's define the matrix $\hat V$ as a matrix of the columns $(\vec V_1,\vec V_2,...,\vec V_{N_x^2})$ and $\vec C = (C_1,...C_{N_x^2})$. The above equation then simplifies to $\vec T(0) = \vec p + \hat V \cdot \vec C = 0$.

The matrix $\hat V$ is easily obtained:

{eigenvals, eigenvecs} = Eigensystem[A];
eigenvecs = Transpose@eigenvecs;


Thankfully

Det[eigenvecs]
(* -1. *)


which means, that the eigenvectors are linearly independent. Had that not been the case, everything would have been a whole lot more complicated. So, back to the matter. To find $\vec C$ we use the following command:

coefs = LinearSolve[eigenvecs,-p] // Chop
(* {0, 0.272672, 2.81034, 11.8769, 0, 0, 9.46032, -15.0092, 0.0503585,
9.14661, -3.92127, -71.2871, 0, -21.2122, 65.0471, -213.123} *)


Now, hopefully, we can construct a correct solution.

T[x_] := p + eigenvecs.(coefs (Exp[# x] & /@ eigenvals))


and test it:

T[0] // Chop
(* {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} *)
D[T[x], x] + b - A.T[x] // Simplify // Chop
(* {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} *)


Great, the solution satisfies both equations. Let's plot it.

Plot[Evaluate@T[x],{x,0,1}]


Unsurprisingly, this is an exponentially decaying solution (all the eigenvalues are negative) which very quickly - as fast as Exp[Max[eigenvals] x] settles off on a constant equal to p.

This is as close to a numeric solution as I can get. It was done with machine-precision arithmetic, not with exact numbers, but the analytical approach is valid in any case.

I haven't the faintest clue, why in the hell Mathematica cannot solve this by itself, also not sure why NDSolve insists on diverging so quickly.

One more attempt with DSolve

Like I've shown before, machine-precision arithmetic solves everything just fine. When using exact arithmetic which is possible if A is redefined as

A = Map[Rationalize, A, {2}];


mathematica fails to do

coefs = LinearSolve[eigenvecs, -p]


If we stay with machine precision numbers

tt = Array[t[#][x]&, 16];
sol = First@DSolve[A.tt - b == D[tt, x] && (tt /. x -> 0) == ConstantArray[0, 16], tt, x];
Plot[tt /. sol, {x, 0, 1}]


works fine and gives a solution that matches the solution from the previous approach to about 10^-10.

Finally, NDSolve started giving me meaningful results too.

tt = Array[t, 16];
sol = First@NDSolve[A.Through[tt[x]] - b == D[Through[tt[x]], x] &&
(Through[tt[0]]) == ConstantArray[0, 16],
tt, {x, 0, 10}]
Plot[Evaluate@(Through[tt[x]] /. sol), {x, 0, 10}]


I guess, the problem vanished on its own.

One last thought

Although Mathematica fails to solve this system exactly, we can deduce quite a lot about the exact solution.

A = Map[Rationalize, A, {2}];
{eigenvals, eigenvecs} = Eigensystem[A];
eigenvecs = Transpose@eigenvecs;
p = LinearSolve[A,b];


These steps are succesfully executed with exact numbers. The exact solution quickly settles at $\vec p$ equal to p. The coefficients for the exponents are exactly $\vec \lambda$ equal to eigenvals (all negative). It's just that $\vec C$ cannot be calculated exactly.

With some slight modification of the code:

Nx = 4; Ny = 4;

A = SparseArray[{{i_,
i_} -> -4, {i_, j_} /; ((i - j) == 1 && Mod[j, Nx] != 0) ->
1, {i_, j_} /; (i - j == -1 && Mod[i, Nx] != 0) ->
1, {i_, j_} /; (i - j == Nx) -> 1, {i_, j_} /; j - i == Nx ->
1}, {Nx Nx, Nx Nx}];

b = SparseArray[{{1, 1} -> 75, {4, 1} -> 50, {5, 1} -> 75, {8, 1} ->
50, {9, 1} -> 75, {12, 1} -> 50, {13, 1} -> 175, {14, 1} ->
100, {15, 1} -> 100, {16, 1} -> 150}];

tbl = (A.Table[T[k][t], {k, 1, Nx Nx}]) - b;

sol = DSolve[
Flatten[{Table[1/(Nx + 1)^2 T[k]'[t] == tbl[[k, 1]], {k, 1, Nx Nx}],
Table[T[k][0] == 0, {k, 1, Nx Nx}]}],
Table[T[k], {k, 1, Nx Nx}], {t, 0, 100}];

GraphicsRow[{Plot[Evaluate@Table[T[i][t], {i, 16}] /. sol, {t, 0, .23}],
Plot[Evaluate@Table[T[i][t], {i, 16}] /. sol, {t, 0, 1}]}]


we get an analytical solution. Ugh.

EDIT

I noticed that

tbl[[k]]


returned things like

{-75 - 4. T[1][t] + 1. T[2][t] + 1. T[5][t]}


which are lists (note the curly braces surrounding it). So I took the first and only part of that list with tbl[[k,1]] which returns the same, but not wrapped in curly braces.

In the definition of the SparseArray A I replaced all numbers like 4. and 1. with 4 and 1 to use exact arithmetic. And then as these are fairly simple first-order differential equations, I figured, Mathematica would solve them analytically just fine. Considering, how quickly they shoot off to infinity, numerical solving would be certainly way off-target. So I replaced NDSolve with DSolve, which returns exact analytic solutions.

EDIT2

On second thoughts, scratch that. Due to an oversight, I forgot to recalculate tbl and inexact arithmetic was still being used.

After I corrected this, Mathematica failed to compute a result for me within any reasonable timeframe. I tried calculating for Nx=2 and Nx=3 taking the first 4 or 9 elements of b, respectively, and that did get solved. As it seems, the solution comes up with expressions such as Exp[-54 t](1 + Exp[18 t])^2 (1 + 2 Exp[18 t]). Solutions like this settle off on a constant value at large t. Of course, tiny errors in the solution when solving numerically will rapidly propagate into very quick divergence to infinity.

I'm almost certain, that as Mathematica has no problems in finding the eigenvectors of the matrix A, with a bit of help it should solve the system for Nx=4, but it will involve a lot of manual labor. I'll let it run for some time and see if it eventually comes up with a solution by itself.

• I am so thankful @LLIAMnYP , I will appreciate if can tell me what is the mistake in my code :) – Akram Ghanem Apr 30 '15 at 14:40
• I was expecting a different graph @LLIAMnYP, also I am looking for a numerical solution – Akram Ghanem Apr 30 '15 at 14:47
• I just realized, that after redefining A with proper integers I forgot to recalculate tbl. Now with exact arithmetic the solution is taking quite a bit longer. – LLlAMnYP Apr 30 '15 at 14:52
• Just a note: you can use Evaluate@Flatten@Table inside of Plot to get the lines in different colors. – 2012rcampion Apr 30 '15 at 15:14
• Yep, it's an oversight on my part. I forgot, that sol` is a list within a list. – LLlAMnYP Apr 30 '15 at 15:25