- #1
kolmog
- 11
- 0
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.
\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
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.
\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