3
$\begingroup$

I am trying to solve the 1st order PDE system

\begin{align} \xi_u^2+\eta_u^2&= \left(1+\frac{\xi^2+\eta^2}{4} \right)^2\\ \xi_v^2+\eta_v^2&=\left(1+\frac{\xi^2+\eta^2}{4} \right)^2 \end{align}

for $\xi=\xi(u,v)$ and $\eta=\eta(u,v)$, with the conditions

\begin{align} &\xi(u,0)=2\tan\frac{u}{2}\:,\:\:&\eta(u,0)=0 \\ &\xi(0,v)=0\:,\:\:&\eta(0,v)=2\tan\frac{v}{2} \end{align}

When trying the following

  reg = Rectangle[{0, 0}, {1, 1}]; mesh = ToElementMesh[reg, MaxCellMeasure -> 0.00001]; 
  eq1 = Derivative[1, 0][\[Eta]][u, v]^2 + Derivative[1, 0][\[Xi]][u, v]^2 == (1 + (\[Eta][u, v]^2 + \[Xi][u, v]^2)/4)^2; 
  eq2 = Derivative[0, 1][\[Eta]][u, v]^2 + Derivative[0, 1][\[Xi]][u, v]^2 == (1 + (\[Eta][u, v]^2 + \[Xi][u, v]^2)/4)^2; 
  bc = {DirichletCondition[\[Xi][u, v] == 2*Tan[u/2], v == 0], DirichletCondition[\[Eta][u, v] == 0, v == 0], DirichletCondition[\[Xi][u, v] == 0, u == 0], DirichletCondition[\[Eta][u, v] == 2*Tan[v/2], u == 0]}; 
  sol = NDSolveValue[{eq1, eq2, bc}, {\[Xi], \[Eta]}, Element[{u, v}, mesh]]

I get just zero and the warnings

  FindRoot::stfail: The method AffineCovariantNewton failed to compute the next step.
  FindRoot::sszero: The step size in the search has become less than the tolerance prescribed by the PrecisionGoal option, but the function value is still greater than the tolerance prescribed by the AccuracyGoal option.

Adding some diffusion helps to get a solution with

  eq1 = (-\[Epsilon])*Laplacian[\[Xi][u, v], {u, v}] +Derivative[1, 0][\[Eta]][u, v]^2 + Derivative[1, 0][\[Xi]][u, v]^2 == (1 + (\[Eta][u, v]^2 + \[Xi][u, v]^2)/4)^2; 
  eq2 =  (-\[Epsilon])*Laplacian[\[Eta][u, v], {u, v}] + Derivative[0, 1][\[Eta]][u, v]^2 + Derivative[0, 1][\[Xi]][u, v]^2 == (1 + (\[Eta][u, v]^2 + \[Xi][u, v]^2)/4)^2; 
  bc = {DirichletCondition[\[Xi][u, v] == 2*Tan[u/2], v == 0], DirichletCondition[\[Eta][u, v] == 0, v == 0], DirichletCondition[\[Xi][u, v] == 0, u == 0], DirichletCondition[\[Eta][u, v] == 2*Tan[v/2], u == 0]}; 
  sol = sol = NDSolveValue[{eq1, eq2,bc} /. {\[Epsilon] -> 1/10}, {\[Xi], \[Eta]}, Element[{u, v}, mesh]]

The problem is that I need to get this $\epsilon$ much smaller to really approach the solution, but when trying for instance $1/50$ I get no answer, with

FindRoot::dfmin: The minimal damping factor of 1/10000 has been reached.
NDSolveValue::fempsf: PDESolve could not find a solution.
   
