# Problem with NDSolve for 1st order non-linear system of PDEs

I wish to solve the PDE system

\begin{align} x_u^2+y_u^2+\left(xx_u+yy_u \right)^2&=1\\ x_v^2+y_v^2+\left(xx_v+yy_v \right)^2&=1 \end{align}

subjected to $$x(0,v)=0$$ and $$y(u,0)=0$$. However, when using

eq1 = D[x[u, v], u]^2 + D[y[u, v], u]^2 + (x[u, v]*D[x[u, v], u] + y[u, v]*D[y[u, v], u])^2 == 1
eq2 = D[x[u, v], v]^2 + D[y[u, v], v]^2 + (x[u, v]*D[x[u, v], v] + y[u, v]*D[y[u, v], v])^2 == 1
bc = {x[0, v] == 0, y[u, 0] == 0};
sol = NDSolve[{{eq1, eq2}, bc}, {x, y}, {u, 0, 1}, {v, 0, 1}];


I get instead of a solution just a bunch of warnings, including NDSolveFEM InitializePDECoefficients and FindRoot::stfail: The method AffineCovariantNewton failed to compute the next step. Any ideas about how to solve it ?

• There is at least one typo: You haven't squared the last term that appears on the left-hand-side of the equations. Sep 21, 2022 at 21:11
• I get instead of a solution just a bunch of warnings this is strange. Using V 13.1 it does give solutions, but they are both zero when plotted. I see the warnings, but you say you get no solution? Screen shot (fixed your ode to add ^2) !Mathematica graphics so for the warnings issue, then need to play with options to get try to get rid of them. There are few questions on this site with similar warnings and some workaround for them to try, Here is one such example .. Sep 21, 2022 at 23:21
• Sep 21, 2022 at 23:25
• @Nasser - with v13.1 on my Mac I don't get an answer just warning messages. Sep 22, 2022 at 1:27
• I am using windows. This is strange. I did not think the OS version would make difference. Sep 22, 2022 at 1:32

There is exact solution to this problem. First, let put $$x=x(u), y=y(v)$$, then we can use DSolve as follows

DSolve[{D[x[u], u]^2+(x[u]*D[x[u], u])^2==1, x == 0}, x, u]
(*{{x -> Function[{u},
InverseFunction[
1/2 (ArcTanh[#1/Sqrt[1 + #1^2]] + #1 Sqrt[
1 + #1^2]) &][-u]]}, {x ->
Function[{u},
InverseFunction[
1/2 (ArcTanh[#1/Sqrt[1 + #1^2]] + #1 Sqrt[1 + #1^2]) &][u]]}}*)


Similar solution we have for $$y(v)$$. We can plot solutions as follows

lst = Table[{u,
Evaluate[
InverseFunction[
1/2 (ArcTanh[#1/Sqrt[1 + #1^2]] + #1 Sqrt[1 + #1^2]) &][
u]]}, {u, 0, 1, .01}];

f = Interpolation[lst];
{Plot3D[f[u], {u, 0, 1}, {v, 0, 1}, Exclusions -> None,
AxesLabel -> Automatic, PlotLabel -> "x",
ColorFunction -> "Rainbow", PlotTheme -> "Marketing",
MeshStyle -> White],
Plot3D[f[v], {u, 0, 1}, {v, 0, 1}, Exclusions -> None,
AxesLabel -> Automatic, PlotLabel -> "y",
ColorFunction -> "Rainbow", PlotTheme -> "Marketing",
MeshStyle -> White]} Please, pay attention that there are 4 solutions in combination of signs x[u], y[v]. It is why numerical solution not unique, and therefore we have a problem with computation. Second solution Third solution Fourth solution We can solve this system with Matematica FEM using special regularization as follows

ClearAll["Global*"]
Needs["NDSolveFEM"]

reg = Rectangle[{0, 0}, {1, 1}]; mesh =
ToElementMesh[reg, MaxCellMeasure -> .001]

eq1 = Inactivate[
d1 Laplacian[x[u, v], {u, v}] + D[x[u, v], u]^2 +
D[y[u, v],
u]^2 + (x[u, v]*D[x[u, v], u] + y[u, v]*D[y[u, v], u])^2,
D | Laplacian];
eq2 = Inactivate[
d2 Laplacian[y[u, v], {u, v}] + D[x[u, v], v]^2 +
D[y[u, v],
v]^2 + (x[u, v]*D[x[u, v], v] + y[u, v]*D[y[u, v], v])^2,
D | Laplacian];
bc = {DirichletCondition[x[u, v] == 0, u == 0],
DirichletCondition[y[u, v] == 0, v == 0]};

sol = NDSolve[{Activate[eq1] == 1, Activate[eq2] == 1,
bc} /. {d1 -> 10^-3, d2 -> 10^-3}, {x, y},
Element[{u, v}, mesh]];


Visualization

{Plot3D[x[u, v] /. sol[], {u, 0, 1}, {v, 0, 1}, Exclusions -> None,
AxesLabel -> Automatic, PlotLabel -> "x",
ColorFunction -> "Rainbow", PlotTheme -> "Marketing",
MeshStyle -> White],
Plot3D[y[u, v] /. sol[], {u, 0, 1}, {v, 0, 1}, Exclusions -> None,
AxesLabel -> Automatic, PlotLabel -> "y",
ColorFunction -> "Rainbow", PlotTheme -> "Marketing",
MeshStyle -> White]} Compare to analytical solution this is second case shown in Figure 2. To compute all cases, we change a sign of d1,d2, for instance, first and fourth solution are given by

sol1 = NDSolve[{Activate[eq1] == 1, Activate[eq2] == 1,
bc} /. {d1 -> -10^-3, d2 -> -10^-3}, {x, y},
Element[{u, v}, mesh]];
sol4 = NDSolve[{Activate[eq1] == 1, Activate[eq2] == 1,
bc} /. {d1 -> 10^-3, d2 -> -10^-3}, {x, y},
Element[{u, v}, mesh]];

• Thank you. I knew the exact solution, but I am trying to understand the numerical problem to be able to use it in more complicated cases. I thus chose this one which would allow me to compare the exact and numerical results. Sep 24, 2022 at 13:18
• @DanielCastro The problem stated above in my answer. There is no unique solution, therefore any numerical method for nonlinear PDEs based on the Newton iterative method will compute several roots. As result solution does not converge. We can try numerical optimization with constrains like $x(u,v)\ge 0, y(u,v)\ge 0$. Sep 24, 2022 at 14:22
• Is there a way to impose those constraints in NDSolve` ? Sep 24, 2022 at 14:50
• @DanielCastro Please, see update to my answer. Sep 24, 2022 at 18:05