A nonlinear system with special structure

Suppose we have an $n \times n$ uniform grid, covering $[-1,1] \times [-1,1]$ (typically, n $\approx$ 500). We have a smooth, differentiable function $z(x,y)$ that we want to determine on the nodes of this grid. At every point, the function $z(x,y)$ satisfies a relation of the form:

$\displaystyle\frac{az + b}{cz + d} = \displaystyle\frac{z_x}{z_y}$

where $z_x$ and $z_y$ are the partial derivatives. Note that the values of $a, b, c, d$ are samples from four known smooth functions (so they are not constant at every point on the grid). How can we solve systems of equations with such structure?

[EDIT 1] (Boundary Condition): The ratio $\displaystyle\frac{z_x}{z_y}$ is known at the boundary. Thus, $z$ is known at the boundary.

[EDIT 2] Am I correct in understanding that this is a quasilinear first-order PDE and the method of characteristics will solve it up to level sets of z? Is there a robust way to solve such PDEs in the presence of noise and isolated singularities in a, b, c, d?

-
If A/B = C/D you can reduce it to an advection equation. –  Taylor Dupuy Aug 30 '11 at 16:49
What are your boundary conditions? Are you sure that the equation has a unique solutions? –  Paul Tupper Sep 1 '11 at 3:41
Hi Paul, we can assume any reasonable boundary conditions. For instance, the value of z can be known at the boundary. –  user17473 Sep 2 '11 at 22:35
I'm not sure I am clear about the set up. Are a,b,c,d, smooth functions whose values you know only at the grid points? And is the relation you describe only satisfied at the grid points? If the latter, then I don't think the method of characteristics makes much sense, since we don't know anything about the invariants away from those points. –  Nilima Nigam Sep 27 '11 at 2:07
Yes, (a,b,c,d) are samples from smooth functions, I have now edited the question to say so. We are looking for a solution that is robust to noise and will be great if we can handle isolated discontinuities in those functions. It is my understanding that the method of characteristics is not well-suited for such situations? –  user17473 Sep 27 '11 at 18:08

(1) Discretize the partial derivatives using finite differences on the grid. At every point on the grid your PDE will give you a nonlinear equation. You will have a total of $n^2$ equations in $n^2$ unknowns. But each equation will only contain a small number of unknowns. You can use a variant Newton's method to solve these equations. But you have to work with the restrictions that you probably don't want to code up the Jacobian for the system, and you could not use a direct method to calculate the solution even if you did. I would recommend looking at a Jacobian-Free Newton-Krylov method.
(2) Exploit the method of characteristics. Rewrite your equation as $$(cz+d) z_x - (az+b) z_y = 0.$$ This gives characteristic equations $$\frac{dx}{dt}=cz+d, \ \ \ \frac{dy}{dt}=a z+b, \ \ \ \frac{dz}{dt}= 0.$$ There are numerical methods for solving nonlinear hyperbolic equations exploiting characteristics. Sethian and Vladimirsky have a nice one. Your problem does not quite fit into their scheme but their paper might help give you ideas.
I'd have to agree- the problem, as stated, may not admit well-defined solutions. Set $c=0, d=1,a=-1,b=0$, and one finds Burger's equations. I can think of boundary data on the edges of the square which would lead to a lack of well-posedness. A numerical attack via finite differences would be reasonable to try. –  Nilima Nigam Sep 28 '11 at 1:41