# Fixed points of coupled NLD involving transcendental functions

Posted 1 month ago
676 Views
|
12 Replies
|
9 Total Likes
|
 I have the Thomas system at hand and want to find out the fixed points of the equation using Mathematica for the whole range of the parameter. The system of equations are $$x'[t]=-b*x[t]+sin[y[t]]$$ $$y'[t]=-b*y[t]+sin[z[t]]$$ $$z'[t]=-b*z[t]+sin[x[t]]$$ where b ranges from 0 to 1. I have tried the following mathematica code,  eqns={x'[t]=-b*x+Sin[y],y'[t]=-b*y+Sin[z],z'[t]=-b*z+Sin[x]}; soln=Solve[eqns==0,{x,y,z}] For this I got an error "This system cannot be solved with the methods available to Solve". How can I solve the above system of equations ? How can I get a full picture of the fixed points?
12 Replies
Sort By:
Posted 1 month ago
 This gives some complex solutions: With[{b = 1, r = 4}, Solve[{0 == Sin[y] - b x, 0 == Sin[z] - b y, 0 == Sin[x] - b z, Abs[x] < r, Abs[y] < r && Abs[z] < r}, {x, y, z}]] % // N 
Posted 1 month ago
 With a plot like this you can make guesses about possible nontrivial real solutions: With[{b = 1, r = 2 Pi}, ContourPlot3D[{0 == Sin[y] - b x, 0 == Sin[z] - b y, 0 == Sin[x] - b z}, {x, -r, r}, {y, -r, r}, {z, -r, r}]] 
Posted 1 month ago
 After running your code I got a warning like this' Solve::incs: Warning: Solve was unable to prove that the solution set found is complete.'. Now the Lyapunov spectrum shows that for b=1 the solutions are such that they converge to a fixed point. But here we have many solutions including the complex ones.
Posted 1 month ago
 Also, I got for fractional values b the following warning ' Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result' what is going on in this situation.
Posted 1 month ago
 Ok, I understood the reason for this error. I can't give fractional values with Solve for example if I want to put 0.5 I need to put it like 1/2.
Posted 1 month ago
 Hello Vijay,I am afraid I don't know what you mean by "fixed points of the equation", but perhaps this could be of interest for you Manipulate[ lsg = NDSolve[ {x'[t] == -b*x[t] + Sin[y[t]], y'[t] == -b*y[t] + Sin[z[t]], z'[t] == -b*z[t] + Sin[x[t]], x[0] == x0, y[0] == y0, z[0] == z0}, {x, y, z}, {t, 0, 50} ] // Flatten; fx = x /. lsg; fy = y /. lsg; fz = z /. lsg; Plot[{fx[t], fy[t], fz[t]}, {t, 0, 50}, PlotRange -> {-4, 4}], {{b, 1/2}, 0, 1}, {{x0, 0}, -2, 2}, {{y0, 0}, -2, 2}, {{z0, 0}, -2, 2}] 
Posted 1 month ago
 Hii, Thank you for the reply..... Well what I meant by 'fixed points' is that I want to find out those values of x ,y, z, say $x^*, y^* ,z^*$ for which x'[t]=0,y'[t]=0 and z'[t]=0.. Of course, your plots give the nature of dynamics for different values b. I am interested in seeing how many fixed points are there and also how they evolve for different values b.