$\endgroup$
4
  • 2
    $\begingroup$ For extremely nonlinear problems you can try to restart the solver with an initial seed from a previous step. See for example here. But the finite element method is not well suited for convection dominant problems, but it's worth a shot. $\endgroup$
    – user21
    Oct 30, 2023 at 17:24
  • 1
    $\begingroup$ You may be able to reduce the complexity of the two PDEs by using independent variables {u + v, u - v} although at the cost of complicating the boundary conditions. It also may be useful to consider the result of DSolve[{eq1, eq2}, {\[Eta][u, v], \[Xi][u, v]}, {u, v}], although it is not a general solution. $\endgroup$
    – bbgodfrey
    Oct 31, 2023 at 4:47
  • $\begingroup$ @bbgodfrey Thank you. The first suggestion do not really help, because it even introduces more terms in the rotated equation. The second is indeed a particular solution that cannot adjust the boundary conditions. $\endgroup$ Oct 31, 2023 at 13:18
  • $\begingroup$ @bbgodfrey In the Eikonal equation the function itself does not appear, just its gradient and the independent variable, which is not the case here. Besides the type of equation, my concern is that Mathematica solves the equation for big values of the parameter $\epsilon$ but not for smaller ones. From the point of view of PDEs the solution exists for any real $\epsilon$, so I can't figure out what goes wrong for $\epsilon<<1.$ $\endgroup$ Oct 31, 2023 at 23:52

1 Answer 1

4
$\begingroup$

Update: Solution accuracy modestly improved and accuracy estimate appended.

The Finite Element Method used in the last section of the question approximates the PDEs by thousands of sparsely coupled nonlinear algebraic equations. The inability of FindRoot to solve these equations, indicated by the error messages near the end of the question, suggests that NDSolve is passing poor initial guesses for the solution to FindRoot. The InitialSeeding option eliminates this difficulty, even for very small artificial viscosity. Replace the last line of code in the question by

sol = NDSolve[{eq1, eq2, bc} /. {\[Epsilon] -> 10^-6}, {\[Xi], \[Eta]}, 
   Element[{u, v}, mesh], 
   InitialSeeding -> {\[Xi][u, v] == 2*Tan[u/2], \[Eta][u, v] == 2*Tan[v/2]}];

(Note that a mesh resolution of 10^-4, used here, provides good accuracy and is much faster than the mesh resolution of 10^-5 specified in the question.) Plots of the solution are given by

Plot3D[sol[[1, 2]][u, v], {u, 0, 1}, {v, 0, 1}, PlotLabel -> \[Xi],
    LabelStyle -> {12, Bold, Black}, AxesLabel -> {u, v}]
Plot3D[sol[[2, 2]][u, v], {u, 0, 1}, {v, 0, 1}, PlotLabel -> \[Eta],
    LabelStyle -> {12, Bold, Black}, AxesLabel -> {u, v}]

enter image description here

enter image description here

Daniel Castro in a comment below requested an estimate of the error introduced by non-zero \[Epsilon]. Such an estimate can be provided by substituting the numerical solutions obtained into the original PDEs, as follows.

Subtract @@@ {eq1, eq2} /. {\[Epsilon] -> 0} /. sol;
FunctionInterpolation[%[[1]], {u, 0, 1}, {v, 0, 1}];
Plot3D[%[u, v], {u, 0, 1}, {v, 0, 1}, PlotLabel -> \[Xi]err, 
    LabelStyle -> {12, Bold, Black}, AxesLabel -> {u, v}, PlotRange -> All]
FunctionInterpolation[%%%[[2]], {u, 0, 1}, {v, 0, 1}];
Plot3D[%[u, v], {u, 0, 1}, {v, 0, 1}, PlotLabel -> \[Eta]err, 
    LabelStyle -> {12, Bold, Black}, AxesLabel -> {u, v}, PlotRange -> All]

enter image description here

enter image description here

It is important to remember that these error estimates include errors both in the solution presented above and in the subsequent calculation of the error itself.

$\endgroup$
2
  • $\begingroup$ Thank you. This is more a mathematical analysis question, but do you have an idea about the convergence of the $\epsilon \rightarrow 0$ limit of the solution? I mean, one should be able to formally prove that adding artificial diffusion does not takes the solution away from the original $\epsilon=0$ problem. $\endgroup$ Nov 1, 2023 at 10:22
  • $\begingroup$ Your comment prompted me to reassess my answer and make some improvements. Even so, the original result was quite accurate. $\endgroup$
    – bbgodfrey
    Nov 1, 2023 at 15:06

Your Answer

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

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