I would like to solve the following problem $$ \vec{v}(x,y)= k\, \nabla \theta(x,y) $$ with respect to the unknown function $\theta$. Parameter $k$ is just a real constant quantity.

I have two matrices where I have stored the values of $v_x(x,y)$ and $v_y(x,y)$, which are the components of the 2D vector field $\vec{v}$.


1 Answer 1


You can convert this equation into a Poisson equation and use standard numerical methods to solve it like finite difference to obtain $\theta(x,y)$. If you take a divergence from both side:

$$\nabla^{2} \theta = \frac{1}{k}\nabla \cdot \vec{v}$$

Since, you know the $\vec{v}$ you can approximate its divergence as:

$$\nabla \cdot \vec{v} = \frac{\partial v_{x}}{\partial x} + \frac{\partial v_{y}}{\partial y} = \frac{v_{x}(x+\Delta x,y) - v_{x}(x-\Delta x,y)}{2 \Delta x} + \frac{v_{y}(x,y+\Delta y) - v_{y}(x,y-\Delta y)}{2 \Delta y}$$

The only concern here is that how you want to handle the boundaries. If it makes sense for simplest case you can assume a periodic boundary condition here. Now you have $\nabla \cdot \vec{v}$ stored in an array and for the rest I assume it is known already. In order to find $\theta$ from the defined Poisson equation, you need to know the boundary conditions of $\theta$ but again I assume it's defined as periodic boundary conditions on all boundaries. So, now you can discretize the Poisson equation as:

$$\frac{\theta(x+\Delta x,y)+\theta(x-\Delta x,y)-2\theta(x,y)}{\Delta x^{2}} + \frac{\theta(x,y+\Delta y)+\theta(x,y-\Delta y)-2\theta(x,y)}{\Delta y^{2}} = \frac{1}{k} \mathrm{div}(x,y)$$

Where $\mathrm{div}(x,y)$ is the array that you stores the values of $\nabla \cdot \vec{v}$ in it beforehand. Now you have a standard discretized Poisson equation that you can convert it to a system of linear equations. In fact, I want to convert above system of discretized equations into a matrix form as:

$$A x = b$$

Matrix $A$ is called mass matrix and if you have a 2D grid of $N \times M$ points, $A$'s dimension is $NM \times NM$. Now, let's say I can map $x$ and $y$ grid into matrix coordinate of $i$ and $j$. In fact, you have $NM$ points in your $x$ and $y$ grid and you can assign an ID from 0 to $NM - 1$ to each of them. This would be a function as:

$$F(x,y) = i$$

In fact you want to flatten a 2D grid into 1D as:

$$F(x,y) = x + N y$$

You give $x$ and $y$ coordinate and it just gives you the id of the point. Then I can assemble $A$ and $b$ as:

$$A_{F(x,y),F(x,y)} = -2 \Bigg ( \frac{1}{\Delta x^{2}} + \frac{1}{\Delta y^{2}} \Bigg)$$

$$A_{F(x+\Delta x,y),F(x,y)} = \frac{1}{\Delta x^{2}}$$

$$A_{F(x-\Delta x,y),F(x,y)} = \frac{1}{\Delta x^{2}}$$

$$A_{F(x,y + \Delta y),F(x,y)} = \frac{1}{\Delta y^{2}}$$

$$A_{F(x,y - \Delta y),F(x,y)} = \frac{1}{\Delta y^{2}}$$

$$b_{F(x,y)} = \frac{1}{k} \mathrm{div}(x,y)$$

Keep in mind you need to treat boundaries specially based on your boundary condition to assemble $A$ and $b$. Finally you solve this matrix equation with standard matrix solvers available in Python or any other programming language that you use for your problem.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.