ijuptilk

85 Reputation

2 Badges

0 years, 329 days

MaplePrimes Activity


These are replies submitted by ijuptilk

@tomleslie 

Please how do I set colors in listcontplot? I wanna know the values that correspond to each color in the contour plot. 
 

  restart:
  doCalc:= proc( xi , u)  #u is the \bar(H): normalize magnetic field magnitude,
                          # where H = bar(H)*H__a
                 # Import Packages
                 uses ArrayTools, Student:-Calculus1, LinearAlgebra,
                      ListTools, RootFinding, plots, ListTools:
                 local gamma__1:= .1093,
                       alpha__3:= -0.1104e-2,
                       k__1:= 6*10^(-12),
                       d:= 0.2e-3,
                       theta0:= 0.1e-3,
                       eta__1:= 0.240e-1, chi:= 1.219*10^(-6),
                       alpha:= 1-alpha__3^2/(gamma__1*eta__1),
                       theta_init:= theta0*sin(Pi*z/d),
                       H__a:= Pi*sqrt(k__1/chi)/d,
                       H := u*H__a,     
                       c:=alpha__3*xi*alpha/(eta__1*(4*k__1*q^2/d^2-alpha__3*xi/eta__1 - chi*H^2)),
                       w := chi*H^2*eta__1*alpha/(4*k__1*q^2/d^2-alpha__3*xi/eta__1 - chi*H^2),
                       n:= 20,
                       g, f, b1, b2, qstar, OddAsymptotes, ModifiedOddAsym,
                       qstarTemporary, indexOfqstar2, qstar2, AreThereComplexRoots,
                       soln1, soln2, qcomplex1, qcomplex2, gg, qq, m, pp, j, i,
                       AllAsymptotes, p, Efun, b, aa, F, A, B, Ainv, r, theta_sol, v, Vfun, v_sol,minp,nstar;

# Assign g for q and plot g, Set q as a complex and Compute the Special Asymptotes

                 g:= q-(1-alpha)*tan(q)- (w*q + c)*tan(q):
                 f:= subs(q = x+I*y, g):
                 b1:= evalc(Re(f)) = 0:
                 b2:= evalc(Im(f)) = 0:
                 qstar:= (fsolve(1/c = 0, q = 0 .. infinity)):
                 OddAsymptotes:= Vector[row]([seq(evalf(1/2*(2*j + 1)*Pi), j = 0 .. n)]);

# Compute Odd asymptote

                 ModifiedOddAsym:= abs(`-`~(OddAsymptotes, qstar));
                 qstarTemporary:= min(ModifiedOddAsym);
                 indexOfqstar2:= SearchAll(qstarTemporary, ModifiedOddAsym);
                 qstar2:= OddAsymptotes(indexOfqstar2);

# Compute Odd asymptote

                 AreThereComplexRoots:= type(true, 'truefalse');
                 try
                      soln1:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = 0 .. infinity});
                      soln2:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = -infinity .. 0});
                      qcomplex1:= subs(soln1, x+I*y);
                      qcomplex2:= subs(soln2, x+I*y);
                 catch:
                       AreThereComplexRoots:= type(FAIL, 'truefalse');
                 end try;

# Compute the rest of the Roots

                 OddAsymptotes:= Vector[row]([seq(evalf((1/2)*(2*j+1)*Pi), j = 0 .. n)]);
                 AllAsymptotes:= sort(Vector[row]([OddAsymptotes, qstar]));
                 if   AreThereComplexRoots
                 then gg:= [ qcomplex1,
                             qcomplex2,
                             op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])),
                             seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)
                          ];
                 elif not AreThereComplexRoots
                 then gg:= [ op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])),
                             seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)
                           ];
                 end if:

# Remove the repeated roots if any & Redefine n

                 qq:= MakeUnique(gg):
                 m:= numelems(qq):

## Compute the `τ_n`time constants

                 for i to m do
                 p[i] := gamma__1*alpha/(4*k__1*qq[i]^2/d^2 - alpha__3*xi/eta__1- chi*H^2):
                 end do:

