This is an extended comment rather than an answer. Results are for version 12.1.0 for Microsoft Windows (64-bit) (March 14, 2020).
NDSolve
attempts to solve this system of PDEs by Method -> "FiniteElement"
, and I surmise that it does so by linearizing the nonlinear terms in the PDEs about some initial guess and then repeating the process until the solution converges, using FindRoot
with Method -> "AffineCovariantNewton"
. The initial guess can be specified using the option InitialSeeding
, and the symbolic solution given in the question should be a good initial guess.
sol = NDSolve[{D[sigmoid[H[x, y]], y] == D[sigmoid[V[x, y]], x],
D[H[x, y], x] + D[V[x, y], y] == rho0[x, y],
H[0, y] == H[1, y] == V[x, 0] == V[x, 1] == 0}, {H, V}, {x, 0, 1}, {y, 0, 1},
InitialSeeding -> {H[x, y] == (x - x^2)/2, V[x, y] == (y - y^2)/2}] // Flatten
NDSolve
twice returns the warning message
NDSolveFEM
InitializePDECoefficients::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help.
but then returns the correct answer, as can be determined from
Plot3D[Evaluate[(D[sigmoid[H[x, y]], y] - D[sigmoid[V[x, y]], x]) /. sol], {x, 0, 1}, {y, 0, 1}]
Plot3D[Evaluate[(D[H[x, y], x] + D[V[x, y], y] - rho0[x, y]) /. sol], {x, 0, 1}, {y, 0, 1}]
which display noise with amplitudes of order 5 10^-15
. Unfortunately, providing an initial guess differing only slightly from that given above
InitialSeeding -> {H[x, y] == (x - x^2)/2.01, V[x, y] == (y - y^2)/2.01}
yields the additional, undocumented error message
FindRoot::stfail: The method AffineCovariantNewton failed to compute the next step.
The resulting solution again gives noise for the first plot, but for the second plot gives the surprisingly smooth
Increasing resolution with
Method -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> .00001}}
is slower and returns essentially the same results. I also tried running
SetOptions[FindRoot, Method -> {"AffineCovariantNewton", "BroydenUpdates" -> False}]
before NDSolve
, which documentation suggested could improve accuracy. Instead, NDSolve
produced numerous error messages before returning unevaluated. I also tried
SetOptions[FindRoot, Method -> "Newton"]
but after about 30 minutes without an answer I terminated the calculation. Perhaps, other readers will find my unsuccessful attempts useful in exploring this deceptively simple question.
pdsolve([pde1,pde2,bc],[H(x,y),V(x,y)])
$$ \left\{ H \left( x,y \right) =-{\frac {{x}^{2}}{2}}+{\frac {x}{2}},V \left( x,y \right) =-{\frac {{y}^{2}}{2}}+{\frac {y}{2}} \right\} $$ current version of MathematicaDSolve
can not solve this analytically. Hopefully V 12.2 it will. $\endgroup$restart; rho0:= (x,y)->1 - x - y; sigmoid:= x->tanh(x); pde1 := diff(sigmoid(H(x, y)), y) = diff(sigmoid(V(x, y)), x); pde2 := diff(H(x, y), x) + diff(V(x, y), y) = rho0(x, y); bc := H(0, y) = 0, H(1, y) = 0, V(x, 0) = 0, V(x, 1) = 0; pdsolve([pde1,pde2,bc],[H(x,y),V(x,y)])
$\endgroup$