acer

32303 Reputation

29 Badges

19 years, 309 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

Sure, you do not have to write out manually all the applications of the individual trig functions to the input values.

There are several ways to do it. You could use the Matrix and Vector commands, with initializers (see Help). You could use map or seq across lists/Vectors of the functions and inputs.

Or you could apply the functions to the inputs using elementwise operations. Below I show an example of how that may construct the Matrix data of the DataFrame. It's not the only way.

The embedded GUI Table doesn't render properly on this site. That's just a Mapleprimes problem. In the actual Maple GUI it looks fine (centered, with appropriate size, etc), like your original. What I changed is how to programmatically construct the DataFrame.

restart

NumericEventHandler('division_by_zero' = proc (f, ops, defVal) NumericStatus('division_by_zero' = false); undefined end proc)

 

NULL

F := [sin, cos, tan]; C := [0, (1/6)*Pi, (1/4)*Pi, (1/3)*Pi, (1/2)*Pi]

df := DataFrame(Matrix(`~`[`@`(Vector, F)](C)), columns = C, rows = F(x))

Tabulate(df, width = 50.)

 

 

0

(1/6)*Pi

(1/4)*Pi

(1/3)*Pi

(1/2)*Pi

sin(x)

0

1/2

(1/2)*sqrt(2)

(1/2)*sqrt(3)

1

cos(x)

1

(1/2)*sqrt(3)

(1/2)*sqrt(2)

1/2

0

tan(x)

0

(1/3)*sqrt(3)

1

sqrt(3)

undefined

 
 

TableauTest_ac.mw

You could also substitute something else for undefined, into the data. Eg,
TableauTest_ac_ne.mw

ps. Your question appears to be about how to programmatically and more easily construct the DataFrame (and its data/Matrix, and its headers). Your question doesn't really seem to be about how to use Tabulate to embed that DataFrame. I've adjusted the Questions's tags accordingly.

You can use the keystrokes Shift-Enter to move the input cursor down to another line within the loop (without causing execution).

In the Maple's documentation (Help) that is called a soft new line here, as one of its keyboard shortcuts.

You'll still need the statement terminators such as a semicolon. (Strictly speaking it's not needed on the last statement before the end do, but I find life easier to always have it.)

I recommend indenting your code involving do-loops, if-then, etc. It helps with legibility. How many spaces to indent is up to you, but being consistent helps.

[edit: I did the following before the OP suggested that the x-axis range he used is incorrect.]

Given the presence of the undefined values in the Matrix, you could do the following,

(minv,maxv):=[min[defined],max[defined]](data[gridji,4])[];

However, it seems that the Float(undefined) values might arise from your relatively poor approach to linear-solving for r, given A.r=B. You are computing the inverse of A and then multiplying. I get the following (Maple2016 the oldest I have at hand right now) if I instead use r:=LinearSolve(A,B).

[edit: It might even be appropropriate to perform a least-squares solve (LeastSquares with QR or SVD method, or LinearSolve with QR method) if the system is underdetermined. I don't know enough about what you're computing to pass judgement on that.]

ps. Every time you ask a Question here you omit the important detail that you are using Maple 18. (I end up adding it to your Question headers, if I remember to look at your .mw attachments -- in the cases that you provide such.) It'd be more helpful to the readership if you could bother to mark your Questions with your particular version number.

restart:

 

Procedure

 

NULL

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*alpha/(4*k__1*q^2/d^2-alpha__3*xi/eta__1 - chi*H^2),
                       n:= 50,
                       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,soln3, Imagroot1, Imagroot2, AreTherePurelyImaginaryRoots;

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

  g:= q-(1-alpha)*tan(q)+  c*tan(q) + w*q:
  f:= subs(q = x+I*y, g):
  b1:= evalc(Re(f)) = 0:
  b2:= y-(1-alpha)*tanh(y) -(alpha__3*xi*alpha/(eta__1*(4*k__1*y^2/d^2+alpha__3*xi/eta__1)))*tanh(y) = 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 complex roots

  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;

AreTherePurelyImaginaryRoots:= type(true, 'truefalse');
  try
# Compute the rest of the Roots
  soln3:= simplify~(fsolve({ b2}, { y = 0 .. infinity}),zero);  
  Imagroot1:=subs(soln3, I*y);
  Imagroot2:= -Imagroot1;
  catch:
   AreTherePurelyImaginaryRoots:= type(FAIL, 'truefalse');
  end try;