## Compute θ_n from initial conditions

                for i to m do
                Efun[i] := cos(qq[i]-2*qq[i]*z/d)-cos(qq[i]):
                end do:

## Compute integral coefficients for off-diagonal elements θ_n matrix

                printlevel := 2:
                for i to m do
                    for j from i+1 to m do
                        b[i, j] := int(Efun[i]*Efun[j], z = 0 .. d):
                        b[j, i] := b[i, j]:
                        aa[i, j] := b[i, j]:
                        aa[j, i] := b[j, i]:
                    end do :
                end do:

## Calculate integral coefficients for diagonal elements in theta_n matrix

                for i to m do
                aa[i, i] := int(Efun[i]^2, z = 0 .. d):
                end do:

## Calculate integrals for RHS vector

               for i to m do
               F[i] := int(theta_init*Efun[i], z = 0 .. d):
               end do:

## Create matrix A and RHS vector B

               A := Matrix([seq([seq(aa[i, j], i = 1 .. m)], j = 1 .. m)]):
               B := Vector([seq(F[i], i = 1 .. m)]):

## Calculate inverse of A and solve A*x=B

              Ainv := 1/A:
              r := MatrixVectorMultiply(Ainv, B):

## Define Theta(z,t)
             theta_sol := add(r[i]*Efun[i]*exp(-t/p[i]), i = 1 .. m):

## Compute v_n for n times constant

             for i to m do
             v[i] := (-2*k__1*alpha__3*qq[i])*(1/(d^2*eta__1*alpha*gamma__1))+ alpha__3^2*xi/(2*(eta__1)^2*qq[i]*alpha*gamma__1)+xi/(2*eta__1*qq[i]) + alpha__3*chi*H^2/(2*eta__1*qq[i]*gamma__1*alpha):
             end do:

## Compute v(z,t) from initial conditions
             for i to m do
             Vfun[i] := d*sin(qq[i]-2*qq[i]*z/d)+(2*z-d)*sin(qq[i]):
             end do:

## Define v(z,t)
             v_sol := add(v[i]*r[i]*Vfun[i]*exp(-t/p[i]), i = 1 .. m):

##
             minp:=min( seq( Re(p[i]), i=1..m) ):
             member(min(seq( Re(p[i]), i=1..m)),[seq( Re(p[i]), i=1..m)],'nstar'):

## Return all the plots
                 return theta_init, theta_sol, v_sol, minp, eval(p), nstar, p[nstar], g, H, H__a;
                 end proc:

# Run the calculation for supplied value of 'xi'

# Put the returned quantities in a simple list

ans:=[doCalc(-0.06, 5)]:
ans1:=[doCalc(0.06, 5)]:
evalf(ans[9]);
evalf(ans[10]);

174.2463358

 

34.84926715

(1)

##

with(plots):
d:= 0.2e-3:

## director Plot negative activity
display( plot(ans[1], z=0..d, color = green),
         plot([seq(eval(ans[2], t = j), j in [seq(.1*j, j = 1 .. 5)])], z=0..d, legendstyle = [font                 = ["ROMAN", 12],location = right], labeldirections = ["horizontal", "vertical"]),                     axis[1]=[tickmarks=[0="0", 0.00005="50", 0.00010="100", 0.00015="150", 0.00020=                         "200"]], axis[2]=[tickmarks=[0="0", 0.00002="2e-05", 0.00004="4e-05", 0.00006=                     "6e-05", 0.00008= "8e-05", 0.0001="1e-04"]], axes=boxed, labels =[typeset(z," (`μm`)                     "), typeset(theta(z, t)," (degrees)")], labelfont = ["TimesNewRoman", 14],                             axesfont = [14, 14]);


