Is the Discretization of this Differential Equation Accurate?

[kolmog]

Hello,

I am trying to solve the following differential equation in a tri-dimensional grid with grid spacing of h_x , h_y , and h_z along the x, y, and z coordinates.
\begin{equation}\nabla\cdot(\epsilon\nabla\phi)=-4\pi\rho\end{equation}
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

\begin{equation}\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)\end{equation}
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
\begin{equation}\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)\end{equation}
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:
\begin{equation}\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)\end{equation}
From now on, we use the following notation:
\begin{equation}\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}
\end{equation}

Then, the left-hand side of the equation that we want to solve is written as:
\begin{equation}\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}\end{equation}

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,
\begin{equation}A\phi=f\end{equation} 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

 

[soarce]

Are you going to solve it numerically ?

 

[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?

 

[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.

 

[soarce]

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

 

[kolmog]

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

 

[kolmog]

With respect to 'Strum' comment. I just use the values at x, x±h_x, x±2h_x, x±3h_x (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.

 

[Chestermiller]

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

 

[Strum]

With respect to 'Strum' comment. I just use the values at x, x±h_x, x±2h_x, x±3h_x (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.

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

 

[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.