## odd and all asymptotes
  OddAsymptotes:= Vector[row]([seq(evalf((1/2)*(2*j+1)*Pi), j = 0 .. n)]);
  AllAsymptotes:= sort(Vector[row]([OddAsymptotes, qstar]));
       
 if (xi > 0) then
  if AreThereComplexRoots and not AreTherePurelyImaginaryRoots
  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 AreThereComplexRoots and AreTherePurelyImaginaryRoots
  then gg:= [qcomplex1, qcomplex2, Imagroot2, 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 and AreTherePurelyImaginaryRoots
  then gg:= [Imagroot2, 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 and not AreTherePurelyImaginaryRoots
  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:
else
if AreThereComplexRoots and not AreTherePurelyImaginaryRoots
  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 AreThereComplexRoots and AreTherePurelyImaginaryRoots
  then gg:= [qcomplex1, qcomplex2, Imagroot2, 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 and AreTherePurelyImaginaryRoots
  then gg:= [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 and not AreTherePurelyImaginaryRoots
  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:
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):
  r := LinearSolve(A,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):

## compute minimum value of tau
  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,qq;
     end proc:

Run the calculation for supplied value of 'xi

 

``

NULL

``

``

``

``

``

``

``

ans:=[doCalc(0.1, 2)]:
evalf(ans[9]);
evalf(ans[10]);

69.69853430

34.84926715

 

Contour plot for pmin (minimum  of tau) for different xi and H values

 

``

``

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

testproc := proc(j, u, k) option remember;
   local calcvals;
   if not [j,u,k]::[numeric,numeric,posint] then
      return 'procname'(args);
   end if;
   calcvals:=doCalc(j,u);
   evalf(calcvals[k]);
end proc:

NN := 8;

8

(gridji,rngj,rngi):=[5,5],0..3,-2.0..-0.0;
PP[gridji,4]:=plot3d(max(-1e2,testproc(i,j,4)), j=rngj, i=rngi, grid=gridji):
PP[gridji,7]:= plot3d(Im(testproc(i,j,7)), j=rngj, i=rngi, grid=gridji):

[5, 5], 0 .. 3, -2.0 .. -0.

## Surface plot using 3d
PP[gridji,4]:

## This is the pmin values
data[gridji,4]:=op([1,3],PP[gridji,4]);

data[[5, 5], 4] := Array(1..5, 1..5, {(1, 1) = -4.5610937267458205, (1, 2) = -16.551730624406634, (1, 3) = -8.369456800130445, (1, 4) = -9.140026175750366, (1, 5) = 2.951847102925982, (2, 1) = -4.407854759516567, (2, 2) = -14.697771941098084, (2, 3) = -7.867574545412465, (2, 4) = -100.0, (2, 5) = 8.746447899601446, (3, 1) = -45.14232589754586, (3, 2) = -11.00108026387169, (3, 3) = -100.0, (3, 4) = -35.20694911216856, (3, 5) = -59.05829652528201, (4, 1) = -16.596794315137668, (4, 2) = -7.751663635971151, (4, 3) = -21.534163510244372, (4, 4) = -15.03671683942269, (4, 5) = -18.171783445443022, (5, 1) = -8.8033427869475, (5, 2) = -39.10547901669723, (5, 3) = -10.022167745982093, (5, 4) = -8.344150581537837, (5, 5) = -9.227858698918745}, datatype = float[8])

## Convert Array to Matrix to allow the use listcontplot
ConMatrix := data[gridji,4];

ConMatrix := Array(1..5, 1..5, {(1, 1) = -4.5610937267458205, (1, 2) = -16.551730624406634, (1, 3) = -8.369456800130445, (1, 4) = -9.140026175750366, (1, 5) = 2.951847102925982, (2, 1) = -4.407854759516567, (2, 2) = -14.697771941098084, (2, 3) = -7.867574545412465, (2, 4) = -100.0, (2, 5) = 8.746447899601446, (3, 1) = -45.14232589754586, (3, 2) = -11.00108026387169, (3, 3) = -100.0, (3, 4) = -35.20694911216856, (3, 5) = -59.05829652528201, (4, 1) = -16.596794315137668, (4, 2) = -7.751663635971151, (4, 3) = -21.534163510244372, (4, 4) = -15.03671683942269, (4, 5) = -18.171783445443022, (5, 1) = -8.8033427869475, (5, 2) = -39.10547901669723, (5, 3) = -10.022167745982093, (5, 4) = -8.344150581537837, (5, 5) = -9.227858698918745}, datatype = float[8])

## Assign the min and max value of pmin to minv and maxv
(minv,maxv):=[min[defined],max[defined]](data[gridji,4])[];

-100., 8.74644789960144564

## GRB color and the contours
(color1,color2):="Red","Yellow":
conts:=[minv, seq(minv+(i-1)*(maxv-minv)/(NN+2-1),i=2..NN+1), maxv]/2.0;

 

[-50.00000000, -43.95853067, -37.91706134, -31.87559202, -25.83412269, -19.79265336, -13.75118404, -7.709714705, -1.668245380, 4.373223950]

colorlist:=ColorTools:-Color~(((c1,c2,N)->[seq(c1+(i-1)*(c2-c1)/(N-1),
                                               i=1..N)])([ColorTools:-Color(color2)[]],
                                                         [ColorTools:-Color(color1)[]],
                                                         nops(conts))):

##

## Use the Conmatrix and the conts to plot contours and transform the axes
subsindets(plots:-listcontplot(ConMatrix, contours=conts,filledregions,thickness=0),
           specfunc(anything,THICKNESS),u->THICKNESS(0.2)):
LCP[gridji,4]:=plottools:-transform((x,y)->[lhs(rngj)+(rhs(rngj)-lhs(rngj))*(x-1)/(gridji[1]-1),
                                            lhs(rngi)+(rhs(rngi)-lhs(rngi))*(y-1)/(gridji[2]-1)])(%):
display(LCP[gridji,4],
        seq(plot(-2.0,1..1,legend=evalf[4](conts[i]),thickness=15,
                 color=colorlist[nops(conts)-i+1]),
            i=1..nops(conts)),
        size=[500,400],legendstyle=[location=right], axes=boxed, labels = ["u", "xi"], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14]);

## This deals with the Imaginary part

## Surface plot using plot3d for Imginary part
PP[gridji,7]:

##

data[gridji,7]:=op([1,3],PP[gridji,7]):

## Convert Array to matrix
AMatrix := Matrix(data[gridji,7]);

AMatrix := Matrix(5, 5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = -.0, (3, 4) = .0, (3, 5) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (5, 1) = .0, (5, 2) = -.0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0}, datatype = float[8])

(minv,maxv):=[min,max](data[gridji,7])[];

0., 0.

(color1,color2):="Red","Yellow":
conts1:=[minv, seq(minv+(i-1)*(maxv-minv)/(NN+2-1),i=2..NN+1), maxv]/~1.5;

[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]

colorlist:=ColorTools:-Color~(((c1,c2,N)->[seq(c1+(i-1)*(c2-c1)/(N-1),
                                               i=1..N)])([ColorTools:-Color(color2)[]],
                                                         [ColorTools:-Color(color1)[]],
                                                         nops(conts1))):

subsindets(plots:-listcontplot(AMatrix,
                               contours=conts1-~1e-9,filledregions,thickness=0),
           specfunc(anything,THICKNESS),u->THICKNESS(0.2)):
LCP[gridji,7]:=plottools:-transform((x,y)->[lhs(rngj)+(rhs(rngj)-lhs(rngj))*(x-1)/(gridji[1]-1),
                                            lhs(rngi)+(rhs(rngi)-lhs(rngi))*(y-1)/(gridji[2]-1)])(%):
display(LCP[gridji,7],
        seq(plot(-2.0,1..1,legend=evalf[4](conts1[i]),thickness=15,
                 color=colorlist[nops(conts)-i+1]),
            i=1..nops(conts1)),
        size=[500,400],legendstyle=[location=right], axes=boxed, labels = ["u", "xi"], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14]);

 

