Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Discretization of equation

  1. Mar 10, 2015 #1
    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
     
  2. jcsd
  3. Mar 11, 2015 #2
    Are you going to solve it numerically ?
     
  4. Mar 11, 2015 #3
    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?
     
  5. Mar 12, 2015 #4
    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.
     
  6. Mar 12, 2015 #5
    Do you have some efficient method/algorithm to solve the resulting linear system ?
     
  7. Mar 12, 2015 #6
    Yes I do. However, my question was related to the discretization of the equation.
     
  8. Mar 12, 2015 #7
    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.
     
  9. Mar 12, 2015 #8
    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
     
  10. Mar 12, 2015 #9
    I was just curious. I don't think that is an often used stencil. How do you handle the boundary conditions?
     
  11. Mar 16, 2015 #10
    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook