# Discretization of equation

1. Mar 10, 2015

### kolmog

Hello,

I am trying to solve the following differential equation in a tri-dimensional grid with grid spacing of hx , hy , and hz along the x, y, and z coordinates.
$$\nabla\cdot(\epsilon\nabla\phi)=-4\pi\rho$$
Here, ε is a scalar quantity that is a function of x, y, and,z, while Φ is a potential.
The left-hand side of the equation can be rewritten as follows

$$\nabla\cdot(\epsilon\nabla\phi)=\frac{\partial\epsilon}{\partial x}\frac{\partial\phi}{\partial x}+\frac{\partial\epsilon}{\partial y}\frac{\partial\phi}{\partial y}+\frac{\partial\epsilon}{\partial z}\frac{\partial\phi}{\partial z}+\epsilon\left(\frac{\partial^2\phi}{\partial x^2}+\frac{\partial^2\phi}{\partial y^2}+\frac{\partial^2\phi}{\partial z^2}\right)$$
According to the weights given in the article entitled "Generation of Finite Difference Formulas on Arbitrarily Spaced Grids" by Bernt Fornberg, the first derivative of the function f(x,y,z) along the x coordinate (up to sixth order) is approximated as
$$\frac{\partial f}{\partial x}=\frac{3}{4}f(x+h_x)-\frac{3}{4}f(x-h_x)-\frac{3}{20}f(x+2h_x)+\frac{3}{20}f(x-2h_x)+\\+\frac{1}{60}f(x+3h_x)-\frac{1}{60}f(x-3h_x)$$
For the y and z coordinates the result would be the same changing 'x' by 'y' or 'z'.
For the second derivative along the x coordinate:
$$\frac{\partial^2f}{\partial x^2}=-\frac{49}{18}f(x,y,z)+\frac{3}{2}f(x+h_x)+\frac{3}{2}f(x-h_x)-\frac{3}{20}f(x+2h_x)-\frac{3}{20}f(x-2h_x)+\\+\frac{1}{90}f(x+3h_x)+\frac{1}{90}f(x-3h_x)$$
From now on, we use the following notation:
\frac{\partial\epsilon}{\partial x}=\epsilon'_x\hspace{1cm}\frac{\partial\epsilon}{\partial y}=\epsilon'_y,\hspace{1cm}\frac{\partial\epsilon}{\partial z}=\epsilon'_z\\\phi(x+h_x,y,z)=\phi_x,\hspace{1cm}\phi(x-h_x,y,z)=\phi_{-x}\hspace{1cm}\\\phi(x+2h_x,y,z)=\phi_{2x},\hspace{1cm}\phi(x-2h_x,y,z)=\phi_{-2x}
\\\phi(x+3h_x,y,z)=\phi_{3x},\hspace{1cm}\phi(x-3h_x,y,z)=\phi_{-3x}
\\\phi(x,y+h_y,z)=\phi_y,\hspace{1cm}\phi(x,y-h_y,z)=\phi_{-y}\hspace{1cm}\\\phi(x,y+2h_y,z)=\phi_{2y},\hspace{1cm}\phi(x,y-2h_y,z)=\phi_{-2y}
\\\phi(x,y+3h_y,z)=\phi_{3y},\hspace{1cm}\phi(x,y-3h_y,z)=\phi_{-3y}
\\\phi(x,y,z+h_z)=\phi_z,\hspace{1cm}\phi(x,y,z-h_z)=\phi_{-z}\hspace{1cm}\\\phi(x,y,z+2h_z)=\phi_{2z},\hspace{1cm}\phi(x,y,z-2h_z)=\phi_{-2z}
\\\phi(x,y,z+3h_z)=\phi_{3z},\hspace{1cm}\phi(x,y,z-3h_z)=\phi_{-3z}

