ZWang

10 Reputation

0 Badges

0 years, 48 days

MaplePrimes Activity


These are questions asked by ZWang

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_end
    then evs(k,1..4):=Array(map(rhs, w));
         solA(eventclear);
    else break;
    fi
od:

interface(warnlevel=3):
M:=DeleteRow(convert(evs,matrix),1):

T:=M[..,1]:

for i from 1 do
tt:=T[i]:
if tt>=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;

HFloat(0.8000000001553702)

(1)

NULL

Rotation_number.mw

Thank you!

Very kind wishes,

Wang Zhe

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

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

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

Dear All, 

I work on the period doubling phenomenon by integrating the dynamical system at some particular parameter values. My main commands are as follows

restart:
with(DEtools): with(plots): with(plottools):
a:=-1: b:=-3: c:=3: d:=1: omega:=1: v1:=1: f:=-4: epsilon:=0.01: 
sys:=diff(u1(t),t)=v1*u1(t)-(omega+u2(t)^2)*u2(t)+u1(t)*(a*(u1(t)^2+u2(t)^2)+b*z(t)^2),
     diff(u2(t),t)=(omega+u1(t)^2)*u1(t)+v1*u2(t)+u2(t)*(a*(u1(t)^2+u2(t)^2)+b*z(t)^2),
     diff(z(t),t)=z(t)*(-v1+c*(u1(t)^2+u2(t)^2)+z(t)^2)+epsilon*z(t)*(v2+f*z(t)^4):

solC:=dsolve({eval(sys, v2=2.0500014987),
              u1(0)=0.6,u2(0)=0.6,z(0)=0.1},
              type=numeric, method=rkf45, maxfun=10^7,
              range=350..750):
p1:=odeplot(solC, [sqrt(u1(t)^2+u2(t)^2),z(t)],
               refine=2,
               color=burgundy,
               thickness=1);
 

The resulting curves are quite blurring and furry. I have tried to increase the "refine" or the "numpoints" values, but with minimal improvement. I would like to know if there a way to obtain clean and sharp integral curves?  

I would be very grateful if someone could help me! 

Thank you in advance.

With kind wishes,

Wang Zhe

1 2 Page 1 of 2