Question: How to clear variables after each iteration -- for loop?

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

Please Wait...