Problem description: I want to numerically solve system of two time-independent partial differential equations (pde) of the following simple form

$$\frac{\partial u(x,y)}{\partial x} = f_1(x,y,u),$$

$$\frac{\partial u(x,y)}{\partial y} = f_2(x,y,u),$$

where $u(x,y)$ is the unknown we want to solve for, $x$ and $y$ are the spatial variables and $f_1$,$f_2$ are some smooth (or at least continuously differentiable) and possibly nonlinear functions satisfying $\frac{\partial f_1}{\partial y}=\frac{\partial f_2}{\partial x}$. More compactly

$$\nabla u(x,y) = f(x,y,u),$$

with $\nabla:=[\frac{\partial~ \cdot}{\partial x}, \frac{\partial~ \cdot}{\partial y}]^T$ and $f:=[f_1,f_2]^T$.

I want to solve it on a rectangular domain (defined let’s say by $x$,$y$ for which $W>x>0$ and $H>y>0,~ W,H \in \mathbb R$) for an initial given boundary condition of $u(0,0)=c,~c \in \mathbb R$.

What I have tried so far: $u(x,y)$ can be imagined like a surface. We stand at the initial point and the partial derivatives tell us how much the surface goes upwards or downwards in the corresponding directions. In 1D case, when we have $\frac{du}{dx}=f(x,u)$, an ODE solver (like ode45 in Matlab) can be used (we can interpret the time $t$ as the spatial variable $x$ without no harm here). I thought, whether there is some way how to extend this approach to the 2 dimensional case. The solution that works to some extent is to solve for $u$ along one of the edges of the domain (e.g. find $u$ on $W>x>0,~y=0$ using ode45 as mentioned above) and then to use points of this solution as starting values for other set of ode's (now going in the perpendicular direction). In other words, if the goal is to find the values of $u$ on a regular grid, then we first find its values along the first row of the grid and then use the result as a starting point for ode's going along each of the columns and thus getting values for the whole grid. Although such a solution gives an idea how the surface $u(x,y)$ looks like, it is, however, not very precise as the individual solutions for rows/columns of the grid get apart from each other as we are moving away from the edge used as a set of initial values.

Note: I have also tried chebfun toolbox for Matlab and pde toolbox, but it seems that they only support time-dependent equations.

  • $\begingroup$ $f_1$ and $f_2$ cannot both be arbitrary, you require $\partial_y f_1 = \partial_x f_2$. $\endgroup$ Commented Feb 4, 2017 at 11:02
  • $\begingroup$ True, I editted the description.They are arbitrary just in the sense that I am not interested in solution only for some specific choice of $f_1$,$f_2$. I don't know how to better describe these functions, but one can imagine to get them by proceeding in opposite direction (from the result) - i.e. by differentiating some smooth surface. $\endgroup$
    – m_eps
    Commented Feb 4, 2017 at 13:01
  • $\begingroup$ Can you give an example $(f_1, f_2)$ pair that gave you bad results? I would like to test it out with some simple method. $\endgroup$ Commented Feb 4, 2017 at 13:51
  • $\begingroup$ Unfortunatelly, I can not do that right now, since I tried to post as abstracted version of the problem as possible. In my original settings things are more complicated - $f_1$ and $f_2$ are not given analytically, but are evaluated based on some quite large dataset. I am pretty sure, however, that all the conditions imposed on these functions are satisfied. Nevertheless, if any of the methods suggested below won't solve my issue, I will try to find a way to construct some minimal example to test the methods on. $\endgroup$
    – m_eps
    Commented Feb 5, 2017 at 11:49
  • $\begingroup$ How about taking the (numerical) divergence of your equations and solving the the resulting poisson-type equation? You could integrate along edges first to get boundary conditions. $\endgroup$ Commented Feb 6, 2017 at 22:57

3 Answers 3


The approach you describe should work, and if you let the step size of your ODE solvers go to zero, you will recover the exact solution of the problem. You can of course also start by first going in $y$ direction instead of the $x$ direction. Or, you can use a marching front algorithm in which you compute, for example, $u(0,h)$ and $u(h,0)$ using one-dimensional steps, and then $u(h,h)$ by averaging the solution you would get by starting from either $u(0,h)$ or $u(h,0)$ and then going in the perpendicular solution. This is likely going to be more accurate. You would then move out from the points you have to the next lattice points on a mesh with points spaces by distances $h$ in each direction, averaging whenever you can.

  • $\begingroup$ Thank you very much for your suggestion, I'll definitely try it. Although I am indeed interested in a solution on some lattice of points, I hoped that there would (for such simply looking problem) exist a method avoiding the use of fixed step size; some algorithm that adaptivelly changes the step size as it is required to achieve some given solution accuracy (like ode45 does in one dimensional case) and thus save computational time. $\endgroup$
    – m_eps
    Commented Feb 5, 2017 at 12:09
  • $\begingroup$ Yes. All I wanted to point out was that you can do a front marching method. You may want to search the literature for the Fast Marching Method which is typically used for the Eikonal equation, but would be equally applicable to your case. I am certain that there are papers for non-uniform meshes. $\endgroup$ Commented Feb 6, 2017 at 16:08

u have a system of differential equations to be solve.

first u are trying to map a surface into a cube. u can break the equation into delta u(x,y)=f(x,y)=phi(u(x,y).

so the set of equations are delta u(x,y)=f(x,y)....................equation 1.. .........................................delta u(x,y)=phi (u(x,y))=f(x,y).......equation 2.............................................with the condition .........f(x,y) where x=f1, y=f2.............. means..................(partial x)/(partial y)=(partial y)/(partial x) ...................................................................equation 3

the system of equations will be { partial u/partial x .{f ={0 partial u/partial y . f = 0

                       partial u/partial x . f                  = 0
                       partial u/partial y . f                  = 0
                       partial x/partial y . partial y/partial x  = 0 
                       partial y/partial x . partial x/partial y} = 0}

that's {delta u(x,y), partial (x,y)}^T dot {f, f, partial (y,x))^T={0}

u just have to solve the equations simultaneously.


How about the following approach?

Step 1: Determine the values on the $x=0$ and $y=0$ lines: $$ u_{i+1,0} = u_{i,0} + f_{i,0}^1\Delta x + \mathcal{O}(\Delta x^2) \\ u_{0,j+1} = u_{0,j} + f_{0,j}^2\Delta y + \mathcal{O}(\Delta y^2) $$

Step 2: Use Taylor series expansion for all the other points: $$ u_{i+1,j+1} = u_{i,j} + f_{i,j}^1\Delta x + f_{i,j}^2 \Delta y + \mathcal{O}( \Delta x^2, \Delta x\Delta y, \Delta y^2) $$

Note that for all the "interior" points, this involves both $f^1$ and $f^2$ contributing to $u_{i,j}$ thus removing directional bias.


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.