Then, the left-hand side of the equation that we want to solve is written as:
\nabla\cdot(\epsilon\nabla\phi)=-\frac{49\epsilon}{18}\left(\frac{1}{h_x^2}+\frac{1}{h_y^2}+\frac{1}{h_z^2}\right)\phi(x,y,z)+\left[\frac{3\epsilon'_x}{4h_x}+\frac{3\epsilon}{2h_x^2}\right]\phi_x+\left[\frac{3\epsilon}{2h_x^2}-\frac{3\epsilon'_x}{4h_x}\right]\phi_{-x}+\\+\left[\frac{3\epsilon'_y}{4h_y}+\frac{3\epsilon}{2h_y^2}\right]\phi_y+\left[\frac{3\epsilon}{2h_y^2}-\frac{3\epsilon'_y}{4h_y}\right]\phi_{-y}+\left[\frac{3\epsilon'_z}{4h_z}+\frac{3\epsilon}{2h_z^2}\right]\phi_z+\left[\frac{3\epsilon}{2h_z^2}-\frac{3\epsilon'_z}{4h_z}\right]\phi_{-z}

\\-\left[\frac{3\epsilon'_x}{20h_x}+\frac{3\epsilon}{20h_x^2}\right]\phi_{2x}+\left[\frac{3\epsilon'_x}{20h_x}-\frac{3\epsilon}{20h_x^2}\right]\phi_{-2x}-\left[\frac{3\epsilon'_y}{20h_y}+\frac{3\epsilon}{20h_y^2}\right]\phi_{2y}+\\+\left[\frac{3\epsilon'_y}{20h_y}-\frac{3\epsilon}{20h_y^2}\right]\phi_{-2y}-\left[\frac{3\epsilon'_z}{20h_z}+\frac{3\epsilon}{20h_z^2}\right]\phi_{2z}+\left[\frac{3\epsilon'_z}{20h_z}-\frac{3\epsilon}{20h_z^2}\right]\phi_{-2z}\\+\left[\frac{\epsilon'_x}{60h_x}+\frac{\epsilon}{90h_x^2}\right]\phi_{3x}+\left[\frac{\epsilon}{90h_x^2}-\frac{\epsilon'_x}{60h_x}\right]\phi_{-3x}+\left[\frac{\epsilon'_y}{60h_y}+\frac{\epsilon}{90h_y^2}\right]\phi_{3y}+\\+\left[\frac{\epsilon}{90h_y^2}-\frac{\epsilon'_y}{60h_y}\right]\phi_{-3y}+\left[\frac{\epsilon'_z}{60h_z}+\frac{\epsilon}{90h_z^2}\right]\phi_{3z}+\left[\frac{\epsilon}{90h_z^2}-\frac{\epsilon'_z}{60h_z}\right]\phi_{-3z}

For the right-hand side, according to the same article by Bernd Fornberg, it is just -4πρ(x,y,z).

Then, one has just to solve the linear system of equations,
$$A\phi=f$$ where 'A' is the matrix given by the discretization shown above for the left-hand side (with the appropriate boundary conditions for the points located at the limits in the x, y, and z coordinates) and f is the vector given -4πρ.

I am checking if everything was well done. Can somebody read this post and tell me if the discretization is ok? I am pretty sure that everything has been calculated well, but I want to be sure.

Thanks

2. Mar 11, 2015

### soarce

Are you going to solve it numerically ?

3. Mar 11, 2015

### Strum

Why haven't you used any of the f(x+h,y+h,z) and sign permutations as well as f(x+h,y+h,z+h) with sign permutations?

4. Mar 12, 2015

### kolmog

I am going to solve it numerically. In fact, I think I have solved it but I am not sure if the discretization was properly done.

5. Mar 12, 2015

### soarce

Do you have some efficient method/algorithm to solve the resulting linear system ?

6. Mar 12, 2015

### kolmog

Yes I do. However, my question was related to the discretization of the equation.

7. Mar 12, 2015

### kolmog

With respect to 'Strum' comment. I just use the values at x, x±hx, x±2hx, x±3hx (and the same for y and z) because I use the formula given in the article by Bernt Fornberg, that just uses the values at these positions for the derivatives of a function. If you wanted, you could have used another definition for the derivatives of a function, but in this case I have adopted this approach.

8. Mar 12, 2015

### Staff: Mentor

Why are you going to 6th order? Isn't 2nd order adequate. It's so much more simple, and application of the boundary conditions is much easier. To check your discretization, why don't you just solve it with 2nd order approximations, and then compare the results with the 6th order results?

Chet

9. Mar 12, 2015

### Strum

I was just curious. I don't think that is an often used stencil. How do you handle the boundary conditions?

10. Mar 16, 2015

### kolmog

Boundary conditions are periodic. They are easy to implement.

With respect to 'Chestermiller' comment. I have also implemented the 2nd order discretization with similar results.