Finite difference methods for multidimensional coupled equations

Question

My knowledge of finite difference is very basic so this could be very trivial. I've seen how multidimensional finite difference works for say fluid equations, but they are also dealing with a single variable. How would one solve a coupled multidimensional equation, for example of the form

$$a \nabla a = \nabla b,$$

in two dimensions? Meaning we have the system of equations

$$a \, \partial_x a = \partial_x b,$$

$$a \partial_y a = \partial_y b.$$

Using a simplistic finite difference method I can solve each equation independently

$$a(x+\Delta x, y) = \frac{b(x+\Delta x,y) - b(x,y)}{a(x,y)} + a(x,y),$$

and

$$a(x, y+\Delta y) = \frac{b(x, y+\Delta y) - b(x,y)}{a(x,y)} + a(x,y).$$

However if I want to get to $a(x+\Delta x, y+\Delta y)$, there are two options either first moving along $x$ and then $y$ or vice versa. In this case I am not guaranteed to get the same result, meaning the answer is path dependent which is not desirable. What is the common practice to deal with this?

Edit: I'll compute the values over a square grid with initial conditions $a(0,0) = a_0$ and $b(x,y)=x-y$.

My purpose is to solve this equation so that I may initialize the grid in $a$. This will be implemented in a program I am writing. The equation given is just a no-frills re-expression of what I'm attempting to do. My true equation is solving an isentropic hydrostatic atmosphere with a complicated potential, $\nabla \Phi$, for the density structure, $\rho$. The equation reduces to $$\gamma \, K \, \rho^{\gamma-2} \, \vec{\nabla} \rho = - \vec{\nabla} \Phi$$, where $K$ and $\gamma$ are some constants.

Answer 1

Since your actual problem is totally different from your initial one I will explain what you should do in what follows, using three different approaches:

Your equation to be solved is: $$\gamma K\rho^{\gamma-2}\vec{\textrm{grad}}\rho=\vec{\textrm{grad}}\Phi \tag{*}$$ $\textbf{Fast solution}$

You can notice from $(*)$ that can be simply rewritten: $$\vec{\textrm{grad}}\left(\frac{\gamma K }{\gamma-1}\rho^{\gamma-1}-\Phi\right)=0\tag{a}$$ Therefore your solution for $\rho$ from $(a)$ is: $$\rho = \left(C+\frac{\gamma -1 }{\gamma K}\Phi\right)^{\frac{1}{\gamma-1}}$$ Where $C$ is a constant that can be adjusted using an arbitrary known value for $\rho$.

