4
$\begingroup$

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 ?

Improve this question
$\endgroup$
10
  • 1
    $\begingroup$ There is at least one typo: You haven't squared the last term that appears on the left-hand-side of the equations. $\endgroup$
    – Diffycue
    Sep 21, 2022 at 21:11
  • $\begingroup$ 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 .. $\endgroup$
    – Nasser
    Sep 21, 2022 at 23:21
  • $\begingroup$ .. mathematica.stackexchange.com/questions/220477/… $\endgroup$
    – Nasser
    Sep 21, 2022 at 23:25
  • 1
    $\begingroup$ @Nasser - with v13.1 on my Mac I don't get an answer just warning messages. $\endgroup$
    – Bob Hanlon
    Sep 22, 2022 at 1:27
  • $\begingroup$ I am using windows. This is strange. I did not think the OS version would make difference. $\endgroup$
    – Nasser
    Sep 22, 2022 at 1:32

1 Answer 1

Reset to default
3
$\begingroup$

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] == 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]}

Figure 1

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 Figure 2
Third solution Figure 3 Fourth solution Figure 4

We can solve this system with Matematica FEM using special regularization as follows

ClearAll["Global`*"]
Needs["NDSolve`FEM`"]

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[[1]], {u, 0, 1}, {v, 0, 1}, Exclusions -> None,
   AxesLabel -> Automatic, PlotLabel -> "x", 
  ColorFunction -> "Rainbow", PlotTheme -> "Marketing", 
  MeshStyle -> White], 
 Plot3D[y[u, v] /. sol[[1]], {u, 0, 1}, {v, 0, 1}, Exclusions -> None,
   AxesLabel -> Automatic, PlotLabel -> "y", 
  ColorFunction -> "Rainbow", PlotTheme -> "Marketing", 
  MeshStyle -> White]}

Figure 5 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]];
Improve this answer
$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – Daniel Castro
    Sep 24, 2022 at 13:18
  • $\begingroup$ @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$. $\endgroup$
    – Alex Trounev
    Sep 24, 2022 at 14:22
  • $\begingroup$ Is there a way to impose those constraints in NDSolve ? $\endgroup$
    – Daniel Castro
    Sep 24, 2022 at 14:50
  • $\begingroup$ @DanielCastro Please, see update to my answer. $\endgroup$
    – Alex Trounev
    Sep 24, 2022 at 18:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.