0
$\begingroup$

I'm trying to get the Bifurcation diagrams for u, v, w, s, f, g and h as function of r

I got some kind of code to do this, but I do n't have any idea of what is going on with code and i'm getting bunch of errors

the problem is like:

Parameters

rf = 20000;r0 = 10.^-7;m = 2;k = 2;rE = 1000;a = 0.001;
a1 = (m (2 + m Sqrt[2 - n] - n + 
 E^((2 m)/Sqrt[2 - n]) (-2 + m Sqrt[2 - n] + n)) )/((1 + E^((2 m)/
 Sqrt[2 - n])) (2 - n)^(3/2));
\[Tau] = 1/(\[Pi]^2 + k^2);
r = Ta;

The system of the algebraic equations

F1 = k^2 \[Tau] Ta v[r] - u[r]/\[Tau] + (n u[r] )/\[Tau] + 
16/(3 \[Pi]) Ta k^2 \[Tau] v[r] s[r];
F2 = -(\[Pi]^2/n)  s[r] ;
F3 = -(\[Pi]^2/2) u[r] v[r] - (4 \[Pi]^2)/n w[r];
F4 = \[Pi] u[r] w[r] + u[r] - v[r]/(n \[Tau]) ;
F5 = rE/2 g[r] + 
16/(6 \[Pi]) rE g[r] s[r] - (4 n)/a f[r] + (4 n^2)/a f[r];
F6 = - (rE/2) f[r] + 
16/(6 \[Pi]) rE f[r] s[r] - (4 n)/a g[r] + (4 n^2)/a g[r];
F7 = - ((4 n)/a) h[r] + (4 n^2)/a h[r] - a1 rE u[r];

getting the maximum of each row

G1 = Table[{D[F1, r]}, {n, 0.2, 0.6, 0.1}]
Dimensions[G1]
G2 = Table[{D[F2, r]}, {n, 0.2, 0.6, 0.1}];
G3 = Table[{D[F3, r]}, {n, 0.2, 0.6, 0.1}];
G4 = Table[{D[F4, r]}, {n, 0.2, 0.6, 0.1}];
G5 = Table[{D[F5, r]}, {n, 0.2, 0.6, 0.1}];
G6 = Table[{D[F6, r]}, {n, 0.2, 0.6, 0.1}];
G7 = Table[{D[F7, r]}, {n, 0.2, 0.6, 0.1}];

Flatten

Difeq = Flatten[{Table[{Simplify[G1[[i, 1]]] == 0, 
  Simplify[G2[[i, 1]]] == 0, Simplify[G3[[i, 1]]] == 0, 
  Simplify[G4[[i, 1]]] == 0, Simplify[G5[[i, 1]]] == 0, 
  Simplify[G6[[i, 1]]] == 0, Simplify[G7[[i, 1]]] == 0}, {i, 1, 5,
   1}]}, 1];
Dimensions[Difeq]

soln = Flatten[{Table[{Solve[{Difeq[[i, 1]], Difeq[[i, 2]], 
    Difeq[[i, 3]], Difeq[[i, 4]], Difeq[[i, 5]], Difeq[[i, 6]], 
    Difeq[[i, 7]]}, {Derivative[1][u][r], Derivative[1][v][r], 
    Derivative[1][w][r], Derivative[1][s][r], Derivative[1][f][r],
     Derivative[1][g][r], Derivative[1][h][r]}]}, {i, 1, 5, 1}]}, 
   3];
Dimensions[soln]
equation = 
Flatten[{Table[{eq[i, 1] = 
   Simplify[Derivative[1][u][r] /. soln[[i, 1]]], 
  eq[i, 2] = Simplify[Derivative[1][v][r] /. soln[[i, 2]]], 
  eq[i, 3] = Simplify[Derivative[1][w][r] /. soln[[i, 3]]],
  eq[i, 4] = Simplify[Derivative[1][s][r] /. soln[[i, 4]]],
  eq[i, 5] = Simplify[Derivative[1][f][r] /. soln[[i, 5]]],
  eq[i, 6] = Simplify[Derivative[1][g][r] /. soln[[i, 6]]],
  eq[i, 7] = Simplify[Derivative[1][h][r] /. soln[[i, 7]]]}, {i, 
  1, 5, 1}]}, 1];
