How to solve a matrix PDE and stop solving when solution becomes singular?

Question

My question consists of two parts:

1. How do I get mathematica to solve a PDE Matrix system and plot the result? See below for the PDE matrix system. (By plot the result I mean plot the region where the solution $\Theta$ is nonsingular.)
2. How do I stop the integration when the solution matrix becomes singular? I know that away from $(x_1,x_2)=(0,0)$ the solution matrix $\Theta$ will become singular how do I stop Mathematica from trying to solve pass this point?

The PDE matrix system I am trying to solve is \begin{align} \lambda \Theta(x_1,x_2)&=\dot{\Theta}(x_1,x_2)+\Theta(x_1,x_2)A &\quad \text{Equation}\\ \Theta(0,0)&=\begin{pmatrix} 1 & -\frac{1}{2} -\frac{\sqrt{3}}{2} \\ 1 & -\frac{1}{2}-\frac{1}{2\sqrt{3}}+\frac{2}{\sqrt{3}} \end{pmatrix} &\quad \text{Initial condition} \end{align} I have specified the values of $\dot{\Theta}(x_1,x_2),A,\lambda$ in the block below:

(* Definitions *)
A = {{-(x1^2 - 1), -2 x2 x1 - 1}, {1, 0}}
lambda = {{1/2, -Sqrt[3]/2}, {Sqrt[3]/2, 1/2}}
(*Value of theta for x1=x2=0 *)
ThetaInit = {{1, -1/2 - Sqrt[3]/2}, {1, -1/2 - 1/(2 Sqrt[3]) + 
    2/Sqrt[3]}}
(*Derative of Theta in terms of t. Note \
\frac{dtheta}{dt}=x1'(t)\frac{\partial theta}{\partial x1}+x2'(t) \
\frac{\partial theta}{\partial x2} *)
ThetaDot = ( -(x1^2 - 1) x2 - x1) D[Theta[x1, x2], x1] + 
  x2 D[Theta[x1, x2], x2]

Note the initial condition is satisfied as

lambda.ThetaInit == ThetaDot + ThetaInit.A /. {x1 -> 0, x2 -> 0}

yields true.

Notes

• If you need any clarification please feel free to ask.
• This problem comes from https://dspace.mit.edu/bitstream/handle/1721.1/9793/42916380-MIT.pdf?sequence=2 I am trying to reproduce the simulation in Mathematica see pages 64-66. I did not share this before because the details are fairly confusing.
• There was an error in the PDE I corrected that error.
• Because people where confused $$\dot{\Theta}(x_1,x_2)=( -(x_1^2 - 1) x_2 - x_1) \frac{\partial \Theta}{\partial x_1}+ x_2 \frac{\partial \Theta}{\partial x_2}$$

Answer 1

This question requests solutions for a first-order, linear, homogeneous partial differential equation. Standard procedures are available in the literature, which we follow here.

Characteristics

Solutions propagate along characteristics of the PDE. As suggested in comments embedded in the code in the question, characteristics satisfy {x1'[t] == -(x1[t]^2 - 1) x2[t] - x1[t], x2'[t] == x2[t]}, where t is a parameterization of the characteristics. (In this case, x2 also could be used as the parameter, but doing so is a bit less convenient.) The characteristic equations can be solved symbolically to yield,

FullSimplify[DSolveValue[{x1'[t] == -(x1[t]^2 - 1) x2[t] - x1[t], x2'[t] == x2[t]}, 
    {x1[t], x2[t]}, t], t \[Element] Reals]
(* {(C[1] (-2 BesselK[1, E^t Sqrt[C[1]^2]] + BesselI[1, E^t Sqrt[C[1]^2]] C[2]))/
    (Sqrt[C[1]^2] (2 BesselK[0, E^t Sqrt[C[1]^2]] + BesselI[0, E^t C[1]] C[2])), E^t C[1]} *)

There are two branches to this expression, the first of which reduces to

s1 = {(-2 BesselK[1, E^t C[1]] + BesselI[1, E^t C[1]] C[2])/
    ( 2 BesselK[0, E^t C[1]] + BesselI[0, E^t C[1]] C[2]), E^t C[1]};
