(* oscillateur de Lorenz *)

eqs={x'[t] ==- sigma x[t] + sigma y[t], y'[t]== r x[t] - y[t] - z[t] x[t], z'[t]== - b z[t] + x[t] y[t]}

NumVal ={sigma-> 10, b-> 8/3, r ->28}

inCond = {x[0]== 2, y[0]== 1, z[0]== 2}

xx = {x[t], y[t], z[t]}

sol=NDSolve[Join[eqs,inCond]/.NumVal,xx, {t, 0,50},AccuracyGoal->2, PrecisionGoal->10] [[1]]

oscgraph=ParametricPlot3D[xx/. sol, {t,0,50},PlotStyle->{RGBColor[0,0,0],Thickness[0.002]}, PlotPoints->500]

(* Etude des points d'equilibre *)

equil=Solve[Thread[Part[#,2]& /@ eqs == {0,0,0}],xx]

equilN =equil/.NumVal

eqs /. NumVal

(* linearisation autour des points d equilibre *)

ff=Part[#,2]& /@ eqs

AA =Outer[D,ff,xx]

AAs=AA /.equil

AAs[[1]]//MatrixForm

AAs[[2]]//MatrixForm

AAs[[3]] //MatrixForm

eigval=Eigenvalues[#]&/@ (AAs/.NumVal //N) //Chop

eigvect=Eigenvectors[#]&/@ (AAs/.NumVal //N) //Chop

(* calcul des plans des variété centre, i.e. ou l'oscillation a lieu pour le systeme linearise *)

(* equilibre 2 *)

equilN[[2]]

xx/.equilN[[2]]//N

eigvect[[2]]

(* on examine les vecteurs complexes et on en deduit par combinaison lineaire deux vecteurs independant reels *)

gen21=eigvect[[2]][[2]]+eigvect[[2]][[3]]//Chop

gen31=eigvect[[3]][[2]]+eigvect[[3]][[3]]//Chop

gen22 =I(eigvect[[2]][[2]]-eigvect[[2]][[3]])//Chop

gen32 =I(eigvect[[3]][[2]]-eigvect[[3]][[3]])//Chop

(* Normale au plan *)

n1=Cross[gen21,gen22]

n2  = Cross[gen31,gen32]

param2=Solve[n1.xx == 0,z[t]][[1]]

param3=Solve[n2.xx == 0,z[t]][[1]]

xxbar2=xx/.equilN[[2]]//N

xxbar3=xx/.equilN[[3]]//N

plan2 = ParametricPlot3D[(xx/.param2)+xxbar2, {x[t],-3,3},{y[t],-3,3},PlotStyle->Opacity[0.8],Mesh->None];

plan3 = ParametricPlot3D[(xx/.param3)+xxbar3, {x[t],-3,3},{y[t],-3,3}, PlotStyle->Opacity[0.8],Mesh->None];

Show[oscgraph,plan2, plan3]

Show[plan2, plan3,oscgraph]

(* integration de le connection heteroclinique *)

equil

eqsBack=Thread[D[xx,t]==(-Part[#,2]&/@ eqs)]

inCondRev2 = Thread[(xx/.{t->0})== 0.0001 eigvect[[2]][[1]]+ xxbar2]

inCondRev3 = Thread[(xx/.{t->0})== 0.0001 eigvect[[3]][[1]]+ xxbar3]

inCondRev2b = Thread[(xx/.{t->0})== -0.0001 eigvect[[2]][[1]]+ xxbar2]

inCondRev3b = Thread[(xx/.{t->0})== -0.0001 eigvect[[3]][[1]]+ xxbar3]

arc1Sol=NDSolve[Join[eqsBack,inCondRev2]/.NumVal,xx,{t,0,1.2}][[1]];

arc2Sol=NDSolve[Join[eqsBack,inCondRev3]/.NumVal,xx,{t,0,1.2}][[1]];

arc1Solb=NDSolve[Join[eqsBack,inCondRev2b]/.NumVal,xx,{t,0,1.2}][[1]];

arc2Solb=NDSolve[Join[eqsBack,inCondRev3b]/.NumVal,xx,{t,0,1.2}][[1]];

arc1=ParametricPlot3D[xx/.arc1Sol,{t,0,1.2}, PlotStyle->Thickness[0.01]];

arc1b= ParametricPlot3D[xx/.arc1Solb,{t,0,1.2}, PlotStyle->Thickness[0.01]];

arc2=ParametricPlot3D[xx/.arc2Sol,{t,0,1.2}, PlotStyle->Thickness[0.01]];

arc2b=ParametricPlot3D[xx/.arc2Solb,{t,0,1.2}, PlotStyle->Thickness[0.01]];

xx/.arc1Sol/.{t->1.2}

xx/.arc1Solb/.{t->1.2}

xx/.arc2Sol/.{t->1.2}

xx/.arc2Solb/.{t->1.2}

Show[oscgraph,plan2, plan3,arc1,arc1b,arc2,arc2b]