NULL

Download Case3_Ijuptilk_Contourplot_pmin_290822_ac.mw

Is this the kind of thing you want to accomplish?

I did not change the way your formulas work, or which kinds of data structure you use. (Eg. it might be easier to create a as an Array, etc.)

[edit: I figured that you are trying to learn Maple coding, and that a direct movement of your original code into a procedure (which is what you asked) would help in your understanding. Some people learn better by focusing on only a few things at once. So I just put it into a reusable procedure, as is.

Look in particular at how the parameters and the locals of the procedure are all just symbols, and not any indexed names as you attempt had tried.]

restart;

Digits := 15;

15

T5Scheme:=proc(x1,y1,a,L)          
local i,k,nk,N,aa,x,y;
   x[1]:=x1; y[1]:=y1;
   N:=numelems(x[1]):         
   for k from 1 to L do              
     nk := 3^(k - 1)*(N - 6) + 6;         
     for i from 3 to nk - 2 do            
       x[k + 1][3*i - 8] := a[-6]*x[k][i + 2] + a[-3]*x[k][i + 1] + a[0]*x[k][i]
                            + a[3]*x[k][i - 1] + a[6]*x[k][i - 2];
       y[k + 1][3*i - 8] := a[-6]*y[k][i + 2] + a[-3]*y[k][i + 1] + a[0]*y[k][i]
                            + a[3]*y[k][i - 1] + a[6]*y[k][i - 2];
       x[k + 1][3*i - 7] := a[-5]*x[k][i + 2] + a[-2]*x[k][i + 1] + a[1]*x[k][i]
                            + a[4]*x[k][i - 1] + a[7]*x[k][i - 2];
       y[k + 1][3*i - 7] := a[-5]*y[k][i + 2] + a[-2]*y[k][i + 1] + a[1]*y[k][i]
                            + a[4]*y[k][i - 1] + a[7]*y[k][i - 2];
       x[k + 1][3*i - 6] := a[-4]*x[k][i + 2] + a[-1]*x[k][i + 1] + a[2]*x[k][i]
                            + a[5]*x[k][i - 1];
       y[k + 1][3*i - 6] := a[-4]*y[k][i + 2] + a[-1]*y[k][i + 1] + a[2]*y[k][i]
                            + a[5]*y[k][i - 1];
     end do;          
   end do;  
   return evalf(convert(x[L+1],list)), evalf(convert(y[L+1],list));
