Numerically solving $\nabla u(x,y) = f(x,y,u)$ on a rectangular domain having initial value of $u$ at some point

Question

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.

Answer 1

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.

Answer 2

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.

Answer 3

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.