Solution of the homeworks

1. How does the stability of the origin depend on the parameter p in the case of the system x'[t]=A x[t] with

A = {{2,p}, {-3,3}}

A = {{2,-5}, {-3,p}}

Bence Tipary

The first system:

ClearAll[A1,A2,A3,A4,A5]
A1[p_]:={{2,p},{-3,3}};

Let us determine the stationary point of the system

Solve[A1[p].{x,y}=={0,0},{x,y}]

We have only one stationary point (in the origin)

The eigenvalues:

Eigenvalues@A1[p]

The sign of the eigenvalues as a fubction of in the function of the parameter:

Make figures but do not believe what you are seeing!

Plot[{(1/2) (5-(Sqrt[1-12 p])),(1/2) (5+(Sqrt[1-12 p]))},{p,-6,(1/12)},AxesLabel->{"p",ev}]

Positive if:

Simplify[Reduce[Re@#>0,p]&/@Eigenvalues@A1[p],Assumptions->p∈Reals]

FullSimplify[%[[1]],Assumptions->p∈Reals]

Negative if:

Simplify[Reduce[Re@#<0,p]&/@Eigenvalues@A1[p],Assumptions->p∈Reals]

Is equal to zero if:

Simplify[Reduce[Re@#==0,p]&/@Eigenvalues@A1[p],Assumptions->p∈Reals]

The eigenvectors:

Eigenvectors@A1[p]

The streamplot around the origin:

Animate[StreamPlot[A1[p].{x,y},{x,-2,2},{y,-2,2},PlotLabel->TraditionalForm["p="~~ToString[p]~~"\neigenvalues="~~ToString[N@Eigenvalues[A1[p]]]],PerformanceGoal->"Quality"],{p,-5,5,0.01},AnimationRate->0.5,AnimationRunning->False,AnimationRepetitions->1]

We have a saddle to unstable node to unstable focus transition here.

Eigenvalues@A1[-3]//N

The second system:

A2[p_]:={{2,-5},{-3,p}};

The stationary point:

Solve[A2[p].{x,y}=={0,0},{x,y}]

The eigenvalues:

Eigenvalues@A2[p]

Plot[{(1/2) (2+p-(Sqrt[64-4 p+(p^2)])),(1/2) (2+p+(Sqrt[64-4 p+(p^2)]))},{p,-5,15}]

The sign of the eigenvalues in the function of the parameter:

Positive if:

Simplify[Reduce[Re@#>0,p]&/@Eigenvalues@A2[p],Assumptions->p∈Reals]

Negative if:

Simplify[Reduce[Re@#<0,p]&/@Eigenvalues@A2[p],Assumptions->p∈Reals]

Is equal to zero if:

Simplify[Reduce[Re@#==0,p]&/@Eigenvalues@A1[p],Assumptions->p∈Reals]

The eigenvectors:

Eigenvectors@A2[p]

The streamplot around the stationary point:

Animate[StreamPlot[A2[p].{x,y},{x,-2,2},{y,-2,2},PlotLabel->TraditionalForm["p="~~ToString[p]~~" eigenvalues="~~ToString[N@Eigenvalues[A2[p]]]],PerformanceGoal->"Quality"],{p,5,10,0.01},AnimationRate->0.5,AnimationRunning->False]

We have a saddle to unstable node transition here.

Tamás Mezei

ClearAll[A, StreamA, Ev, check, Stab];
A[p_]:=<|1->{{2, p}, {-3, 3}}, 2->{{2, -5}, {-3, p}}|>;

Ev[p_, number_]:=Eigenvectors[N@A[p][number]];

Ev[][]

check[p_, number_]:={AllTrue[Flatten[Im@Ev[p, number]], #==0&],AllTrue[Flatten[Re@Ev[p, number]], #>0&]};

AllTrue[{2,3,4,6,8}, EvenQ]

Stab[p_, number_]:=Which[
check[p, number]=={False,True}, "Unstable Focus", 
check[p, number]=={False,False}, "Stable Focus", 
check[p, number]=={True,False}, "Saddle"]

StreamA[p_, number_]:=Show[StreamPlot[A[p][number].{x, y}, {x, -2, 2}, {y, -2, 2}, 
PlotLabel->Grid[{{"Eigenvalues: ", Eigenvalues[N@A[p][number]], SpanFromLeft}, {"Eigenvectors: ", MatrixForm@First@Ev[p, number], MatrixForm@Last@Ev[p, number]}, {Stab[p, number], SpanFromLeft, SpanFromLeft}}],PlotRange->2, Mesh->Full],
If[AllTrue[Flatten[Im@Ev[p, number]], #==0&], Graphics[{Arrow[{{0, 0}, {Ev[p, number][[1, 1]], Ev[p, number][[1,2]]}}], Arrow[{{0, 0}, {Ev[p, number][[2,1]], Ev[p, number][[2, 2]]}}]}], Graphics[]]];

Manipulate[StreamA[p, number], 
{{number, 1, "Case"}, {1-> Grid[{{Text[Style["A =", Italic]], MatrixForm[{{2, Text[Style["p", Italic]]}, {-3, 3}}]}}], 2-> Grid[{{Text[Style["A =", Italic]], MatrixForm[{{2, -5}, {-3, Text[Style["p", Italic]]}}]}}]}}, {{p, 0}, -50, 50, Appearance->"Open"}]

Tamás Pölöskei

A1[p_]:=(({{2, p}, {-3, 3}}));

A2[p_]:=(({{2, -5}, {-3, p}}));

hw121[fun_]:=Module[{},
Manipulate[
GraphicsRow[
{Graphics[
{Thickness[0.01],
Arrow[{{0,0},{Re[Eigenvalues[fun[p]][[1]]],Im[Eigenvalues[fun[p]][[1]]]}}],
Arrow[{{0,0},{Re[Eigenvalues[fun[p]][[2]]],Im[Eigenvalues[fun[p]][[2]]]}}]
},
PlotRange->{{-6,6},{-6,6}},
Axes->True
],
StreamPlot[fun[p].{x,y},{x,-5,5},{y,-5,5},PlotLabel->Eigenvalues[N@fun[p]]]},
ImageSize->500
]
,{p,-3,3,0.5}
]
]

hw121[A1]

hw121[A2]

2. Let us draw the direction field for the systems below (use StreamPlot and VectorPlot)

Bence Tipary

The first system:

A3={{0,1},{-1,0}};

Stationary point:

Solve[A3.{x,y}=={0,0},{x,y}]

Stream and vectorplot around the stationary point:

StreamPlot[A3.{x,y},{x,-1,1},{y,-1,1}]

VectorPlot[A3.{x,y},{x,-.1,.1},{y,-.1,.1}]

We have a center here.

The second system:

A4={{-x,1},{-1,-y}};

The two stationary points are:

Solve[A4.{x,y}=={0,0},{x,y},Reals]

The streamplot:

StreamPlot[A4.{x,y},{x,-3,2},{y,-2,3}]

The vectorplots around the stationary points:

VectorPlot[A4.{x,y},{x,-.1,.1},{y,-.1,.1}]

In the small vicinity o {0, 0} there is a center.

VectorPlot[A4.{x,y},{x,-1.1,-.9},{y,.9,1.1}]

In the small vicinity of {-1, 1} there is a saddle.

The third system:

A5={{x,1},{-1,y}};

We have again two stationary points:

Solve[A5.{x,y}=={0,0},{x,y},Reals]

The streamplot:

StreamPlot[A5.{x,y},{x,-2,3},{y,-3,2}]

The vectorplots around the stationary points:

VectorPlot[A5.{x,y},{x,-.1,.1},{y,-.1,.1}]

In the small vicinity of {0, 0} there is a center.

VectorPlot[A5.{x,y},{x,.9,1.1},{y,-1.1,-.9}]

In the small vicinity of {1, -1} there is a saddle.

Tamás Pölöskei

StreamPlot[{y,-x},{x,-2,2},{y,-2,2},StreamScale->0.06,StreamStyle->Red]

StreamPlot[{y-(x^2),-x-(y^2)},{x,-2,2},{y,-2,2},StreamScale->0.06,StreamStyle->Red]

StreamPlot[{y+(x^2),-x+(y^2)},{x,-2,2},{y,-2,2},StreamScale->0.06,StreamStyle->Red]

3. Write a program to decide if a given function is the first integral of a system or not (for autonomous ordinary differential equations). Let us apply the program for the case of harmonic oscillator and total energy (potential and kinetic).

V[p_,q_]:=q+((p^2)/4)
f[x_,y_]:={-x,2y+(x^2)}

W[p_,q_,r_]:=r+((p^2)/3)
g[x_,y_,z_]:={-x,-y+(x^2),z+(x^2)}

This works to any number of variables.

FirstIntegralQ[fival_,rhsval_,vars_:{x,y}]:=D[fival,{vars,1}].rhsval

FirstIntegralQ[V[x,y],f[x,y]]

FirstIntegralQ[V[p,q],f[p,q]]//Trace

Erroneous calculation! The gradient is zero if one uses the wrong variables.

FirstIntegralQ[V[p,q],f[p,q],{p,q}]

FirstIntegralQ[W[p,q,r],g[p,q,r],{p,q,r}]

Tamás Pölöskei

How can I determine the number of variables of a function?

This will work for two variable functions.

hw123twovar[V_,f_]:=Module[{x,y,t},
Print[If[FullSimplify[D[V[x[t],y[t]],t]/.Thread[{x'[t],y'[t]}->f[x[t],y[t]]]]===0,"It is a first integral.","It isn't a first integral."]];
]

V1[x_,y_]:=y+((x^2)/4)
f1[x_,y_]:={-x,2y+(x^2)}

hw123twovar[V1,f1]

V2[x_,y_]:=(x^2)+(y^2)
f2[x_,y_]:={y,-x}

hw123twovar[V2,f2]

This will work for three variable functions.

hw123threevar[V_,f_]:=Module[{x,y,z,t},
Print[If[FullSimplify[D[V[x[t],y[t],z[t]],t]/.Thread[{x'[t],y'[t],z'[t]}->f[x[t],y[t],z[t]]]]===0,"It is a first integral.","It isn't a first integral."]];
]

W[x_,y_,z_]:=z+((x^2)/3)
g[x_,y_,z_]:={-x,-y+(x^2),z+(x^2)}

hw123threevar[W,g]

On the qualitative theory of ordinary differential equations 2.

Sufficient conditions of existence: repetition

The existence of first integrals with appropriate level sets

The trajectories go along the level curves of the first integral, if this is a closed curve, then the system has a periodic solution. See the harmonic oscillator or the Volterra Lotka system:

ClearAll[vlrhs,W];
vlrhs[x_,y_]={x-x y,x y-y};
W[p_,q_]:=Log[p]+Log[q]-p-q

D[W[p,q],{{p,q},1}].vlrhs[p,q]//Simplify

ContourPlot[W[p,q],{p,0,3},{q,0,3}]

Bendixson-sack

rhsb={x-y-(x^3),x+y-(y^3)};

Nullclines

p1=Plot[x-(x^3),{x,-2,2}]

p2=ParametricPlot[{-y+(y^3),y},{y,-2,2},PlotStyle->Red]

Show[p1,p2,PlotLabel->"The origin is the only stationary point."]

Show[p1,p2,PlotLabel->"The origin is the only stationary point.",Epilog->First@Graphics[{Thick,Brown,Circle[{0,0},1],Circle[{0,0},(Sqrt[2])]}],AspectRatio->Automatic]

StreamPlot[rhsb,{x,-2,2},{y,-2,2}]

Show[%%,%]

Inequalities, Poincaré map

*Construct the Poincaré map numericlally in this case. Use WhenEvent

Hopf-bifurcation

Manipulate[Module[{nds,x,y},nds=NDSolve[{x'[t]==1+a ((x[t])^2)y[t]-(b+1)x[t],y'[t]==a(b x[t]-a ((x[t])^2)y[t]),x[0]==8,y[0]==3},{x,y},{t,0,100}][[1]];
ParametricPlot3D[Evaluate[{{90,x[t],y[t]},{t,x[t],0},{t,4,y[t]}}/.nds],{t,90,100},BoxRatios->{1,1,1},PlotPoints->100,PlotRange->{{90,100},{0,4},{0,8}},AxesLabel->{Style["t",Italic],Row[{Style["x",Italic],"(",Style["t",Italic],")"}],Row[{Style["y",Italic],"(",Style["t",Italic],")"}]}, PlotStyle->Directive[Thick,RGBColor[.49,0,0]],ImageSize->{400,400},PlotLabel->"Trajectory and coordinate functions of the solutions\n\n"]],
{a,1.9,2.03,0.001,Appearance->"Labeled"},
{b,4.8,6.0,0.001,Appearance->"Labeled"}]

Necessary conditions of existence

Monotonicity

eq={x'[t]==-((x[t])^2),y'[t]==((x[t])^2)+2((y[t])^4)}

StreamPlot[(Last/@eq)/.{cc_[t]->cc},{x,-2,-3},{y,-4,5},PlotLabel->Style["Strictly monotonous components",Bold,18,Red]]

Bendixson´s theorem

x'=x-x y     y'=x y-y

The periodic solution intersects the bisectrix.

x'=x+y-((x^3)/3)     y'=-x (van der Pol)

The limit cycle must cross the vertical lines  x=±1.

f[x_,y_]:=x+y-((x^3)/3)
g[x_,y_]:=-x

div[Fval_]:=Tr[D[Fval,{{x,y},1}]]

div[{f[x,y],g[x,y]}]

x'=x-y-(x^3)     y'=x+y-(y^3)

div[{x-y-(x^3),x+y-(y^3)}]

Implies[1<(x^2)+(y^2)<2,2-3 (x^2)-3 (y^2)<0]//Simplify(*How do you prove this by hand?*)

RegionPlot[1<(x^2)+(y^2)<2,{x,-1.5,1.5},{y,-1.5,1.5},PlotLabel->"The div is negative,\nstill there is a periodic solutions in the annulus."]

Bendixson-Dulac´s theorem

Theorem. (Bendixson-Dulac) Suppose Ω⊂(ℝ^2) is a simply connected domain and let f∈(C^1)(Ω,(ℝ^2)); B∈(C^1)(Ω,ℝ) . If the divergence of B f does not change sign in Ω then the eq. has no closed orbit lying entirely in   Ω. If A⊂Ω is an annular region where the divergence does not change sign then there is at most one limit cycle.

Póta-Hanusse-Tyson-Light

x'=(y/(1+(x^2)))     y'=((-x+y(1+x+(x^2)))/(1+(x^2)))

An appropriate Dulac function is (x,y)↦(1+(x^2))

div[(1+(x^2)){(y/(1+(x^2))),((-x+y(1+x+(x^2)))/(1+(x^2)))}]

x'=y     y'=-a x-b y+α (x^2)+β (y^2)

B[x_,y_]:=b (E^(-2β x))

div[{y b (E^(-2β x)),(-a x-b y+α (x^2)+β (y^2))b (E^(-2β x))}]//Simplify

x'=x-x y     y'=x y-y

How do you choose a Dulac function to prove that closed trajectory intersects the hyperbola y=(1/x).

Theorem (Póta, Hanusse, Tyson, Light): Two-component bimolecular reactions cannot have a limit cycle.

An easy to read application of the Bendixson-Dulac theorem:

Póta, Gy. (1983), Two-component bimolecular systems cannot have limit cycles: A complete proof. J. Chem. Phys. 78:1621–1622

Homeworks

The system

eq={x'==x+(y^2)+(x^3),y'==-x+y+(x^2)y}

Has no periodic solution.

The system

eq={x'==x-x (y^2)+(y^3),y'==3y-(x^2)y+(x^3)}

Has no periodic solution in the disk with radius 2 centered at the origin.

Study the sensitivities of the parameters in the equation

x'[t]==-p x[t],y'[t]==p x[t]-q y[t],z'[t]==q y[t]

as a function of time at different parameter values.

Patterns

One can also use

RegularExpression

Defining functions

f1[x_]:=(x^2)

f1/@{1,2,3.4,π}

f2[x_Integer]:=(x^2)

f2/@{1,2,3.,π}

f3[x_]:=ppp[x]/;x>0

f3/@{2,-4}

The earlier, the better.

fac[n_/;n>0]:=n!

Condition

fac[6]+fac[-4]

fac[6.2]+fac[-4]

f4[x_]:=Module[{a},a=Prime[(10^x)];
(FactorInteger[a+1])/;a<(10^6)]

f4[4]

Prime[(10^4)]

%+1

f4[4]//Trace

f4[5]

Prime[(10^5)]

v[x_,1-x_]:=p[x]

v[(a^2),1-(a^2)]

v[4,-3]

v[4,1-4]

w[x_,y_]:=p[x]/;y==1-x

w[4,-3]

MemberQ[{x,(x^n)},n]

MemberQ[{x,(x^n)},n,Infinity]

FreeQ[{x,(x^n)},n]

ClearAll[x,h,a,b];h[a_ b_,x_]:=a h[b,x]/;FreeQ[a,x]

h[a b x,x]+h[2 (1+x) (x^2),x]

q[{x_Integer,y_Integer}?(Function[v,v.v>4])]:=qp[x+y]

{q[{3,4}],q[{1,1}],q[{-5,-7}]}

ClearAll[y,q];{1,x,(x^2),(x^3),(y^2)}/.(x|x^_)->q

ClearAll[y,q];{1,x,(x^2),(x^3),(y^2)}/.(x^_)->q

ClearAll[y,q];{1,x,(x^2),(x^3),(y^2)}/.(x^_.)->q

ClearAll[y,q];{1,x,(x^2),(x^3),(y^2)}/.x^_->q

ClearAll[y,q];{1,x,(x^2),(x^3),(y^2)}/.x^_.->q

{a[2],b[3],c[4],a[5]}/.(f:(a|b))[x_]->r[f,x]

Pattern sequences

Repeated patterns

ClearAll[f,a,b];
Cases[{f[a],f[a,a,b],f[a,b,a],f[a,b,b]},f[a..,b..]]

vn[x:{{_,n_}..}]:=Transpose[x]

vn[{{a,2},{b,2},{c,2}}]

vn[{{a,2},{b,2},{c,3}}]

vn[{{a,2,1},{b,2,1},{c,2,1}}]

vn[{{a,{2,1}},{b,{2,1}},{c,{2,1}}}]

vni[x:{{_,n_Integer}..}]:=Transpose[x]

f5[x_Symbol, p__Integer]:= Apply[Plus, x^{p}]

(x^({1,2}))

{f5[y,4],f5[uu,π],f5[x,1,2],f5[0,1,2],f5[x,1,2.]}

ClearAll[f,p,x];
f[x___Integer] := p[x,Plus[x]]

{f[],f[1],f[1,2,3],f[1,b,3]}

ClearAll[f];
f[x_Integer, y___] := Plus[y]^x

{f[1],f[1,2.],f[1,2.,x],f[x,1,y]}

Selections

Cases[{{1,2,3},a,{4,5},6},t:{__Integer}]

Cases[{{1,2,3},a,{4,5},6},t:{__Integer}:>(t^2)]

Full form:

Cases[{{1,2,3},a,{4,5},6},Pattern[t,{__Integer}]:>(t^2)]

Cases[{z[1,1+(Sqrt[Gamma[π/E]])],z[-1,1],z[-2,2]},z[x_/;x<0,y_]]

Cases[{a,b,0,1,2,x,y},Except[0,_Integer]]

Count[{3,4,x,x^2,x^3},x^_.]

DeleteCases[{3,4,x,2+x,3+x},_Integer,Infinity]

DeleteCases[{3,4,x,2+x,3+x},_Integer]

Cases: pattern, Select:criterion

Select

Select[{1,2,4,7,6,2},#>2&]

Select[{1,2,4,7,6,2},#>2&,1]

ClearAll[f,a,b];Select[f[1,a,2,b,3],IntegerQ]

ClearAll[f,a,b];Select[Plus[1,a,2,b,3],IntegerQ]

Tuples[{a,b},4]

Select[Tuples[{a,b},4],#==Reverse[#]&]

Select[Tuples[{0,1},{3,3}],Det[#]==1&,4]

Select[Names["*"],StringLength[#]<3&]

Head[7 (π^2) (x^2) (y^2)]

(7 (π^2) (x^2) (y^2))[[0]]

Select[7 (π^2) (x^2) (y^2),NumericQ]

Select[7 (π^2) (x^2) (y^2),NumberQ]

Local assignments/evaluation

f[a]+f[b]/.f[x_]->(x^2)

f[a]

f[{a,b}]+f[c]/.f[{x_,y_}]->p[x+y]

{1,x,(x^2),(x^3)}/.(x^(n_))->r[n]

{1,x,(x^2),(x^3)}/.(x^(n_.))->r[n]

{a/b,1/(b^2),2/(b^2)}/.(b^(n_))->d[n]

FullForm[{a/b,1/(b^2),2/(b^2)}]

{6,-7,3,2,-1,-2}/.x_/;x<0->w

{a,4,5,b}/.x_Integer->p[x]

{2.3,4,7/8,a,b}/.(x_/;NumberQ[x])->x^2

?*Q

Finding the position

Position[{f[a],g[b],f[c]},f[x_]]

Position[{4,4+x^a,x^b,6+x^5},x^_,Infinity,2](*first two powers ofxappearing at any level*)

FullForm[4+x^a]

TreeForm[4+x^a]

Homeworks

1. Consider the differential equation x´=(Subscript[a, 11])x+(Subscript[a, 12])y   y'=(Subscript[a, 21])x+(Subscript[a, 22])y. Construct a figure in the trace-det plane showing the type of the origin (node, saddle, center, focus). (Small phase plots may also be inserted.)

Filling, RegionPlot

2. Construct 12 matrices of the size 10×10 with random entries 0 or 1. Select those matrices which are adjacency matrices of connected graphs and make a figure of each of these graphs.

RandomInteger,Select, GraphPlot, Graph, AdjacencyGraph, AdjacencyMatrix,ConnectedGraphQ