Solving for Potential due to charge using Poisson Equation using FDM in Python (1D)

Question

Im trying to get electric potential distribution due to charge (or charges) in one plane (1D). To get it I'm solving Poisson Equation using FDM in Python. In my example, the Poisson equation looks like this: A times C = Rho. Where A =

[[-2  1  0 ...  0  0  0]
 [ 1 -2  1 ...  0  0  0]
 [ 0  1 -2 ...  0  0  0]
 ...
 [ 0  0  0 ... -2  1  0]
 [ 0  0  0 ...  1 -2  1]
 [ 0  0  0 ...  0  1 -2]]

Rho =

 [0.0e+00  0.0e+00  0.0e+00  0.0e+00 -1.6e-19  0.0e+00  1.6e-19  0.0e+00
  0.0e+00  0.0e+00]

When using C = np.linalg.solve(A, Rho) my output plotted in graph looks like this: graph I tried lowering dX and increasing X range to receive higher number of results. The only thing that helps me is changing the A matrix type to integer, but it changes the matrix values (-199 instead of 200, and 99 instead of 100, because 1 and -2 multiplied by 1/dX times dX)

A = 
[[-199   99    0 ...    0    0    0]
 [  99 -199   99 ...    0    0    0]
 [   0   99 -199 ...    0    0    0]
 ...
 [   0    0    0 ... -199   99    0]
 [   0    0    0 ...   99 -199   99]
 [   0    0    0 ...    0   99 -199]]

graph2 (Here i changed dX so the values are higher because they are multiplied by 1/dX**2) Can anyone give me direction where should i look for fault? I think its Rho Matrix problem, but i don't know what's wrong. Cheers.

Im uploading whole code here: (I'm beginner in Programming so it may look bad)

import numpy as np
import matplotlib.pyplot as plt

X = 30
dX = 0.1
N = int(X / dX)

x = np.arange(0, N)
x = x * dX

Rho = np.zeros(N)
Rho[int(10 / dX)] = -(1.6e-19) / dX ** 2
Rho[int(20 / dX)] = -(-1.6e-19) / dX ** 2

A = np.zeros((N, N), int)
np.fill_diagonal(A, -2 / dX ** 2)
A.ravel()[1::(2 * N + 2)] = 1 / dX ** 2
A.ravel()[2 + N::(2 * N + 2)] = 1 / dX ** 2
A.ravel()[N::(2 * N + 2)] = 1 / dX ** 2
A.ravel()[2 * N + 1::(2 * N + 2)] = 1 / dX ** 2
print(A)
print(Rho)
C = np.linalg.solve(A, Rho)

fig, ax = plt.subplots()
ax.plot(x, C, linewidth=1.0)
plt.show()

I've tried changing the dX, X and N values. I was expecting that program needs more points to calculate and it will resolve my problem. The only thing that "fixes" my problem is changing A matrix type to integer.

Answer 1

One comment asked "why do you think the 2nd graph is wrong"; My question would rather be "why do you think the second graph is correct". I mean, you saw the problem yourself : rounding error makes the matrix wrong. 100 -200 100, so 1 -2 1, to a 100 factor, is just a classical second derivative scheme. So, with a rho that is a pair of dirac (0 everywhere, but for 2 places), what you are solving is the equation "second derivative is 0 everywhere, but for two points".

And second derivative 0 everywhere means "lines everywhere". But at the non-0 places. So, the 1st graph.

Sure, the equation, with floats, is playing a bit too closely with floating point precision. It would be safer to solve AX = Rho* 10**19, and then C=X*10**-19. But that is not really the problem here. You get straight lines, because that is what you are supposed to get.

But with the rounding error, with ints, you change the equation

[99, -199, 99] means, in terms of derivation scheme -f(x) + 99*(f(x-h)-2*f(x)+f(x+h)) = f(x) + 99h² f''(x)

So what you are solving is ODE -y + αy'' = 0 everywhere, but in 2 places.

with α>0, that is y''=βy, (with β=1/α positive). So basic solution for that is y=exp(√βt) (times some line, to account for limit condition). Which is what you see in your second graph.

But that exponential is there because of a rounding error, that changes a y''=0 into a y''=βy in your derivative scheme.

That answers not to your (imho wrong) question "why I can't get the false answer with floats", but to the question "why to I get wrong answer with int".

But, well, in short, you are indeed dangerously close to floating points limit with that rho and equation. But that is not what explain differences here. What explain differences is that you are solving two different ODE.