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.
{x1, x2}
do you wish a solution, and what are the boundary conditions for this region? It seems unlikely that specifyingTheta[0, 0]
is sufficient. $\endgroup$ – bbgodfrey Dec 3 '19 at 6:24{0, 0}
. Please clarify. $\endgroup$ – bbgodfrey Dec 3 '19 at 6:36lambda.ThetaInit == (ThetaDot + A) /. {x1 -> 0, x2 -> 0}
yieldsFalse
. In other words,ThetaInit
does not satisfy the matrix PDE at{0, 0}
. $\endgroup$ – bbgodfrey Dec 3 '19 at 15:11