end proc:

## Initial polygon
x[1] := [-1, 1/2, 1, 0, -sqrt(2)/2, -1, 1/2, 1, 0, -sqrt(2)/2, -1, 1/2, 1, 0, -sqrt(2)/2];

[-1, 1/2, 1, 0, -(1/2)*2^(1/2), -1, 1/2, 1, 0, -(1/2)*2^(1/2), -1, 1/2, 1, 0, -(1/2)*2^(1/2)]

y[1] := [0, sqrt(3)/2, 0, -1, -sqrt(2)/2, 0, sqrt(3)/2, 0, -1, -sqrt(2)/2, 0, sqrt(3)/2, 0, -1, -sqrt(2)/2];

[0, (1/2)*3^(1/2), 0, -1, -(1/2)*2^(1/2), 0, (1/2)*3^(1/2), 0, -1, -(1/2)*2^(1/2), 0, (1/2)*3^(1/2), 0, -1, -(1/2)*2^(1/2)]

s2 := [-1, 1/2, 1, 0, -sqrt(2)/2, -1];

[-1, 1/2, 1, 0, -(1/2)*2^(1/2), -1]

t2 := [0, sqrt(3)/2, 0, -1, -sqrt(2)/2, 0];

[0, (1/2)*3^(1/2), 0, -1, -(1/2)*2^(1/2), 0]

assign(a[-6] = 13/1296, a[-5] = -11/648, a[-4] = -1/16, a[-3] = -107/1296,
       a[-2] = 179/1296, a[-1] = 9/16, a[0] = 137/144, a[1] = 137/144,
       a[2] = 9/16, a[3] = 179/1296, a[4] = -107/1296, a[5] = -1/16,
       a[6] = -11/648, a[7] = 13/1296);

s1,t1 := T5Scheme(x[1],y[1],a,3):

plot([<<s1> | <t1>>, <<s2> | <t2>>], linestyle = [1, 3], color = [black, red]);

Download zz123_ac.mw

It works for me, in the following sense of referencing a file on my local machine.

In the Hyperlink properties I set the Type to URL and I set the Target to,

    file:///home/acer/mapleprimes/collectexample.html

and clicking on that link in my Maple GUI document makes my web-browser open that local file. (Files with extension .html extension are associated with my web-browser application.)

The key thing above was to use the  file:  protocol in the URL.

It's not entirely clear to me what you want to happen by linking to a geogebra file. Something with your web-browser, or do you expect some 3rd party app to launch it, or...?

Look at the output of your assignment of an operator to the name f.

Clearly something has gone wrong in your 2D Input, since the A*sin(omega*t + varphi)*I term is missing.

