## 0 Badges

0 years, 107 days

## How to plot a prettier bifurcation diagr...

Maple

Dear All,

The following code plots the bifurcation diagram for a three-dimensional continuous dynamical system as a variable Re varies. However, the resulting plot (by pointplot command) is rather ugly, comparing with other bifurcation diagrams, see attached. Could anyone point me out how to improve its looking?

 > restart: with(plots): with(DEtools): with(plottools):with(LinearAlgebra): doSol:=proc( R )              local t_start:=2500,                    t_end:=5000,                    b:=-3,                    c:=3,                    v1:=1,                    f:=-4,                    v2:=2.0,                    omega:=0.1*sqrt(R),                    epsilon:=evalf(1/R),                    k:=0.1,                    kH:=(c+1)/(b-1),                    sys, evs, w, M, T, i, tt, solA, N_alpha,                    r_center, z_center, Zshift, alpha, alpha1,                    alpha2, del_alpha, m, Z, Rshift, RR;              sys:=diff(u1(t),t)=v1*u1(t)-(omega+k*u2(t)^2)*u2(t)-(u1(t)^2+u2(t)^2-b*z(t)^2)*u1(t),                   diff(u2(t),t)=(omega+k*u1(t)^2)*u1(t)+v1*u2(t)-(u1(t)^2+u2(t)^2-b*z(t)^2)*u2(t),                   diff(z(t),t)=z(t)*(kH*v1+c*u1(t)^2+c*u2(t)^2+z(t)^2)+epsilon*z(t)*(v2+f*z(t)^4):              solA:=dsolve({sys, u1(0)=0.6, u2(0)=0.6, z(0)=0.1},                           {u1(t),u2(t),z(t)},                           numeric, method=rkf45, maxfun=10^7,                           events=[[[u1(t)=0, u2(t)>0], halt]]                          );              evs:=Array():              evs(1,1..4):=Array([t,u1(t),u2(t),z(t)]);              interface(warnlevel=0):              for i from 2 do                  w:=solA(t_end):                  if   rhs(w[1])=t_start                  then break;                  end if;              end do:              RR:=M[i..,3]: Z:=M[i..,4]:              N_alpha:=numelems(RR):              r_center:=add(RR[i],i=1..N_alpha)/N_alpha:              z_center:=add(Z[i],i=1..N_alpha)/N_alpha:              Zshift:=Z-~r_center: Rshift:=RR-~z_center:              alpha:=(arctan~(Zshift/~Rshift)+(Pi/2)*sign~(Rshift))/~(2*Pi)+~0.5:              return alpha;   end proc:
 > i_start:=125: i_end:=250: i_delta:=0.1: M:=[seq( [j, doSol(j)], j=i_start..i_end, i_delta)]:
 > S:=seq([Vector(numelems(M[i,2]),fill=M[i,1]),M[i,2]], i=1..numelems(M),1): display(seq( pointplot(S[i],              symbolsize=1,              symbol=point) ,              i=1..numelems(M)),              axes=boxed,              view=[i_start..i_end,0..1],              size=[650,400],              axesfont= ["TimesNewRoman", 16],              labels=["Re",phi[n]],              labelfont = ["TimesNewRoman", 16],              labeldirections=[horizontal, vertical]);

Bifurcation_diagram.mw

Thank you.

Very kind wishes,

Wang Zhe

## How to clear variables after each iterat...

Maple

Dear All,

