[PYTHON] Numerical analysis of electric field: Finite difference method (FDM) -Practice-
Preface
Continuation of the previous "Numerical analysis of electric field: Finite difference method (FDM) -Basics-".
This time, I will finally write not only mathematical formulas but also finite difference method with python. Even so, I'm a little embarrassed to write a program on the net because it seems that I'm not good at it. Especially in python, the cleanliness of the script directly leads to speeding up and optimization, so it seems that various things will be rushed in.
In that sense, C, C ++, Fortran, etc. are faster even if you turn a for loop etc. into the brain muscle, and recently I have come to like this one. ...... The story went awry. Let's start numerical analysis.
The program I actually made has github in the summary section, so if you want to run it once, please download it first and run it.
Problem setting
It is difficult to write a program that can solve all problems from the beginning, so let's set the problem first.
Poisson equation
The calculation area is defined as the following isosceles trapezoid, the upper and lower sides are Dirichlet boundary conditions, and the left and right are Neumann boundary conditions, and are set as follows.
** Figure 1 **: Computational area
Let's write it all together.
\begin{align}
\nabla \cdot \nabla u &= 0 \ \ \text{(In the area)}\\
u &= 1 \ \ \text{(Top side)} \\
u &= -1 \ \ \text{(Bottom side)} \\
\frac{\partial u}{\partial n} &= 0 \ \ \text{(Left and right sides)}
\end{align}
Let's write a program that solves this problem. However, it's boring to write a program that solves only this, so let's implement it while considering some generality.
The libraries used this time are only numpy and matplotlib.
Implementation
Before saying the program ..., let's first write a simple class diagram. Even if you say a class diagram, it is not shared by everyone, it is written for the purpose of organizing your own brain, so it is not written in a sloppy manner, it seems that there is such a class. I think (prevention line).
First, let's check the FDM flow.
--Step 1: Problem setting --Area setting / Source item setting --Step 2: Grid division --Step 3: Solve with FDM --Generate simultaneous equations and solve them --Step 4: Result
I honestly divided these four into large classes.
Figure 2: Class diagram
Let me explain briefly. The PoissonEquation class is the so-called main one. Then, Problem was set as a problem, and it was divided into two, the domain domain class and the source term Source class. Boundary conditions may be separated, but I decided to include them in the Domain class.
Next, in the grid division class, I put it in the Grid class. Grid is made up of classes (Node) that handle grid points. I purposely created the Node class because I was thinking of dealing with line segments (Edge) and areas (Cell) connecting each grid point, which will be used in the finite element method, which I plan to do next time. ..
Next, the Solver class also puts FDM under it, because it also puts other methods such as FEM under it. Although it is not specified by the line, should it be a child class? (Honestly, I don't know because I haven't made it yet)
Finally, the result was divided into potential (Potential) and electric field (Flux Density).
Alright, we've got a direction and let's write a program! However, there is a high possibility that it will change little by little while making it.
By the way, I'm the type who writes in detail from the end of the class, but is this common ...? I think it's better to read a book like "The royal way of writing a program". Well, this time I will proceed according to my own method, so I'm sorry if it is a strange method.
Step 1: Problem setting
First, let's set the problem programmatically. The first thing to think about at this time is how to save it as data.
Domain setting (Domain class)
For example, when considering that it corresponds only to the above-mentioned isosceles trapezoid problem, the shape can be uniquely determined as long as there is "bottom side length", "upper side length", and "height", which is sufficient as data. .. However, with such data, isosceles trapezoids, rectangles (upper side = lower side), and squares (upper side = lower side = height) can be realized, but all other shapes cannot be supported, and there is no generality. Not interesting. On the contrary, it is difficult to find too much generality, and it is almost impossible to express all vertices and curves with a finite number of data. You could do it by passing a function, but it's obvious that the program becomes quite complicated. In this way, seeking generality tends to complicate the program, and conversely, trying to simplify the program tends to impair generality (just my rule of thumb, I admit disagreement).
This time, let's keep the "coordinates of the vertices of the polygon" as an array. By doing this, it can handle all polygons, maintain a certain degree of generality, and it seems to be easy as data. It can't handle circles, but ... well, if you make it a regular icosagon, you can approximate it to some extent, and you should expand it when necessary. The source is as follows.
domain.py
import numpy as np
import matplotlib.pyplot as plt
class Domain(object):
def __init__(self, shape, bc):
self.shape=shape
self.bc = bc
class Polygon(Domain):
def __init__(self, vertexes, bc, priority=None):
super().__init__("Polygon", bc)
self.nVertexes=len(vertexes)
self.vertexes = np.array(vertexes) #If the argument vertexes is a list, np.Replace in array.
The reason why I created the Domain class and made it the Polygon class as a child class is that there is a possibility that I will create a Circle class so that it can handle circles in the future. On the contrary, it is also possible to create a rectangle class and process it at high speed. Anyway, I did this for extensibility. The author is not a professional, but an amateur, so I'm not saying that this is correct (prevention line). I think it's a good idea to write your own code.
As you can see from the above code, the data we have has shape, bc in the parent class and nVertexes, vertexes in the small class.
As you can see from the variable name and code, it is defined as follows.
--shape: The name of the shape of the area to be used. shape = "Polygon"
--bc: Boundary condition
--nVertexes: Number of vertices in Polygon
--vertexes: Coordinates of vertices. If you connect them in order from vertexes [0], you can make a polygon.
By the way, how should the boundary conditions be expressed? This time, I will express it in the following dictionary type.
bc = {"bc":[{"bctype":"Dirichlet"Or"Neumann" , "constant":constant}, ...], "priority":[...]}
Specify the constant on the right side of the boundary condition in " constant " and the Diricre / Neumann boundary condition in " bctype ".
" bc " is made as a list of " constant " and " bc type ". It means that bc ["bctype "] [0] is the boundary condition of the boundary connecting vertexes [0] and vertexes [1].
As you can see, it only supports constants ... this is just a hassle. I think that you can easily realize other than constants by passing a function type as an argument.
Which boundary condition does " priority " take precedence? I added it in the sense that. For example, which of the boundary conditions on the vertices of vertexes [1] is bc ["bctype"] [0] orbc ["bctype"] [1]? It is for specifying that.
The " priority " ruled that the priorities should be stored in a list in order. I don't explain in detail because I think that the rules that should be implemented differ depending on the person.
In addition, we added calcRange to calculate right angles including polygons and plot to plot for debugging. There is no need to write it, so I will omit it.
As a result, the following figure was obtained.
Figure 3: Plot of Domain class
Alright, the calculation area is now set.
Source term specification (Problem class)
Next, let's create the Source class. This is straightforward and can take an expression (with a two-dimensional array as an argument) as an argument. The source of the Source class itself is as follows.
Source.py
import numpy as np
class Source(object):
def __init__(self, source):
self.source=source
…… That, it fits in one line. I don't think this kind of classification is good ... Let's be careful not to do this! It's a bit off topic, but I think python doesn't need getters and setters because it can change and get class member variables directly from the outside. I have the impression that python is for easy use, so I think this is wonderful. I personally use python for the prototype and C ++ for the production.
By the way, if this is left as it is, for example, even if the source term is set to $ 0 $, function = lambda x: 0 must be given as an argument using the anonymous function lambda.
It's a bit annoying, so if the source term is a constant term, edit it so that you can pass a constant.
The code is modified as follows.
Source.py
import numpy as np
class Source(object):
def __init__(self, source):
print("-----Source-----")
if isinstance(source, (int, float)):
self.source = lambda x: source
else:
self.source = source
Now, when the argument source is int type or float type, it seems to work as a constant function.
Domain, Source call (Problem class)
Let's briefly write about the Problem class that summarizes the Domain and Source classes, which also serves as a summary.
Problem.py
from . import Domain
from . import Source
class Problem(object):
def __init__(self, problem):
self.setDomain(**problem["domain"])
self.setSource(problem["source"])
def setDomain(self, **domain):
if domain["shape"]=="Polygon":
self.domain=Domain.Polygon(domain["vertexes"], domain["bc"])
self.domain.calcRange()
def setSource(self, source):
self.source=Source.Source(source)
Let me briefly explain the arguments. An example of the argument is written on the next page.
setDomain(
shape = "Polygon"
vertexes = np.array([[-1,-1], [1,-1], [0.4,1],[-0.4,1]])
bc = {"bc": [{"bctype":"Dirichlet", "constant": 1},
{"bctype":"Neumann" , "constant": 0},
{"bctype":"Dirichlet", "constant":-1},
{"bctype":"Neumann" , "constant": 0}]
"priority":[0,2,1,3]}}
)
#source is a constant function f=When setting to 0
setSource(
source = 0
)
#source sin(x)*sin(y)When doing like
setSource(
source = lambda x: np.sin(x[0])*np.sin(x[1])
)
You can call these functions from the PoissonEquation class.
Step 2: Lattice generation
Next, the area is divided into grids. To be honest, this was the most difficult ... This section does not touch on detailed algorithms, but focuses on what kind of intention and what kind of specifications were used. Because [Computational Geometry](https://ja.wikipedia.org/wiki/%E8%A8%88%E7%AE%97%E5%B9%BE%E4%BD%95%E5%AD Is it a little off the subject because it is too close to% A6)? Because I thought. Of course, the method of grid division is also one of the deep points of FDM, but it is too deep, so I will leave it for a moment and make it a minimum story. If you are interested, study from various books and dissertations.
I made it in the following 4 steps.
- Arrangement of grid points (nodes) horizontal to the $ x $ axis and $ y $ axis
- Put a node on the boundary
- Erase nodes outside the area
- Definition of the next node
Manipulating grid points (Node class)
Step 1. Place the node
At first, I was thinking of separating it into a Node class and a NodeManagement class, but if you think about it carefully, python has a dictionary type. Store the parameters of each grid point (node) in a dictionary type, and make the Node class a class that manipulates the set.
The node parameters are as follows.
node={
"point":np.array(point),
"position":"nd",
"nextnode":[
{"no":None, "position":"nd"},
{"no":None, "position":"nd"},
...
]
}
Each explanation is summarized below
-"point": Node coordinates -"position": Node positional relationship ("nd": undefined "in": within the area "b (number n)": on the boundary n "v (number n)": on the vertex n) --"nextnode": List of information for the next node -"no": Index of the next node -"position": Positional relationship of adjacent nodes ("n": undefined "l": left "r": right "d": bottom "u": top)
Cartesian's __init__ (constructor) argument requires a Domain instance and the number of divisions (named variable name div).
For example, if you pass [10,10], let's set it so that it is divided into 10 equal parts vertically and horizontally. Then, it becomes like this.
Figure 4: Equally spaced grid
You may also want to try something that isn't evenly divided (eh? No?).
For example, [[-10, -7, -3, -1,0,1,3,7,10], [-10,-7,-3, -1,0,1,3,7,10 ]] If you pass it like , it will handle it well,
Figure 5: Non-equidistant grid
I tried to divide the grid like this. It may not be used in practice, but it is good for research. Now, this can be achieved as follows.
class Cartesian(Node):
def __init__(self, domain, div):
np.array(div)
if isinstance(div[0], (int)):
self.nDivX = div[0]
self.xs = np.linspace(domain.left, domain.right, self.nDivX)
else:
self.nDivX = len(div[0])
self.xs = np.array(div[0])
self.xs = np.sort(self.xs)
d = (domain.right - domain.left) / (self.xs[-1] - self.xs[0])
self.xs = d * self.xs
if isinstance(div[1], (int)):
self.nDivY = div[1]
self.ys = np.linspace(domain.down, domain.up , self.nDivY)
else:
self.nDivY = len(div[1])
self.ys = div[1]
self.ys = np.sort(self.ys)
d = (domain.up - domain.down) / (self.ys[-1] - self.ys[0])
self.ys = d * self.ys
self.nDiv = self.xs * self.ys
self.nodes = [{"point":np.array([self.xs[i],self.ys[j]]), "position":"nd", "nextnode":{"no":None, "position":"nd"} }
for j,i in itertools.product(range(self.nDivX), range(self.nDivY))]
It's a hassle, so I won't bother to explain it. Readers will be able to make their own in their own way. For the time being, just think that you made something that creates a node like the one shown above.
Step 2. Place a node on the boundary
Next, place the nodes on the boundary.
It doesn't matter if it overlaps with an existing node. Anyway, place it where it intersects the horizon of the $ x $, $ y $ axes of the existing node.
I made this as getBorderPoint in the Domain class. This time it's Polygon, but I thought the algorithm changed depending on the shape.
After placing the nodes on the boundary, if you delete the overlapping nodes, the result will be as shown in the following figure.
Figure 5: Boundary placement
... I made a mistake in setting the color. It may be hard to see because the red border and the magenta node overlap, but please be patient. You can see that the second node from the left and the second node from the top are very close to the boundary. I noticed this the last time I tested it, but it causes a lot of error, so it's reasonable to remove it. How close do you want to erase? Let's try to give something as an argument.
I wrote it lightly, but it was unexpectedly troublesome ... If you want to see the source, please see it on github.
Step 3. Erase the outer node
Next, the internal determination process is performed to delete the nodes outside the calculation area.
First of all, we have to judge the inside and outside of each, but various algorithms have already been proposed for this.
By the way, this time the node on the boundary is saved on the data as " position ":" b "+ (number n) in the above process, so edit only the " position ": nd. Just do it.
This is very helpful because the process of internal / external judgment at a node on the boundary is quite troublesome.
There are two famous methods. One is to draw a half-line in any direction from the point to be judged and count the number of crossings with the boundary (Crossing Number Algorithm). If the number of intersections is even, it is outside, and if it is odd, it is inside.
Figure 6: Crossing Number Algorithm
This method is fast because it is only four arithmetic operations, but there are some things to be careful of. As you can see from the example of ③ in the above figure, what to do when the vertex touches the line segment? You have to think about that. Actually, it is sufficient to impose the condition "Which direction is the vector of the line segment of the collision boundary?", But it seems to be troublesome.
So, let's adopt another method called Winding Number Algorithm this time. From the point of view of judgment, if the sum of the angles between each vertex is $ 2π $, it is inside, and if it is $ 0 $, it is outside.
Figure 7: Winding Number Algorithm
Since this method uses inverse trigonometric functions, it takes some time, but it is not necessary to consider the above conditions.
Now, I'm writing a program, but it's tricky to run a for loop on each node.
Since we are using numpy, let's manage to reduce the number of lines and reduce the calculation cost.
The program was as follows.
class Polygon(Domain):
#abridgement
def deleteOutDomain(self, node, thetaerr=0.1):
pts = np.array([node[i]["point"] for i in range(len(node))]) #["point"]Replaced in numpy
pos = np.array([node[i]["position"] for i in range(len(node))]) #["position"]Replaced in numpy
thetas = 0
for i in range(len(self.vertexes)):
p1 = self.vertexes[i] #End point of line segment 1 p1
p2 = self.vertexes[(i + 1) % self.nVertexes] #End point of line segment 2 p2
v1 = (pts - p1) #Vector from the point you want to judge to p1
v2 = (pts - p2) #Vector from the point you want to judge to p2
dot = v1[:,0] * v2[:,0] + v1[:,1] * v2[:,1] #inner product(dot,I didn't aim for inner ...)
cross = v1[:,0] * v2[:,1] - v1[:,1] * v2[:,0] #Cross product(Same as above)
vn1 = np.linalg.norm(v1, axis=1) #Get v1 distance
vn2 = np.linalg.norm(v2, axis=1) #v2 distance acquisition
theta = np.arccos(np.clip(dot / (vn1 * vn2), -1, 1)) #Clip because it can be more than 1 due to numerical error
thetas += np.sign(cross)*np.array(theta) #Add the angle every time for(cross is a plus / minus judgment)
inx = np.where( ((pos == "nd") & ~(thetas < thetaerr)))[0] # ["position"]Is nd and thetas are not almost 0, but get the index
#abridgement
If done well, it may be possible to eliminate even the vertexes loop, but I think that if you do this much, it will forgive you.
The result is as follows.
Figure 8: The outside is erased
4. Let's define what the next node is
Finally, let's define node [next node]. Defining this simplifies matrix generation for the next step.
I was worried about various things, but decided to judge the top, bottom, left, and right by the coordinates.
Floating point == (actually np.is close) comes out, so it may not be really good, but I couldn't think of anything else.
Maybe I should have had an integer index as data first, but please forgive me.
And before that, let's sort the node list by coordinates. I think it's easier.
You don't have to bother to raise it, so please refer to github.
(Actually, when I make a program, I say "Don't use for! "And it's embarrassing because it's duplicated ...)
So far, let's output the result of the node list.
Even if the number of data is too large, it cannot be seen comprehensively, so let's put it out around div = [4,4].
node0 {'point': array([-1., -1.]), 'position': 'c0', 'nextnode': [{'no': 1, 'position': 'r'}]}
node1 {'point': array([-0.333, -1. ]), 'position': 'b0', 'nextnode': [{'no': 5, 'position': 'u'}, {'no': 0, 'position': 'l'}, {'no': 2, 'position': 'r'}]}
node2 {'point': array([ 0.333, -1. ]), 'position': 'b0', 'nextnode': [{'no': 6, 'position': 'u'}, {'no': 1, 'position': 'l'}, {'no': 3, 'position': 'r'}]}
node3 {'point': array([ 1., -1.]), 'position': 'c1', 'nextnode': [{'no': 2, 'position': 'l'}]}
node4 {'point': array([-0.8 , -0.333]), 'position': 'b3', 'nextnode': [{'no': 5, 'position': 'r'}]}
node5 {'point': array([-0.333, -0.333]), 'position': 'in', 'nextnode': [{'no': 1, 'position': 'd'}, {'no': 9, 'position': 'u'}, {'no': 4, 'position': 'l'}, {'no': 6, 'position': 'r'}]}
node6 {'point': array([ 0.333, -0.333]), 'position': 'in', 'nextnode': [{'no': 2, 'position': 'd'}, {'no': 10, 'position': 'u'}, {'no': 5, 'position': 'l'}, {'no': 7, 'position': 'r'}]}
node7 {'point': array([ 0.8 , -0.333]), 'position': 'b1', 'nextnode': [{'no': 6, 'position': 'l'}]}
node8 {'point': array([-0.6 , 0.333]), 'position': 'b3', 'nextnode': [{'no': 9, 'position': 'r'}]}
node9 {'point': array([-0.333, 0.333]), 'position': 'in', 'nextnode': [{'no': 5, 'position': 'd'}, {'no': 13, 'position': 'u'}, {'no': 8, 'position': 'l'}, {'no': 10, 'position': 'r'}]}
node10 {'point': array([0.333, 0.333]), 'position': 'in', 'nextnode': [{'no': 6, 'position': 'd'}, {'no': 14, 'position': 'u'}, {'no': 9, 'position': 'l'}, {'no': 11, 'position': 'r'}]}
node11 {'point': array([0.6 , 0.333]), 'position': 'b1', 'nextnode': [{'no': 10, 'position': 'l'}]}
node12 {'point': array([-0.4, 1. ]), 'position': 'c3', 'nextnode': [{'no': 13, 'position': 'r'}]}
node13 {'point': array([-0.333, 1. ]), 'position': 'b2', 'nextnode': [{'no': 9, 'position': 'd'}, {'no': 12, 'position': 'l'}, {'no': 14, 'position': 'r'}]}
node14 {'point': array([0.333, 1. ]), 'position': 'b2', 'nextnode': [{'no': 10, 'position': 'd'}, {'no': 13, 'position': 'l'}, {'no': 15, 'position': 'r'}]}
Yeah, it feels good ... oh! I forgot.
I wanted to represent the nodes on the boundary as f (forward) and b (back). Since it is diagonal on the boundary, it cannot be expressed only by up, down, left, and right.
So I checked the program again.
node0 {'point': array([-1., -1.]), 'position': 'c0', 'nextnode': [{'no': 1, 'position': 'f'}, {'no': 4, 'position': 'b'}, {'no': 1, 'position': 'r'}]}
node1 {'point': array([-0.333, -1. ]), 'position': 'b0', 'nextnode': [{'no': 0, 'position': 'b'}, {'no': 2, 'position': 'f'}, {'no': 0, 'position': 'l'}, {'no': 2, 'position': 'r'}]}
node2 {'point': array([ 0.333, -1. ]), 'position': 'b0', 'nextnode': [{'no': 1, 'position': 'b'}, {'no': 3, 'position': 'f'}, {'no': 1, 'position': 'l'}, {'no': 3, 'position': 'r'}]}
node3 {'point': array([ 1., -1.]), 'position': 'c1', 'nextnode': [{'no': 2, 'position': 'b'}, {'no': 7, 'position': 'f'}, {'no': 2, 'position': 'l'}]}
node4 {'point': array([-0.8 , -0.333]), 'position': 'b3', 'nextnode': [{'no': 0, 'position': 'f'}, {'no': 8, 'position': 'b'}, {'no': 5, 'position': 'r'}]}
node5 {'point': array([-0.333, -0.333]), 'position': 'in', 'nextnode': [{'no': 4, 'position': 'l'}, {'no': 6, 'position': 'r'}]}
node6 {'point': array([ 0.333, -0.333]), 'position': 'in', 'nextnode': [{'no': 5, 'position': 'l'}, {'no': 7, 'position': 'r'}]}
node7 {'point': array([ 0.8 , -0.333]), 'position': 'b1', 'nextnode': [{'no': 3, 'position': 'b'}, {'no': 11, 'position': 'f'}, {'no': 6, 'position': 'l'}]}
node8 {'point': array([-0.6 , 0.333]), 'position': 'b3', 'nextnode': [{'no': 4, 'position': 'f'}, {'no': 12, 'position': 'b'}, {'no': 9, 'position': 'r'}]}
node9 {'point': array([-0.333, 0.333]), 'position': 'in', 'nextnode': [{'no': 8, 'position': 'l'}, {'no': 10, 'position': 'r'}]}
node10 {'point': array([0.333, 0.333]), 'position': 'in', 'nextnode': [{'no': 9, 'position': 'l'}, {'no': 11, 'position': 'r'}]}
node11 {'point': array([0.6 , 0.333]), 'position': 'b1', 'nextnode': [{'no': 7, 'position': 'b'}, {'no': 15, 'position': 'f'}, {'no': 10, 'position': 'l'}]}
node12 {'point': array([-0.4, 1. ]), 'position': 'c3', 'nextnode': [{'no': 13, 'position': 'b'}, {'no': 8, 'position': 'f'}, {'no': 13, 'position': 'r'}]}
node13 {'point': array([-0.333, 1. ]), 'position': 'b2', 'nextnode': [{'no': 12, 'position': 'f'}, {'no': 14, 'position': 'b'}, {'no': 12, 'position': 'l'}, {'no': 14, 'position': 'r'}]}
node14 {'point': array([0.333, 1. ]), 'position': 'b2', 'nextnode': [{'no': 13, 'position': 'f'}, {'no': 15, 'position': 'b'}, {'no': 13, 'position': 'l'}, {'no': 15, 'position': 'r'}]}
node15 {'point': array([0.4, 1. ]), 'position': 'c2', 'nextnode': [{'no': 11, 'position': 'b'}, {'no': 14, 'position': 'f'}, {'no': 14, 'position': 'l'}]}
Alright, it's OK. Let's change the shape a little and try it.
Figure 9: When the shape is changed
Good vibes.
Step 3: FDM solver
Matrix Generator
I'm finally here. All I have to do is make a matrix ... but there are two problems here.
Discrete expression inside the calculation area
The first is the discrete equation derived last time,
As shown in the figure below, set the node spacing of left, right, top, and bottom as $ \ Delta_l, \ Delta_r, \ Delta_u. \ Delta_d $. At this time, if the value of $ u $ at each point is $ u_l, u_r, u_u. U_d $ and the center $ u_0 $ is used for Taylor expansion, the result is as follows.
\begin{align}
u_l &= u_0 - \Delta_l \frac{\partial u_0}{\partial x} + \frac{\Delta_l^2}{2} \frac{\partial^2 u_0}{\partial x^2} -\frac{\Delta_l^3}{6} \frac{\partial^3 u_0}{\partial x^3} + \mathcal{O}(\Delta_x^4) \\
u_r &= u_0 + \Delta_r \frac{\partial u_0}{\partial x} + \frac{\Delta_r^2}{2} \frac{\partial^2 u_0}{\partial x^2} +\frac{\Delta_r^3}{6} \frac{\partial^3 u_0}{\partial x^3} + \mathcal{O}(\Delta_x^4) \\
u_d &= u_0 - \Delta_u \frac{\partial u_0}{\partial y} + \frac{\Delta_u^2}{2} \frac{\partial^2 u_0}{\partial y^2} -\frac{\Delta_u^3}{6} \frac{\partial^3 u_0}{\partial y^3} + \mathcal{O}(\Delta_y^4) \\
u_u &= u_0 + \Delta_d \frac{\partial u_0}{\partial y} + \frac{\Delta_d^2}{2} \frac{\partial^2 u_0}{\partial y^2} +\frac{\Delta_d^3}{6} \frac{\partial^3 u_0}{\partial y^3} + \mathcal{O}(\Delta_y^4)
\end{align}
After that, if you add and subtract well and erase the term of the first derivative, you can find it as follows.
\begin{align}
\frac{\partial^2 u}{\partial x^2} &= \frac{2}{\Delta_l \Delta_r(\Delta_l + \Delta_r)}\left( \Delta_r(u_l - u_0) + \Delta_l (u_r - u_0) \right) + \mathcal{O}(\Delta_x)\\
\frac{\partial^2 u}{\partial y^2} &= \frac{2}{\Delta_d \Delta_u(\Delta_d + \Delta_u)}\left( \Delta_u(u_d - u_0) + \Delta_d (u_u - u_0) \right) + \mathcal{O}(\Delta_y)
\end{align}
Unfortunately, the accuracy has become the primary accuracy, but it can't be helped. After that, by substituting it into Poisson's equation, the following equation can be obtained.
\begin{multline}
\frac{2}{\Delta_l \Delta_r(\Delta_l + \Delta_r)}\left( \Delta_r(u_l - u_0) + \Delta_l (u_r - u_0) \right) \\+ \frac{2}{\Delta_d \Delta_u(\Delta_d + \Delta_u)}\left( \Delta_u(u_d - u_0) + \Delta_d (u_u - u_0) \right) = f + \mathcal{O}(\Delta_{x,y})
\end{multline}
It looks awkward, but a computer can calculate it right away.
One thing to note is the size of the $ \ mathcal {O} (\ Delta) $ term. As you can intuitively understand, the larger the difference between the left and right (up and down) intervals, the larger the error. When you actually calculate it, you can see that if the difference between $ \ Delta_l $ and $ \ Delta_r $ is large, then $ \ mathcal {O} (\ Delta_x) $ will increase accordingly. It's not that difficult, and if you're free, try calculating it (those who can understand it by mental arithmetic).
To put it the other way around, it is better for the grid spacing to be the same as much as possible. Things like Figure 5 aren't really good. Actually, it is better to increase the number of nodes or decrease them.
Discrete expression on the boundary
The other problem is more troublesome. Approximate formula of Neumann boundary condition derived last time
\begin{align}
\frac{u_{m+1,n} - u_{m,n}}{\Delta_x} = g \\
\frac{u_{m,n+1} - u_{m,n}}{\Delta_y} = g
\end{align}
Is just an approximation when the boundary is perpendicular to the $ x, y $ axis (that is, partial differential of $ x $ or $ y $), and cannot be used in a slightly oblique case like this time.
Figure 10: Diagonal boundary conditions
It would be nice if there was a node in the direction of the arrow in the above figure, but that is not the case. So, let's go back to the basics and think from the Taylor expansion.
And before that
In other words, somehow the following discretization should be performed.
\begin{align}
u_{r, s}&=\sum_{p=0}^\infty \frac{1}{p!}\left( r \frac{\partial}{\partial x} + s \frac{\partial}{\partial y} \right)^p u_0 \\
&=u_{0,0} + r \frac{\partial u_0}{\partial x} + s \frac{\partial u_0}{\partial y} + \cdots
\end{align}
The second and third terms should be the same as Eq. (1). This can be done by combining two equations with Taylor expansion at two points. If you dare to write it redundantly,
\begin{align}
u_{r_1, s_1}=u_{0,0} + r_1 \frac{\partial u_0}{\partial x} + s_1 \frac{\partial u_0}{\partial y} + \cdots \\
u_{r_2, s_2}=u_{0,0} + r_2 \frac{\partial u_0}{\partial x} + s_2 \frac{\partial u_0}{\partial y} + \cdots
\end{align}
It is sufficient to add or subtract 3 items and 4 items to formula (1). In other words, the following equation should be solved for $ c_1 $ and $ c_2 $.
\begin{align}
r_1c_1 + r_2c_2 = a \\
s_1c_1 + s_2c_2 = b
\end{align}
If we solve for $ c_1 $ and $ c_2 $, we can approximate it as follows.
$ \ mathcal {E} $ represents the error. I didn't know how to write it properly, so I wrote it properly.
Well, it's easy if you can do this. For $ u_ {r_1, s_1} $, $ u_ {r_2, s_2} $, use the two near the boundary. Create a matrix about this and solve it with numpy.linalg.solve.
However, there are 3 nodes near the boundary, and it is a waste to solve with 2 nodes. Of course, if you can solve with two, you can solve with more accuracy if you have three.
In other words
Is established. If $ g $ is a constant, the right side is $ 0 $. As $ g = 0 $
Is. If this equation is made with three points, secondary accuracy can be achieved. For the time being, let's write down the three formulas below.
\begin{align}
u_{r_1, s_1}=u_{0,0} + r_1 \frac{\partial u_0}{\partial x} + s_1 \frac{\partial u_0}{\partial y} + \frac{r_1^2}{2} \frac{\partial^2 u_0}{\partial x^2} + r_1s_1 \frac{\partial^2 u_0}{\partial x \partial y} + \frac{s_1^2}{2} \frac{\partial^2 u_0}{\partial y^2} + \cdots \\
u_{r_2, s_2}=u_{0,0} + r_2 \frac{\partial u_0}{\partial x} + s_2 \frac{\partial u_0}{\partial y} + \frac{r_2^2}{2} \frac{\partial^2 u_0}{\partial x^2} + r_2s_2 \frac{\partial^2 u_0}{\partial x \partial y} + \frac{s_2^2}{2} \frac{\partial^2 u_0}{\partial y^2} + \cdots \\
u_{r_3, s_3}=u_{0,0} + r_3 \frac{\partial u_0}{\partial x} + s_3 \frac{\partial u_0}{\partial y} + \frac{r_3^2}{2} \frac{\partial^2 u_0}{\partial x^2} + r_3s_3 \frac{\partial^2 u_0}{\partial x \partial y} + \frac{s_3^2}{2} \frac{\partial^3 u_0}{\partial y^2} + \cdots
\end{align}
You may wonder if equation (2) can be solved by equation 3 while it is 5 terms, but $ k_1 $ and $ k_2 $ are arbitrary values, so it can be managed. Let's write it down as a system of equations. First, the following terms are created by multiplying each term by a coefficient and adding them.
\begin{equation}
\left\{
\begin{alignedat}{5}
r_1 &c_1& + r_2 &c_2 &+ r_3 c_3&=& a \\
s_1 &c_1& + s_2 &c_2 &+ s_3 c_3&=& b \\
r_1^2&c_1& + r_2^2 &c_2 &+ r_3^2 c_3&=& 2 k_1a \\
r_1 s_1 &c_1& + r_2 s_2 &c_2 &+ r_3 s_3 c_3&=& k_1 b + k_2 a \\
s_1^2&c_1& + s_2^2 &c_2 &+ s_3^2 c_3&=& 2 k_2b
\end{alignedat}
\right.
\end{equation}
If you look at these $ k_1 $ and $ k_2 $ as variables, you will have 5 variables and 5 equations. By rewriting it into a matrix, we can make the following formula:
\begin{equation}
\left[
\begin{matrix}
r_1 & r_2 & r_3 & 0 & 0 \\
s_1 & s_2 & s_3 & 0 & 0 \\
r_1^2 & r_2^2 & r_3^2 & -2a & 0 \\
r_1s_1 & r_2s_2& r_3s_3& -b & -a \\
s_1^2 & s_2^2 & s_3^2 & 0 & -2b \\
\end{matrix}
\right]\left[
\begin{matrix}
c_1 \\
c_2 \\
c_3 \\
k_1 \\
k_2
\end{matrix}
\right] =
\left[
\begin{matrix}
a \\
b \\
0 \\
0 \\
0
\end{matrix}
\right]
\end{equation}
This seems to be solved somehow.
Step 4: Result output
Now, finally, the result is output.
It's the most important thing in terms of showing results to people, and it's a lot of fun for developers to understand. However, if bad results are obtained here, the feeling of despair is great. No, it's unlikely that you'll succeed the first time ...
Potential
Before assuming the result, let's imagine the answer for a moment. The problem was that the upper and lower sides of the trapezoid were 1 and -1, respectively, under the Dirichlet boundary condition. What if it was a square? This is actually exactly the same as a capacitor in the ideal state, that is, the electric field is constant from top to bottom, and the potential is proportional. Then, if the upper side becomes smaller ... So, I wrote an image diagram below.
Figure 10: Image of the solution
The green arrows are the electric fields, and the light rainbow lines represent the contour lines of the potential. If it looks like this, I was able to get the correct answer! For the time being, I tried it with 100x100.
I wrote it as if it was done at once, but it is a result of crushing bugs many times and fixing it with some ingenuity.
Electric field (Flux Density)
Let's also output the electric field. This time, we will calculate the electric field on the grid points from the four surrounding potentials. In other words
\begin{align}
\left. \frac{\partial u}{\partial x} \right|_{m,n} = \frac{1}{2\Delta_x}(u_{m+1, n} - u_{m-1, n}) \\
\left. \frac{\partial u}{\partial y} \right|_{m,n} = \frac{1}{2\Delta_y}(u_{m, n+1} - u_{m, n-1})
\end{align}
Calculate as follows. The result is as follows.
It feels very good. I also wanted to compare it with the theoretical solution, but this time it has become long, so I will stop it. If I feel like it someday ...
Poisson equation class
I created a PoissonEquation class that summarizes all of the above.
#abridgement
if __name__ == "__main__":
#Step 1: Problem
domain = {"shape":"Polygon",
"vertexes":[[-1,-1],[1,-1], [0.4,1],[-0.4,1]],
"bc":{
"bc": [
{"bctype":"Dirichlet", "constant":-1},
{"bctype":"Neumann", "constant":0},
{"bctype":"Dirichlet", "constant":1},
{"bctype":"Neumann", "constant":0}
],
"priority":[0,2,1,3]
}
}
source = {"source": 0}
poisson = PoissonEquation(domain,source)
# Step 2: generateGrid
grid = {"type":"Cartesian", "div":[25,25]}
poisson.generateGrid(grid, "Node")
# Step 3: Solve
method = "FDM"
poisson.solve(method)
# Step 4: Result
poisson.result()
poisson.plot("Potential")
poisson.plot("DensityFlux")
Summary
This time, we designed and developed the FDM solver for Poisson's equation in the following four parts.
--Problem --Grid class (Grid) --Solver class --Result class
The program I made this time is on the following URL.
[https://github.com/atily17/PoissonEquation]
I think it works if you have numpy and matplotlib. I haven't done much testing, so I don't guarantee that it will work properly. Perhaps, if you change the problem conditions a little, it will collapse at once. If you would like to fix a bug, please let me know. I'll fix it when I feel like it.
TODO The next article will talk about the finite element method (FEM). However, there are still some areas where FDM is not enough, so let's write only the keywords here. I will write an article when I feel like it someday.
--Coordinate systems other than Cartesian coordinate systems --Polar coordinate system --Curvilinear coordinate system --Analysis method of boundary conditions at infinity - Boundary Relaxation Method --High-order precision FDM --Formulation of quaternary accuracy --Compact Stencil Fourth order accuracy --Unstructured grid FDM --Subgrid
gallery
U-shaped example
narrow
(It feels a little unstable ...)
Charge in the middle