I deleted and retyped that term,

omega := 5

5

`&varphi;` := (1/3)*Pi

(1/3)*Pi

A := 4

4

"f(t):=A*cos(omega*t+`&varphi;`)+A*I*sin(omega*t+`&varphi;`)"
 

proc (t) options operator, arrow, function_assign; A*cos(omega*t+varphi)+I*A*sin(omega*t+varphi) end proc

f(0)

2+(2*I)*3^(1/2)

A*exp(I*`&varphi;`)

2+(2*I)*3^(1/2)

NULL

Download Mapleprimes_Book_2_Question_3.6_ac.mw

For the purpose of plotting you may not need high accuracy in the numeric integration.

Also, if you forcibly take only the (significant) real part of calls to N, while also "disabling" overflow exceptions under evalhf and the _d01amc quadrature method, then the float evaluation of calls to J appears to speed up.

(Also, check that might simplified revision of L is correct...)

I didn't experiment with using a finite range of integration. (Maybe that could be faster, while still reasonably accurate? I didn't look...)

Perhaps you could tell us how long your original plot took to compute?

restart;

with(plots):

F:=kappa->kappa;

proc (kappa) options operator, arrow; kappa end proc

f:=(alpha,delta)->exp(-abs(F(kappa))^2*(1+delta^2)/2-abs(F(kappa))*alpha)/abs(F(kappa));

proc (alpha, delta) options operator, arrow; exp(-(1/2)*abs(F(kappa))^2*(1+delta^2)-abs(F(kappa))*alpha)/abs(F(kappa)) end proc

L:=(alpha,delta,Lambda) ->
  (lambda^2*exp(-alpha^2/2)/4)*(2*Int(f(alpha,delta),kappa= -infinity..-Lambda));

proc (alpha, delta, Lambda) options operator, arrow; (1/2)*lambda^2*exp(-(1/2)*alpha^2)*(Int(f(alpha, delta), kappa = -infinity .. -Lambda)) end proc

#0.0008209373770*lambda^2
forget(evalf);
CodeTools:-Usage( evalf(L(4,1,0.001)) );

memory used=1.57MiB, alloc change=0 bytes, cpu time=29.00ms, real time=30.00ms, gc time=0ns

0.8209373770e-3*lambda^2

g:=(beta,delta)->exp(-I*kappa*beta-abs(F(kappa))^2*(1+delta^2)/2)/abs(F(kappa));

proc (beta, delta) options operator, arrow; exp(-I*kappa*beta-(1/2)*abs(F(kappa))^2*(1+delta^2))/abs(F(kappa)) end proc

E:=(omega,gg)->exp(I*omega*gg)*(1-erf((gg+I*omega)/sqrt(2)));

proc (omega, gg) options operator, arrow; exp(I*omega*gg)*(1-erf((gg+I*omega)/sqrt(2))) end proc

MyHandler := proc(operator,operands,default_value)
               NumericStatus( overflow = false );
               return 10^10;
             end proc:
ig1template := proc(kappa)
              NumericEventHandler(overflow = MyHandler);
              evalhf(__dummy); end proc:

igdum:='Re'(simplify(exp(-alpha^2/2)/8*(g(beta,delta)*(E(abs(F(kappa)),gg)+E(abs(F(kappa)),-gg))))):
J := subs(__igdum=igdum, proc(alpha,delta,Lambda,beta,gg) local ig;
  ig := subs(__dummy= __igdum, eval(ig1template));
  evalf(lambda^2*abs(Int(ig,-infinity..-Lambda,epsilon=1e-4,method=_d01amc)
                     +Int(ig,Lambda..infinity,epsilon=1e-4,method=_d01amc)));
end proc):

forget(evalf);
CodeTools:-Usage( J(4,1,0.001,8,3) );

memory used=49.92MiB, alloc change=0 bytes, cpu time=319.00ms, real time=319.00ms, gc time=0ns

0.7304272433e-3*lambda^2

N := (beta,alpha) -> (J(alpha,1,0.001,beta,3)-L(alpha,1,0.001))/lambda^2;

proc (beta, alpha) options operator, arrow; (J(alpha, 1, 0.1e-2, beta, 3)-L(alpha, 1, 0.1e-2))/lambda^2 end proc

forget(evalf);
CodeTools:-Usage(contourplot(N, 0..10, 0..10, grid=[25,25]));