$\textbf{Equation as "it is''}$

Compute explicitly the gradient of your potential $\Phi$ and denote by $\Phi_x$ and $\Phi_y$ its partial derivatives.

Your discretised equations, following the simplified equation $C(\rho)\vec{grad}{\rho}=(\Phi_x,\Phi_y)^{T}$ would be: $$C_{i,k}\left(\rho_{i+1,k}-\rho_{i,k}\right)=\Delta x\, \Phi_{x\,i,k} \quad i=1,...,N_x-1 \quad k=1,...,N_y\tag{b1}$$ $$C_{i,k}\left(\rho_{i+1,k}-\rho_{i,k}\right)=\Delta y\, \Phi_{x\,i,k}\quad i=1,...,N_x \quad k=1,...,N_y-1\tag{b2}$$

Since we have the equation $(*)$ we know that $\rho$ is compatible with the right hand side, $i.e$ rho is a true potential and coincides with $\Phi$.

Since we have two equations: $(b1)$ and $(b2)$ (that are not independent) to solve one variable $\rho$, we have an overdetermined system of equations. But we know from the fact that $\rho$ is a true potential that the integration along $x$ (equation $(b1)$) or along $y$ (equation $(b2)$) must result in the same function. The number of equations given by $(b1)$ is $(N_x-1)N_y$ and the given by $(b2)$ is $N_x(N_y-1)$. Therefore the total number of equations $2N_xN_y-N_x-N_y$ is bigger that the true number of unknowns $N_xN_y$.

Build the total system and drop off the equations that give nothing new. Do not forget to impose the given value that confers uniqueness to your problem replacing an additional equation.

$\textbf{Well posed problem}$

For the kind of problem: find the potential $\phi$ that verifies $\vec{\textrm{grad}}\phi = \vec{f}$ one can set the Poisson problem that solves:

$$\triangle \phi = \textrm{div}\vec{f} \tag{**}$$

With the boundary condition:

$$\frac{\partial\phi}{\partial n} = \vec{f}\cdot\vec{n}\tag{***}$$ Where $\vec{n}$ is the boundary unitary outer normal.

The problem $(**)$ and $(***)$ is numerically preferred compared with the second procedure, given by $(b1)$ and $(b2)$ and $\phi$ can be obtained with better precision choosing an appropiate numerical scheme, $e.g$ in your case central finite differences (compute analitically the derivatives of the RHS):

$$\frac{\phi_{i+1,k}-2\phi_{i,k}+\phi_{i-1,k}}{\Delta x^2}+\frac{\phi_{i,k+1}-2\phi_{i,k}+\phi_{i,k-1}}{\Delta y^2}=\left(\frac{\partial f_1}{\partial x}\right)_{i,k}+\left(\frac{\partial f_2}{\partial y}\right)_{i,k} $$

With the discretised boundary conditions

$$\frac{\phi_{N_x+1,k}-\phi_{N_x-1,k}}{2\Delta x} = \left(f_1\right)_{N_x,k}$$ $$\frac{\phi_{i,N_y+1}-\phi_{i,N_y-1}}{2\Delta y} = \left(f_1\right)_{i,N_y}$$

Answer 2

Next time, provide please the complete problem, i.e. with boundary conditions. This simplify everything and allow us to provide you the best answer at once.

Supposing your mesh has the same points in $x$ and $y$ axes, namely $n$, we have $2n^2$ unknowns for a discrete vector $\vec{w}$ in which, for the sake of simplicity, the unknowns are ordered by pairs $w_i = (a,b)_i$.

The goal of the method is to find the solution of the following system: $$M(\vec{w})\vec{w} = \vec{b}$$ Where $M$ is a matrix that depends in a nonlinear way on $\vec{w}$ (your PDE is nonlinear in $a$ and therefore $M=M(a)$), and $b$ is the vector in which the Boundary conditions would be applied.

First provide an easy ordering for your points through the mesh, i.e. begin to count from the lower left corner, from left to right and after the row is competed, keep going to the upper one, reading from left to right...

Therefore the point indices: $$(1,1) = 1, (1,2) = 2,..., (n,1) = n, (2,1) = n+1, (2,2) = n+2, ...(n,n) = n^2$$

To index $1$ corresponding to the point $(1,1)$ we would have a pair of values $(a_1,b_1)$ and so on...

The notation for a variable $\phi$ would be considered here: $$\phi(x,y) = \phi_{i,k}\qquad \phi(x+\Delta x,y) = \phi_{i+1,k} \qquad \phi(x,y-\Delta y) = \phi_{i,k-1} $$ Your discretised equations read in index notation: $$a_{i,k}\left(a_{i+1,k}-a_{i,k}\right)-\left(b_{i+1,k}-b_{i,k}\right) = 0 \tag{1}$$ $$a_{i,k}\left(a_{i,k+1}-a_{i,k}\right)-\left(b_{i,k+1}-b_{i,k}\right) = 0 \tag{2}$$ You must arrange with care the matrix $M$ for the vector $\vec{w}$ given by: $$\vec{w} = (a_1,b_1,a_2,b_2,...,a_n, b_n, ..., a_{n^2}, b_{n^2})^{T}$$

For the $2n$th and $2n+1$th row of $M$ corresponding to $(1)$ and $(2)$ respectively you'll have:

$$M(2n,2n) = -a_{n},\quad M(2n,2n+1) = 1, \quad M(2n,2n+2) = a_n,\quad M(2n,2n+3) = -1$$ $$M(2n+1,2n) = -a_{n},\quad M(2n+1,2n+1) = 1, \quad M(2n+1,4n+2) = a_n,\quad M(2n+1,4n+3) = -1$$

Hope I have not made any mistake in the indices.

After that, you will must provide an initial guess for $a(x,y)$ to evaluate the matrix $M$ at the first iteration, i.e., you will solve the following (simplest) iterative numerical scheme: $$M(\vec{w}^{n})\vec{w}^{n+1}=\vec{b}$$

OF COURSE it is extremely important to set an appropiate $\vec{b}$ and modified $M$ entries in order to fulfill the Boundary conditions.