3
$\begingroup$

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}$$
||improve this question
$\endgroup$
  • 1
    $\begingroup$ How do you calculate $x_1'(t)\frac{\partial \theta}{\partial x_1}+x_2'(t) \ \frac{\partial \theta}{\partial x_2}$? $\endgroup$ – xzczd Dec 3 '19 at 6:00
  • 2
    $\begingroup$ Over what region in {x1, x2}do you wish a solution, and what are the boundary conditions for this region? It seems unlikely that specifying Theta[0, 0] is sufficient. $\endgroup$ – bbgodfrey Dec 3 '19 at 6:24
  • 1
    $\begingroup$ Unless, that is, you wish only solutions lying along a characteristic curve passing through {0, 0}. Please clarify. $\endgroup$ – bbgodfrey Dec 3 '19 at 6:36
  • 1
    $\begingroup$ lambda.ThetaInit == (ThetaDot + A) /. {x1 -> 0, x2 -> 0} yields False. In other words, ThetaInit does not satisfy the matrix PDE at {0, 0}. $\endgroup$ – bbgodfrey Dec 3 '19 at 15:11
  • $\begingroup$ @xzczd $x'_1(t),x'_2(t)$ are given in my post. $\frac{\partial \Theta}{\partial x_1}$ is a matrix. $\Theta$ is what needs to be solved. $\endgroup$ – AzJ Dec 3 '19 at 16:03
3
+300
$\begingroup$

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]

enter image description here

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

enter image description here

enter image description here

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.

||improve this answer
$\endgroup$
  • $\begingroup$ Thanks for your answer it is very comprehensive. I do have a few questions. 1) what is the difference between the p1,p2,p3,p4 plots? Different initial conditions or something else? I want to reproduce the plot in your answer to make sure I understand them, but I cannot reproduce show[p1,p2,p3,p4] if I don't know the difference between different p. $\endgroup$ – AzJ Dec 7 '19 at 19:55
  • $\begingroup$ 2) If you change tmax to a sufficiently large value (in my testing tmax \approx 20) $\Theta_{12},\Theta_{22}$ diverges to infinity. From pg 65 in the text: "The nonlinear eigenvector field Θ can be computed recursively by solving numerically [PDE system] along system trajectories. The result is illustrated in Figure... Also note that in the plot we have saturated Θ (which becomes infinite) at the boundary of the eigenvector field." My question is can I tell mathematica to stop integration when Θ becomes saturated? $\endgroup$ – AzJ Dec 7 '19 at 20:08
  • 1
    $\begingroup$ @AzJ p3 is obtained by replacing s1 by -s1, and p4 by replacing s2 by -s2. I shall add s2 to the answer for clarity. I too can obtain solutions for theta that diverge to infinity, but only when x1 or x2 also diverge to infinity. I presumed that you were not seeking such cases. If you have found a solution that diverges for finite x1 and x2. please send me the initial conditions for x1 and x2, so that I can take a look. Thanks for thoroughly examining my answer. $\endgroup$ – bbgodfrey Dec 7 '19 at 21:11
  • 1
    $\begingroup$ @AzJ You can obtain thetaDot from D[theta, t] by the chain rule, namely D[theta, t] = D[theta, x1] x1'[t] + D[theta, x2] x2'[t]. With respect to your second comment, NDSolve stops automatically when it encounters a singularity. To stop it sooner, use WhenEvent. Run NDSolveValue[{{x1'[t] == -(x1[t]^2 - 1) x2[t] - x1[t], x2'[t] == x2[t], x1[0] == -1.8, x2[0] == 1}, eqth, inth}, vars, {t, -5, 7}] to see NDSolve stop by itself. Note that this case corresponds to x1` approaching infinity. $\endgroup$ – bbgodfrey Dec 8 '19 at 6:56
  • 1
    $\begingroup$ @AzJ Thank you very much. It is an interesting problem. $\endgroup$ – bbgodfrey Dec 8 '19 at 19:10

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.