I need to solve for function H[x,y] and V[x,y] on square {x,0,1},{y,0,1}, obeying the following PDEs and boundary conditions: enter image description here

wherein $\rho_0(x,y)$ is simply defined as 1-x-y, and $\sigma(t)=tanh(t)$ is the sigmoid function.

I use the following code to NDSolve it

rho0[x_, y_] := 1 - x - y
sigmoid[x_] := Tanh[x]
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}]

but it produces no solution:

enter image description here enter image description here

However, this simple instance actually has an analytical solution: $H(x,y)=F(x)$ and $V(x,y)=F(y)$ where $F(t) := (t - t^2)/2$. (This solution is very likely the unique one.) I'm wondering why the numerical solver gets nowhere near this theoretical solution.

The ultimate goal is to solve the system given more nested (but smooth at least) $\rho_0$'s. (An exogenous condition for $\rho_0$ is $\int_0^1 \int_0^1 \rho_0(x,y)dxdy=0$.) But I get stuck with the result above.

Any help appreciated! Thanks in advance.

Mathematica Version 12.0.0 student.

  • 1
    $\begingroup$ fyi, Tried it on Maple. It gives gives 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 Mathematica DSolve can not solve this analytically. Hopefully V 12.2 it will. $\endgroup$
    – Nasser
    Apr 27, 2020 at 3:48
  • $\begingroup$ full Maple code of the above if someone is interested 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$
    – Nasser
    Apr 27, 2020 at 3:49

1 Answer 1


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

NDSolveFEMInitializePDECoefficients::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

enter image description here

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.