# Problem with non-linear system of PDEs

I wish to solve the following system of hyperbolic PDEs: $$\frac{\partial y}{\partial u}=\alpha \sqrt{2y-y^2}$$, $$\frac{\partial y}{\partial v}=-\frac{1}{\alpha}\frac{\partial\alpha}{\partial v}\:\left(1-y\right)$$, for the functions $$\alpha=\alpha(u,v)$$ and $$y=y(u,v)$$ on a rectangle aligned with the axes on the first quadrant, and boundary conditions $$\alpha (u,0)=1$$, $$y(u,0)=\sin u+1$$, $$y(0,v)=1$$. However, I posed the problem as follows

eq1 = D[y[u, v], u] == \[Alpha][u, v] Sqrt[2 y[u, v] - y[u, v]^2];
eq2 = D[y[u, v], v] \[Alpha][u,v] == -D[\[Alpha][u, v], v] (1 - y[u, v]);
bc = {\[Alpha][u, 0] == 1, y[u, 0] == 1 + Sin[u], y[0, v] == 1};
region = Rectangle[{0, 0}, {1.5, 1.5}];
sol = NDSolve[{{eq1, eq2}, bc}, {\[Alpha], y}, {u, v} \[Element]region]



but I am not obtaining any solution, just a string of errors starting with

NDSolveFEMInitializePDECoefficients::femcnsd: The PDE coefficient -(((2-2 y[u,v]) \[Alpha][u,v])/(2 Sqrt[2 y[u,v]-y[u,v]^2])) does not evaluate to a numeric scalar at the coordinate {1.2075,1.2075}; it evaluated to Indeterminate instead.


and ending with

NDSolve::fempslf: The linearization process in PDESolve failed.


Anybody could please guide me about what is going wrong?

• Seems like a blow up of your PDE. Oct 31, 2021 at 21:44
• Please, pay attention that bc at u=0, v=0 are formulated outside of region. Also your region is not "a rectangle aligned with the axes on the first quadrant". Nov 1, 2021 at 3:44
• @AlexTrounev Thanks, I fixed it. The errors however persist. Nov 1, 2021 at 6:58
• Using the FEM is not a good idea for this problem, as it is convection dominated. The tensor product grid method is better in this case. Nov 2, 2021 at 9:54

This problem can be solve with using option "MethodOfLines". Note, that in a case of boundary condition $$\alpha (u,0)=1$$ the solution is trivial $$\alpha (u, v)=1$$. Nevertheless, NDSolve generates solution with some message.

eq1 = D[y[u, v], u] == \[Alpha][u, v] Sqrt[2 y[u, v] - y[u, v]^2];
eq2 = D[y[u, v], v] \[Alpha][u,
v] == -D[\[Alpha][u, v], v] (1 - y[u, v]);
bc = {\[Alpha][u, 0] == 1, y[u, 0] == 1 + Sin[u], y[0, v] == 1};

sol = NDSolve[{{eq1, eq2}, bc}, {\[Alpha], y}, {u, 0, 1.5}, {v, 0,
1.5}, Method -> {"IndexReduction" -> Automatic,
"EquationSimplification" -> "Residual",
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 41, "MaxPoints" -> 41,
"DifferenceOrder" -> 2}}}] // Quiet;


Visualization

{DensityPlot[\[Alpha][u, v] /. sol[[1]], {u, 0, 1.2}, {v, 0, 1.5},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotLabel -> "\[Alpha]", FrameLabel -> {"u", "v"}],
DensityPlot[y[u, v] /. sol[[1]], {u, 0, 1.2}, {v, 0, 1.5},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotLabel -> "y", FrameLabel -> {"u", "v"}]}