## director angle Plot postive activity
display( plot(ans1[1], z=0..d, color = green, legend = theta[0]),
         plot([seq(eval(ans1[2], t = j), j in [seq(.1*j, j = 1 .. 5)])], z=0..d, legend = [t[1] =                  .2, t[2] = .4, t[3] = .6, t[4] = .8, t[5] = 1.0], legendstyle = [font = ["ROMAN", 12],                location = right], labeldirections = ["horizontal", "vertical"]), axis[1]=[tickmarks=                [0="0", 0.00005="50", 0.00010="100", 0.00015="150", 0.00020="200"]], axis[2]=                        [tickmarks=[0="0", 0.00002="2e-05", 0.00004="4e-05", 0.00006="6e-05", 0.00008=                       "8e-05", 0.0001="1e-04"]], axes=boxed, labels =[typeset(z," (`μm`)"), typeset(theta(z, t)," (degrees)")], labelfont =["TimesNewRoman", 14], axesfont = [14, 14]);

## v Plot negative activity
plt2:= plot([seq(eval(ans[3], t = j), j in [seq(.2*j, j = 1 .. 5)])], z = 0 .. d, legendstyle = [font = ["ROMAN", 12],                location = right], labeldirections = ["horizontal", "vertical"]):

display( plt2, axis[1]=[tickmarks=[0="0", 0.00005="50", 0.00010="100", 0.00015="150", 0.00020="200"]          ], axes=boxed, labels =[typeset(z," (`μm`)"), typeset(v(z, t)," (m/s)")], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14]);
 
## v Plot postive activity
plt21:= plot([seq(eval(ans1[3], t = j), j in [seq(.2*j, j = 1 .. 5)])], z = 0 .. d, legend = [t[1] =                .2, t[2] = .4, t[3] = .6, t[4] = .8, t[5] = 1.0], legendstyle = [font = ["ROMAN", 12],                location = right], labeldirections = ["horizontal", "vertical"]):

display( plt21,axis[1]=[tickmarks=[0="0", 0.00005="50", 0.00010="100", 0.00015="150", 0.00020="200"]          ], axes=boxed, labels =[typeset(z," (`μm`)"), typeset(v(z, t)," (m/s)")], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14]);

  

 

 

 

 

# q plot

#

plot(ans[8], q=0..5*Pi,view=[DEFAULT,-10..10], labels =[q, r(q)], labelfont = ["TimesNewRoman", 14], labeldirections = ["horizontal", "vertical"] );

 

# p values

evalf(ans[5]):

 # Director angle plot for at d/2 as a funtion of time

theta_f := subs(z=d/2, Re(ans[2])):
plot(theta_f, t=0..10, view=[DEFAULT, DEFAULT], labels =[typeset(t," (seconds)"), typeset(theta('d/2', t),"(rad)")], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14], labeldirections = ["horizontal", "vertical"], axes=boxed);

theta_f := subs(z=d/2, ans1[2]):
plot(theta_f, t=0..1, ``view=[DEFAULT,DEFAULT], labels =[typeset(t," (seconds)"), typeset(theta('d/2', t),"(rad)")], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14], labeldirections = ["horizontal", "vertical"], axes=boxed);
 

 

 

#Convert and store xi and doCalc into an array
# ans1:= CodeTools:-Usage(Array([seq(seq( [j, u, doCalc(j,u)], j=-2..0, 1), u =0..2,1)])):

#

with(ListTools):

testproc := proc(j, u)
   global calcvals:=doCalc(j,u);
   evalf(calcvals[4]), evalf(Im(calcvals[7]));
end proc;

proc (j, u) global calcvals; calcvals := doCalc(j, u); evalf(calcvals[4]), evalf(Im(calcvals[7])) end proc

(2)

with(plots):
setcolors(["BlueViolet","Coral"]):
whatjs:=[op([seq(1..4)]),op([seq(6..10)])]:

CodeTools:- Usage(listcontplot([seq([seq(testproc(i,j)[1], i=-2..0,0.1)],j=1..3)],filledregions=true));

CodeTools:- Usage(listcontplot([seq([seq(testproc(i,j)[2], i=-2..0, 0.1)],j=1..3)],filledregions=true));

memory used=20.42GiB, alloc change=32.00MiB, cpu time=2.50m, real time=4.84m, gc time=13.31s

 

 

CodeTools:- Usage(listcontplot3d([seq([seq(testproc(i,j)[1], i=-2..0, 0.1)],j=1..3)],filledregions=true));
 