p1 = ParametricPlot[Evaluate@Table[s1 /. {C[1] -> 1, C[2] -> 10^i}, {i, -3, 5}], {t, -4, 4}, 
    PlotRange -> {-3, 3}, AspectRatio -> 1, AxesLabel -> {x1, x2},
    ImageSize -> Large, LabelStyle -> {Black, Bold, 15}];

and the second to

s2 = {(2 BesselK[1, E^t C[1]] + BesselI[1, E^t C[1]] C[2])/
      (-2 BesselK[0, E^t C[1]] + BesselI[0, E^t C[1]] C[2]), E^t C[1]};
p2 = ParametricPlot[Evaluate@Table[s2 /. {C[1] -> 1, C[2] -> 10^i}, {i, -3, 5}], {t, -4, 4}, 
    PlotRange -> {-3, 3}, AspectRatio -> 1, AxesLabel -> {x1, x2},
    ImageSize -> Large, LabelStyle -> {Black, Bold, 15}];

Solutions are symmetric about the origin, as is evident from the characteristics equations. Therefore, solutions with x2 < 0 are obtained by replacing s1 and s2 by their negatives in the code above. The four plots together provide the total set of characteristics.

Show[p1, p2, p3, p4]

Several separatrices are evident, and {0, 0} is an equilibrium point (also as noted in the report referenced in the question). The characteristics can, of course, also be computed by numerical integration of the characteristic equations. For instance,

NDSolveValue[{x1'[t] == -(x1[t]^2 - 1) x2[t] - x1[t], x2'[t] == x2[t],
    x1[0] == 3, x2[0] == -.12}, {x1[t], x2[t]}, {t, 0, 5}];
tst = ParametricPlot[%, {t, 0, 5}, PlotRange -> {{-3, 3}, {-3, 3}}, 
    AspectRatio -> 1, AxesLabel -> {x1, x2}, PlotStyle -> Black, 
    ImageSize -> Large, LabelStyle -> {Black, Bold, 15}]

yields a curve corresponding almost precisely the the Red curve mostly in the lower right quadrant of the plot above.

Solutions along Characteristics

Solutions of the PDE can be obtained from lambda.theta == D[theta, t] + theta.a, integrated along characteristics. (thetaDot in the question can be obtained from D[theta, t] by the chain rule, if desired.) If theta is represented as

theta = {{theta11[t], theta12[t]}, {theta21[t], theta22[t]}}

then the equations needed by NDSolve are obtained by

eqth = Flatten[Thread /@ Thread[lambda.theta == D[theta, t] + theta.a]];
inth = Flatten[Thread /@ Thread[(theta /. t -> 0) == thetaInit]]

Expressions for the characteristics obtained above then can be inserted into a, and the equations integrated. However, it is simpler to integrate the characteristic equations and the theta equations simultaneously. For instance,

vars = Flatten[{{x1[t], x2[t]}, theta}];

tmax = 5.7; 
sol = NDSolveValue[{{x1'[t] == -(x1[t]^2 - 1) x2[t] - x1[t], x2'[t] == x2[t], 
    x1[0] == 0.01, x2[0] == .01}, eqth, inth}, vars, {t, 0, tmax}];
ParametricPlot[sol[[1 ;; 2]], {t, 0, tmax}, PlotRange -> {{-3, 3}, {-3, 3}}, 
    AspectRatio -> 1, AxesLabel -> {x1, x2}, PlotStyle -> Black, 
    ImageSize -> Large, LabelStyle -> {Black, Bold, 15}]
Plot[Evaluate@sol[[3 ;; 6]], {t, 0, tmax}, PlotRange -> All, 
    PlotLegends -> {theta11, theta12, theta21, theta22}, 
    AxesLabel -> {t, th}, ImageSize -> Large, LabelStyle -> {Black, Bold, 15}]

There is no sign of a singularity in the solution, either here or for any other parameters I have examined. Indeed, one would not expect singularities to arise from this PDE, because it is linear in theta and its coefficients are well behaved.