memory used=20.12GiB, alloc change=76.00MiB, cpu time=2.50m, real time=2.32m, gc time=19.72s

 

Download Negativity_v1_acc.mw

It might even be possible to improve the timing even further by parallelizing computation over a 2D Array (over beta and alpha values), and then using listcontplot.

@nm You might try it with a different web browser.

A month ago I was unable to submit postings with Chromium, but Firefox worked. That problem went away after a while, as mysteriously as it arrived.

You could use the parametric calling sequence of the plot command to interchange the roles of the axes.

Below, I don't have to create the first plot (in order to transform ie, say). I include it only for illustration of the effect of the second plot, which is easily and directly constructed with a single command.

restart

NULL

y := 3*x^2-4*x+5

3*x^2-4*x+5

plot(y, x, axes = box, labels = ["x", "y"])

plot([y, x, x = -10 .. 10], axes = box, labels = ["y", "x"]);

 

Download plot_graph_ac.mw

Have you tried uploading your .mw file to the Maple Cloud, and accessing it online (visually) from there? (See also here.)

You should even be able to login to that cloud from your Maple version in which you create such .mw files, and upload directly from that. You can upload to a public area, or make a private group for you and your students.

I do not know whether there is an emulator or similar within the ipad that can alloow it to run the (free) MaplePlayer. That's why I suggested the Maple Cloud, since you can access that from your webbrowser. However the quality of reproduction in the Maple Cloud is not always as good as that in the Maple Player.

The filled-regions version of LineChart is AreaChart.

For example,

restart;

with(Statistics):

A := <1, 2, 3, 4, 2, 3, 1>:
B := <1, 5, 1, 0, 2, 1, 1>:
C := <1, 2, 1, 3, 1, 4, 2>:

AreaChart([A, B, C], format = stacked,
          color=[white,pink,orange], transparency=0.3);

plots:-display(
  AreaChart([A, B, C], format = stacked,
          color=[white,pink,orange], transparency=0.3),
  LineChart([A, B, C], format = stacked));

 

Download AreaChart.mw

Transposing a column Vector seems an awkward way to construct a row Vector.

A:=`<|>`(2,7);

Vector[row](2, {(1) = 2, (2) = 7})

Vector[row]([1,A]);

Vector[row](3, {(1) = 1, (2) = 2, (3) = 7})

Vector[row]([A,1]);

Vector[row](3, {(1) = 2, (2) = 7, (3) = 1})

Download rowV_augment.mw

You could also do it in these alternate ways:

`<|>`(1,seq(A));
`<|>`(seq(A),1);

or,

<1|seq(A)>;
<seq(A)|1>;

I didn't bother to figure out which is most efficient if done millions of times.

If you do it as,

  B := A[1..4, 2..2];

instead of your,

  B := A(1..4, 2);

then you get B as a 4x1 Matrix instead of your 4x1 Column Vector.

restart

A := Matrix(4, 4, {(1, 1) = m[1, 1], (1, 2) = m[1, 2], (1, 3) = m[1, 3], (1, 4) = m[1, 4], (2, 1) = m[2, 1], (2, 2) = m[2, 2], (2, 3) = m[2, 3], (2, 4) = m[2, 4], (3, 1) = m[3, 1], (3, 2) = m[3, 2], (3, 3) = m[3, 3], (3, 4) = m[3, 4], (4, 1) = m[4, 1], (4, 2) = m[4, 2], (4, 3) = m[4, 3], (4, 4) = m[4, 4]})

A := Matrix(4, 4, {(1, 1) = 0, (1, 2) = m[1, 2], (1, 3) = m[1, 3], (1, 4) = m[1, 4], (2, 1) = 0, (2, 2) = m[2, 2], (2, 3) = m[2, 3], (2, 4) = m[2, 4], (3, 1) = 0, (3, 2) = m[3, 2], (3, 3) = m[3, 3], (3, 4) = m[3, 4], (4, 1) = 0, (4, 2) = m[4, 2], (4, 3) = m[4, 3], (4, 4) = m[4, 4]})

A(1 .. 4, 1) := Matrix(4, 1):

A

Matrix([[0, m[1, 2], m[1, 3], m[1, 4]], [0, m[2, 2], m[2, 3], m[2, 4]], [0, m[3, 2], m[3, 3], m[3, 4]], [0, m[4, 2], m[4, 3], m[4, 4]]])