memory used=20.42GiB, alloc change=0 bytes, cpu time=3.92m, real time=6.55m, gc time=19.97s

 

 

CodeTools:- Usage(listcontplot3d([seq([seq(testproc(i,j)[2], i=-2..0, 0.1)],j=1..3)],filledregions=true));

memory used=20.41GiB, alloc change=0 bytes, cpu time=2.51m, real time=4.99m, gc time=13.38s

 

 

CodeTools:- Usage(matrixplot([seq([seq(testproc(i,j)[2], i=-2..0, 0.1)],j=1..3)],filledregions=true));

#for i from -2 to 0 by 0.1 do
#   for j from 1 to 10 do
#       testproc(i,j);
#       print(i, j);
# end do;
# end do;

NULL

NULL


 

Download Case_3_IjuptilK_300522.mw

@tomleslie 

 

Thank you very much.

@tomleslie

I really understand what is happening. In case, I have to use contourplot instead of listconplot.  Is there any way I can change.

CodeTools:- Usage(listconplot([seq([seq(testproc(i,j)[2], i=-2..0,0.1)],j=1..3)],filledregions=true));

 to something say

Usage(contourplot([seq([seq(testproc(i,j)[2], i=-2..0,0.1)],j=1..3)],filledregions=true));

Because, I want the x-axis to be -2..0 and y 1..3. 

I tried it but it returned an error.

@tomleslie 

Okay, thank you.

@Carl Love

Okay, thank you.

@tomleslie

Hi Tomleslie,

## please why do I have repeated values for j and u? And how can solve this problem without using Makeunique function.
## I also want to do contourplot for u vs j for complex values of p_min. 

In the attached file you will see what I'm talking about at the end of the file.

Case_3_IjuptilK_300522.mw

@tomleslie

Okay, thank you very much. 

@tomleslie

 

Okay, thanks 

@tomleslie 

 

Hi Tom,

Please, how can I place the legend in the center of the curve?

@tomleslie 

 

Okay, thanks 

@tomleslie

Hi Tom,

Please how can I change the scale of a plot: For instance, on the theta_plot, I wanna change 0.00005 to 50 micrometer on the x-axis and convert the 0.0001 to 1*10^-4(standard form) on the y-axis.

Scaling.mw

@tomleslie 

Thank you for your assistance.

@tomleslie

Hi Tom, 

You are right, I want to plot the imaginary parts of the two complex numbers (- 23.3886061400*I and 23.3886061400*I) against the xi in the list.

[-10, -9.6000000000, -9.2000000000, -8.8000000000, -8.0000000000, -7.6000000000, -7.2000000000, -6.8000000000, -6.0000000000, -5.6000000000, -5.2000000000, -4.8000000000, -4.4000000000, -4.0000000000, -3.6000000000, -2.8000000000, -2.4000000000, -0.4000000000]

The idea is that once we know the xi that gives us the minimum value of Re(p_n), then we can plot the Im(p_n) for that xi values.

Note that imaginary for the other xi values are zeros except for the two above.  The reason for doing this is that I want to find out if the Im(p_n) is ever non-zero for these values of xi listed above.

@tomleslie 

Sorry again, how do I plot the Im(p_n) for the values of xi we got from this: nstart := [ seq( `if`( member( complex(extended_numeric), whattype~({entries(ans[j,5], nolist)})), ans[j,1], NULL), j=1..op([2,1,2], ans))];

Because I added the seq( Im(gamma1*alpha/(4*k__1*qq[i]^2/d^2-alpha3*xi/eta__1)), i=1..m)  to the return list, and I got 5 entries in the ans. 

And if I want to plot nstart vs convert(ans[..,5], they don't have the same length.

xi = [-10, -9.6, -9.2, -8.8, -8.0, -7.6, -7.2, -6.8, -6.0, -5.6, -5.2, -4.8, -4.4, -4.0, -3.6, -2.8, -2.4, -.4] ( i.e the xi  for real p_n).

@tomleslie 

Thank you very much. I really appreciate this.

2 3 4 5 6 7 Page 4 of 7