Posted 1 month ago
 Vijay Arjun,the main part of my approach was already done above by the approach of @Hans Dolhaine: This nice presentation of the dynamics implies that a fixed point (always?) has the form {x,x,x}, i.e. all coordinates have the same value. Assuming this to be true then the equations reduce to $$x'(t) = -b *x(t) +\sin(x(t))$$At a fixed point we have $x'(t)=0$, which gives a simple relation between the fixed point coordinate $x$ and the parameter $b$:$$b=\frac{\sin(x)}{x}$$Next I calculate the eigenvalues of the respective Jakobian and put the above relation into: ClearAll["Global*"]; jackobian = {{-b, Cos[y], 0}, {0, -b, Cos[z]}, {Cos[x], 0, -b}}; (* eigenvalues of Jakobian: *) ev = Refine[Re@Eigenvalues[jackobian /. {b -> Sinc[x], y -> x, z -> x}], x \[Element] Reals]; isneg = Boole[And @@ (Sign[#] < 0 & /@ ev)]; Plot[{ev, -isneg, Sinc[x]}, {x, 0, 23}, PlotLabels -> {"ev1", "ev2", "ev3", None, "b"}, PlotStyle -> {Blue, Blue, Blue, Green, Red}, ImageSize -> 700, Filling -> {4 -> Axis}] Here $b$ and the real parts of the eigenvalues (depending on $x$) are shown; the green bars are indicating the intervals where all these reals parts are negative - which is required for a fixed point to be stable/attractive.And finally these intervals can be calculated: (* Product of eigenvalues: *) evProd = Times @@ ev; (* reasonable starting points: *) xstart = {2.1272, 7.73540, 8.04482, 14.1171, 14.3105, 20.34412, 20.49883}; (* intervals for possible fixed points: *) xIntervals = BlockMap[Interval, Prepend[x /. FindRoot[evProd == 0, {x, #}] & /@ xstart, 0], 2]; (* respective intervals for 'b' *) bIntervals = Sinc[xIntervals]; GraphicsColumn[{NumberLinePlot[xIntervals, GridLines -> Automatic, PlotLabel -> "Fixed point intervals"], NumberLinePlot[bIntervals, GridLines -> Automatic, PlotLabel -> "Island of stability in 'b'"]}, ImageSize -> 700] I hope this is helpful - but in any case you should check my reasoning carefully!Regards -- Henrik
Posted 1 month ago
 @ Henrik: Coool.
Posted 1 month ago
 Hii... Thank you for the reply... But I think you wrote the equation with x. We need y. There is a cyclic symmetry among the equations...I am talking about the sine term. My understanding of the equation is that there is a fixed point between b=0.329(appro) and b=1. At the former value, oscillations will begin showing a bifurcation. When we decrease the value of b, the system will execute the chaotic motion. This we can see clearly from plots of @Hans Dolhaine:
 But I think you wrote the equation with x. We need y. Well, if you read my answer you find that I concentrated on fixed points where all coordinate values are identical - you may call it as you like, 'x' or 'y' or whatever. At the former value [b=0.329(appro)], oscillations will begin showing a bifurcation. When we decrease the value of b, the system will execute the chaotic motion. This we can see clearly from plots of @Hans Dolhaine No, if you look carefully you just see "islands of stability" - as I call it. I explicitly calculated the intervals for stable fixed points ('xIntervals') and the respective intervals for 'b' ('bIntervals'). E.g. a possible fixed point candidate is {8.,8.,8.,}, and then 'b' has to be Sinc[8.] (which gives a value as low as 0.12367). Just try it using the code of @Hans Dolhaine : {x0, y0, z0} = {8., 8., 8.}; b = Sinc[8.]; lsg = NDSolve[{x'[t] == -b*x[t] + Sin[y[t]], y'[t] == -b*y[t] + Sin[z[t]], z'[t] == -b*z[t] + Sin[x[t]], x[0] == x0, y[0] == y0, z[0] == z0}, {x, y, z}, {t, 0, 2000}] // Flatten; fx = x /. lsg; fy = y /. lsg; fz = z /. lsg; Plot[{fx[t], fy[t], fz[t]}, {t, 0, 2000}, PlotRange -> {0, 10}] You can go further using the next interval, e.g. {x0, y0, z0} = {14.1, 14.1, 14.1}; b = Sinc[14.1] (* Out: 0.07087 *) and again you will find a stable system.
 For a given vlue of b you can use Solve (or NSolve) provided the domain is suitably restricted. Solve[{x == 2 Sin[y], y == 2 Sin[z], z == 2 Sin[x], 0 <= x <= 2 Pi, 0 <= y <= 2 Pi, 0 <= z <= 2 Pi}, {x, y, z}] (* Out[582]= {{x -> 0, y -> 0, z -> 0}, {x -> 2 Sin[2 Sin[ Root[{-2 Sin[2 Sin[2 Sin[#1]]] + #1 &, 1.89549426703398094714}]]], y -> 2 Sin[ Root[{-2 Sin[2 Sin[2 Sin[#1]]] + #1 &, 1.89549426703398094714}]], z -> Root[{-2 Sin[2 Sin[2 Sin[#1]]] + #1 &, 1.89549426703398094714}]}} *) `