B := A[1 .. 4, 2 .. 2]

B := Matrix(4, 1, {(1, 1) = m[1, 2], (2, 1) = m[2, 2], (3, 1) = m[3, 2], (4, 1) = m[4, 2]})

A.B

4

F := Matrix(4, 1, {(1, 1) = 0.4279668887e-7, (2, 1) = -0.3901148183e-7, (3, 1) = 0.3900685346e-7, (4, 1) = 0.})

A.B-F

Matrix(4, 1, {(1, 1) = m[1, 2]*m[2, 2]+m[1, 3]*m[3, 2]+m[1, 4]*m[4, 2]-0.4279668887e-7, (2, 1) = m[2, 2]^2+m[2, 3]*m[3, 2]+m[2, 4]*m[4, 2]+0.3901148183e-7, (3, 1) = m[3, 2]*m[2, 2]+m[3, 3]*m[3, 2]+m[3, 4]*m[4, 2]-0.3900685346e-7, (4, 1) = m[2, 2]*m[4, 2]+m[3, 2]*m[4, 3]+m[4, 2]*m[4, 4]})

NULL

Download soalzarb_ac.mw

Hopefully I have coded "King's" method as you intend. (You can check that.)

The first worksheet's way is simpler, but more expensive. It uses plots:-densityplot. It colors (hue) the points by the iteration count. But it does not illustrate to which root each input point converges. It's not very fast, even at a relatively small maximal-iteration limit.  King_iter.mw

The second worksheet colors by the root to which the input point converges.
King_vs_Newt_basin.mw

restart;


This is not a fast approach, for modest resolution.

It consists of creating a density plot of the iteration count
to achieve a given forward-error.

The color shading is by iteration count.

__FF := proc(z) option inline; z^3-1 end proc:

 

KingCount := subs(g=3, F=__FF, dF=D(__FF),
  proc(x0,y0) local u,i,v,A;
  u[1] := evalf(x0+y0*I);
  for i from 1 to 25 do
    v[i] := u[i] - F(u[i])/dF(u[i]);
    A[i] := F(u[i])+(g-2)*F(v[i]);
    u[i+1] := v[i] - 1/A[i]*(F(u[i])+g*F(v[i]))*F(v[i])/dF(u[i]);
    if abs(F(u[i+1]))<1e-3 then return i; end if;
  end do:
  return i-1;
end proc):

NewtCount := subs(g=3, F=__FF, dF=D(__FF),
  proc(x0,y0) local u,i,v,A;
  u[1] := evalf(x0+y0*I);
  for i from 1 to 25 do
    u[i+1] := u[i] - F(u[i])/dF(u[i]);
    if abs(F(u[i+1]))<1e-3 then return i; end if;
  end do:
  return i-1;
end proc):

 

PK := CodeTools:-Usage(
        plots:-densityplot(KingCount, -2..2, -2..2, grid=[201,201],
                           title="King's method", size=[300,300],
                           colorstyle=HUE, style=surface) ):

memory used=2.05GiB, alloc change=36.00MiB, cpu time=12.19s, real time=12.19s, gc time=1.48s

PN := CodeTools:-Usage(
        plots:-densityplot(NewtCount, -2..2, -2..2, grid=[201,201],
                           title="Newton's method", size=[300,300],
                           colorstyle=HUE, style=surface) ):

memory used=0.74GiB, alloc change=0 bytes, cpu time=4.74s, real time=4.75s, gc time=615.31ms

plots:-display(Array([PK, PN]));

 

 

 

 

 