The following code calculates the Rotation_number for given set of parameters. I would like to plot the value of Rotation_number against R for R from 500 to 1000. Since there are many variables for each iteration, I would like to know if there is a way to create a for loop for the following code while automatically clear variables except for Rotation_number to speed it up?

 > restart: with(plots): with(DEtools): with(plottools):with(LinearAlgebra): t_start:=2500: t_end:=5000: b:=-3: c:=3: v1:=1: f:=-4: v2:=2.0: omega:=1: epsilon:=evalf(1/R): R:=100: k:=0.25: kH:=(c+1)/(b-1): sys:=diff(u1(t),t)=v1*u1(t)-(omega+k*u2(t)^2)*u2(t)-(u1(t)^2+u2(t)^2-b*z(t)^2)*u1(t),      diff(u2(t),t)=(omega+k*u1(t)^2)*u1(t)+v1*u2(t)-(u1(t)^2+u2(t)^2-b*z(t)^2)*u2(t),      diff(z(t),t)=z(t)*(kH*v1+c*u1(t)^2+c*u2(t)^2+z(t)^2)+epsilon*z(t)*(v2+f*z(t)^4): solA:=dsolve({sys, u1(0)=0.6, u2(0)=0.6, z(0)=0.1},              {u1(t),u2(t),z(t)},               numeric, method=rkf45, maxfun=10^7,               events=[[[u1(t)=0, u2(t)>0], halt]]): evs:=Array(): evs(1,1..4):=Array([t,u1(t),u2(t),z(t)]): interface(warnlevel=0): for k from 2 do     w:=solA(t_end):     if rhs(w[1])=t_start then break; end if; end do: R:=M[i..,3]: Z:=M[i..,4]: N_alpha:=numelems(R): r_center:=add(R[i],i=1..N_alpha)/N_alpha: z_center:=add(Z[i],i=1..N_alpha)/N_alpha: Zshift:=Z-~r_center: Rshift:=R-~z_center: alpha:=(arctan~(Zshift/~Rshift)+(Pi/2)*sign~(Rshift))/~(2*Pi)+~0.5: alpha2:=alpha[2..N_alpha]: alpha1:=alpha[1..N_alpha-1]: del_alpha:=alpha2-alpha1: m:=numelems(del_alpha): for i from 1 to m do if del_alpha[i]<0 then del_alpha[i]:=del_alpha[i]+1; else del_alpha[i]:=del_alpha[i]; fi: od: Rotation:=add(del_alpha[i],i=1..m)/m;
 (1)

Rotation_number.mw

Thank you!

Very kind wishes,

Wang Zhe

## Plot3d for data imported from an excel f...

Maple

Dear All,

The excel file consists of data (v2, Re, t) and I would like plot the variable "t" aginst the variables "v2" and "Re". Could anyone point me out?

Attached is the excel data.Collapsetime.xlsx

Thank you.

Very kind wishes,

Wang Zhe

## How to probability density function of a...

Maple

Dear All,

I would like to plot the probability density function of a state variable obtained from solving differential equations. I have found that there are functions called "PDF" and "KernelDensityPlot" in the Statistics package, but they really confuse me. Could you please point me out? My code is as follows.

Ps. Is it possible to plot the PDF directly from the solution of dsolve() without discretizing the results?

restart:
with(plots): with(DEtools): with(plottools):with(LinearAlgebra): with(Statistics):

v1:=1: f:=-4: v2:=2.515: omega:=1: epsilon:=0.001: k:=0:
sys:=diff(u1(t),t)=v1*u1(t)-(omega+k*u2(t)^2)*u2(t)-(u1(t)^2+u2(t)^2+3*z(t)^2)*u1(t),
diff(u2(t),t)=(omega+k*u1(t)^2)*u1(t)+v1*u2(t)-(u1(t)^2+u2(t)^2+3*z(t)^2)*u2(t),
diff(z(t),t)=z(t)*(-v1+3*u1(t)^2+3*u2(t)^2+z(t)^2)+epsilon*z(t)*(v2+f*z(t)^4):

t_start:=50: t_end:=300: dt:=0.05: fs:=1/dt:

solA:=dsolve({sys, u1(0)=0.6, u2(0)=0.6, z(0)=0.1},
{u1(t),u2(t),z(t)},
type=numeric, method=rkf45, maxfun=0,
output=Array([seq(i,i=t_start..t_end, dt)])):

u1:=solA[2,1][..,2]:
u2:=solA[2,1][..,3]:
z:=solA[2,1][..,4]:

u0:=sqrt~(u1^~2+u2^~2+z^~2):

Phi:=z/~u0:

# How could I plot the probability density of Phi (y-axis) against Phi(x-axis)?

Probability_density_function_plot.mw

Thank you!

Very kind wishes,

Wang Zhe

## Frequency analysis for dsolve()...

Maple 2016

Dear All,

I am working on ODEs and have obtained the plot for "variable vs time". I would like to know if it is possible and how to analyze those data in the frequency domain.ODEs.mw

Thank you.

Very kind wishes,

Wang Zhe

 1 2 Page 1 of 2
﻿