Dimensions[equation]
sol = Flatten[{Table[{NDSolve[{Derivative[1][u][r] == 
     Simplify[Derivative[1][u][r] /. soln[[i, 1]]], 
    Derivative[1][v][r] == 
     Simplify[Derivative[1][v][r] /. soln[[i, 2]]], 
    Derivative[1][w][r] == 
     Simplify[Derivative[1][w][r] /. soln[[i, 3]]], 
    Derivative[1][s][r] == 
     Simplify[Derivative[1][s][r] /. soln[[i, 4]]], 
    Derivative[1][f][r] == 
     Simplify[Derivative[1][f][r] /. soln[[i, 5]]], 
    Derivative[1][g][r] == 
     Simplify[Derivative[1][g][r] /. soln[[i, 6]]], 
    Derivative[1][h][r] == 
     Simplify[Derivative[1][h][r] /. soln[[i, 7]]], 
    u[r0] == 0.0001, v[r0] == 0.0001, w[r0] == 0.0001, 
    s[r0] == 0.0001, f[r0] == 0.0001, g[r0] == 0.0001, 
    h[r0] == 0.0001}, {u[r], v[r], w[r], s[r], f[r], g[r], 
    h[r]}, {r, r0, rf}, 
   Method -> {"Extrapolation", 
     Method -> "LinearlyImplicitEuler"}]}, {i, 1, 5, 1}]}, 3];
Dimensions[sol]

Plotting

q1 = Plot[Evaluate[ u[Ta] /. sol], {Ta, r0, rf}, 
AxesLabel -> {r,  Subscript[U, s][ r ]}, PlotLabel -> "(A)" , 
Frame -> True, 
FrameLabel -> {Style["r", Bold, Blue], 
Style["\!\(\*SubscriptBox[\(U\), \(s\)]\)[ r ]", Bold, Blue]}, 
PlotStyle -> {Black, Blue, Red, Thick, Dashed}]
$\endgroup$
2
  • 1
    $\begingroup$ There are a few minor issues to fix first. It's a good idea not to use capitals for your function and variable names. Your def of Difeq runs a table up to i=7, but your Gi only have 5 elements. Your def of soln, equation and sol have Tables indexed by {i, 1, 1, 1}, which is "i runs from 1 to 1 in steps of 1". Is this a placeholder for something else, or can you remove the Table entirely? The curly brackets outside your {Table[...]}s are superfluous. You don't appear to need the levelspec in Flatten, eg: in your def of Difeq, Flatten[stuff, 7] === Flatten[stuff]. $\endgroup$ Commented Aug 8, 2017 at 0:19
  • $\begingroup$ @aardvark2012 thank u for ur comment. I have edited the code. fortunately, I don't get any errors now, but the plot that I receive is blank $\endgroup$ Commented Aug 8, 2017 at 10:57

1 Answer 1

3
$\begingroup$

A couple of changes can get the plot working. First, some expressions are undefined at r = 0 (try equation /. Ta -> 0). But that can easily be fixed by setting r0 = 10.^-7. Second, since you set r = Ta you need to replace r with Ta in your plot function:

q1 = Plot[Evaluate[u[Ta] /. sol], {Ta, r0, rf}, 
  AxesLabel -> {r, Subscript[U, s][r]}, PlotLabel -> "(A)", 
  Frame -> True, 
  FrameLabel -> {Style["r", Bold, Blue], 
    Style["\!\(\*SubscriptBox[\(U\), \(s\)]\)[ r ]", Bold, Blue]}, 
  PlotStyle -> {Black, Blue, Red, Thick, Dashed}]

which gives

enter image description here

Hope that helps.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.