An alternative way to get a similar image of the iteration count
(using Newton's method) is as follows. This is much faster.
The hue spread is slightly differently utilized.

M := CodeTools:-Usage(
       Fractals:-EscapeTime:-Newton(300, -2 - 2*I, 2 + 2*I, z^3-1,
                                    iterationlimit=50, tolerance=1e-3) ):
plot(-2..2,-2..2,background=M);

memory used=19.73MiB, alloc change=19.62MiB, cpu time=342.00ms, real time=248.00ms, gc time=25.96ms

 

King_iter.mw

And here is the second:

restart:


The following uses the Escape command from the IterativeMaps package.

It is reasonable fast (even with some relatively ineffient post-processing.

I have chosen to post-process the output data so that the coloring is
made according to the root to which the given x+y*I input point converges.

Optionally we can scale the saturation by the iteration count.

ee := z^3 - 1;

z^3-1

## For ee a polynomial we can compute all the roots.
## This is used to specify the colors, below.
rts:=[fsolve(ee,z,complex)]:
evalf[7](%);

[-.5000000-.8660254*I, -.5000000+.8660254*I, 1.]

#(minr,maxr) := (min,max)(Re~(rts)):
#(mini,maxi) := (min,max)(Im~(rts)):
(minr,maxr) := -1, 1;
(mini,maxi) := -1, 1;

-1, 1

-1, 1

# One iteration of Newton's method
dee := diff(ee, z):
Newt1 := eval( z - ee/dee, z=zr+zi*I);

zr+I*zi-(1/3)*((zr+I*zi)^3-1)/(zr+I*zi)^2

# One iteration of King's method
F := unapply(ee,z): dF:=D(F):
v := u - F(u)/dF(u):
A := F(u)+(g-2)*F(v):
uu := v - 1/A*(F(u)+g*F(v))*F(v)/dF(u):
King1 := eval(uu, [g=3, u=zr+zi*I]);

zr+I*zi-(1/3)*((zr+I*zi)^3-1)/(zr+I*zi)^2-(1/3)*((zr+I*zi)^3-4+3*(zr+I*zi-(1/3)*((zr+I*zi)^3-1)/(zr+I*zi)^2)^3)*((zr+I*zi-(1/3)*((zr+I*zi)^3-1)/(zr+I*zi)^2)^3-1)/(((zr+I*zi)^3-2+(zr+I*zi-(1/3)*((zr+I*zi)^3-1)/(zr+I*zi)^2)^3)*(zr+I*zi)^2)

H := proc(Z,N,maxiter,{scalesaturation::truefalse:=false})
  local final,G,h,w,img,temp,temp2,newt,fzi,fzr;
  uses IterativeMaps, ImageTools;
  fzi := evalc( Im( Z ) );
  fzi := convert(simplify(fzi),horner);
  fzr := evalc( Re( subs(zi=zit, Z) ) );
  fzr := convert(simplify(fzr),horner);
  h,w := N, floor(N*(maxr-minr)/(maxi-mini));
  newt := Escape( height=h, width=w,
                       [zi, zr, zit],
                       [fzi, fzr, zi],
                       [y, x, y], evalc(abs(eval(ee,[z=zr+zit*I])))^2<1e-10,
                       2*minr, 2*maxr, 2*mini, 2*maxi,
                       iterations = maxiter,
                       redexpression = zr, greenexpression = zit, blueexpression=i );
  img:=Array(1..h,1..w,1..3,(i,j,k)->1.0,datatype=float[8]);
  temp:=ln~(map[evalhf](max,copy(newt[..,..,3]),1.0));
  temp2:=Threshold(temp,0.0001,low=1);
  img[..,..,3]:=FitIntensity(zip(`*`,temp,temp2));
  ## Done inefficiently, set the hue by position in the list of
  ## all roots
  G:=Array(zip((a,b)->min[index](map[evalhf](abs,rts-~(a+b*I))),
               newt[..,..,1],newt[..,..,2]),datatype=float[8],order=C_order):
  img[..,..,1]:=240*FitIntensity(G):
  if scalesaturation then
    img[..,..,2]:=1-~img[..,..,3]:
  end if;
  final := HSVtoRGB(img):
end proc:

 

Newtimg := H(Newt1,400,200):

# scale down the saturation by the time to converge
NewtimgS := H(Newt1,400,200,scalesaturation):

ImageTools:-Embed([Newtimg,NewtimgS]);

    

Kingimg := H(King1,400,200):

# scale down the saturation by the time to converge
KingimgS := H(King1,400,200,scalesaturation):

ImageTools:-Embed([Kingimg,KingimgS]);

   

King_vs_Newt_basin.mw

The names Column and Row are also exports of the LinearAlgebra package.

So if you load LinearAlgebra after loading DocumentTools:-Layout then you rebind those two names to LinearAlgebra -- but your code was written earlier to utilize them from DocumentTools:-Layout, and so you've broken your code.

There are many ways to fix up this muddle. Too many for me to show them all.

Here's one way, with minor editing of your attempt: loop_Alger_ac.mw

Personally, I'd rather use procedures (and uses), instead of all those loose top-level with calls to begin your Document. But that might be too much of a change (to your preferred methodology) for you to accept.

First 62 63 64 65 66 67 68 Last Page 64 of 336