acer

32333 Reputation

29 Badges

19 years, 322 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

Maple has no facilities for creating animation files (eg. mp4, mov, avi, etc) from multiple images.  Animated gif from plots is all it has.

Yes, this is a sorry state of affairs. It's been requested many times over the past two decades. (Sister product MapleSim can do mpeg-2 export of some simulations.)

The ImageTools:-Embed command can only insert once per Execution Group or Document Block. Later instances within the same Group will overwrite the contents of the embedded Task region. That GUI limitation is shared with Explore, Tabulate, InsertContent, and a few other things. IIRC it is documented in a few places. 

Did you want to export it an a single column Matrix? For example,

restart

Digits := 8:

with(plots):

with(CurveFitting):

with(plottools):

with(ExcelTools):

v := .7:

Disp := 20:

esp := 1000000:

k := 0:

E := proc (x, t) options operator, arrow; Int(exp((-esp*w^4+Disp*w^2+k)*t)*cos(w*(x+v*t))/Pi, w = 0 .. infinity, epsilon = 0.1e-6) end proc;

proc (x, t) options operator, arrow; Int(exp((-esp*w^4+Disp*w^2+k)*t)*cos(w*(x+v*t))/Pi, w = 0 .. infinity, epsilon = 0.1e-6) end proc

f := proc (x) options operator, arrow; 15.5*exp(-(1/3710000)*(x-12590)^2)+14.55*exp(-(1/3000000)*(x-16100)^2) end proc:

NULL

u := proc (x, t) options operator, arrow; Int(E(x-xi, t)*f(xi), xi = 0 .. 0.2e5, epsilon = 0.1e-4) end proc;

proc (x, t) options operator, arrow; Int(E(x-xi, t)*f(xi), xi = 0 .. 0.2e5, epsilon = 0.1e-4) end proc

NULL

NULL

uu := evalf(Int(E(0-xi, i)*f(xi), xi = 0 .. 20000, method = _NCrule, epsilon = 10^(-6)))

Int((Int(.31830988*exp((-1000000.*w^4+20.*w^2)*i)*cos(w*(-1.*xi+.7*i)), w = 0. .. Float(infinity)))*(15.5*exp(-0.26954178e-6*(xi-12590.)^2)+14.55*exp(-0.33333333e-6*(xi-16100.)^2)), xi = 0. .. 20000.)

M := [seq(evalf(Int(E(0-xi, 39000-i)*f(xi), xi = 0 .. 20000, method = _NCrule, epsilon = 10^(-6))), i = 10000 .. 20000, 250)]:

plots:-listplot(M);

ExportMatrix(cat(kernelopts(homedir), "/Documents/mapleprimes/M.txt"), `<|>`(`<,>`(M)));

435

 

``

Download getD1_ac.mw

The GUI does not send individual commands in the same session to separate kernel instances. A "session" implies a kernel, communicating input and output back and forth with the interface.

The GUI does not control or run try...catch, the kernel does.

If the kernel instance crashes then there's nothing else (on top of that, controlling the computation, with access to the session's computational history).

You need to wait for MapleSim 2020, if you want it to work properly with Maple 2020.

MapleSim 2019 works properly with Maple 2019, and not Maple 2020.

MapleSim 2020 won't be available for another few weeks.

You don't have to suppress output (using full colons, say), if running this in the GUI.

Instead you could either,
1) Set  interface(typesetting=standard)
or,
2) Retain interface(typesetting=extended)  which is default, while
    also setting  interface(prettyprint=1)

writeto_gui.mw

The second argument passed to plots:-animate should be a list, the entries of which get passed to the first argument (here, the plot command).

But you are trying to use the parametric calling sequence of plot, ie, a list of component formulas and a range.

So the second argument passed to plots:-animate should be a list whose first entry is that list (of component fomulas and range).

I've adjusted the range of the animating parameter. You may need to adjust the lines' component formulas or its own range so that the lines are the appropriate lengths (I'm guessing you want them to end at the outer curve intersection points).

restart;

with(plots):

rec:=(fx,fy)->arctan(-diff(fy,s),diff(fx,s)):

x[0]:=0.5*(s-2)^2:
y[0]:=-s^3-0.5*(s-0.5)^2:
a[0]:=-1.5:
b[0]:=1.5:
w[0]:=0.2:

x[1]:=s:
y[1]:=s^3+2*s^2:
a[1]:=-2.3:
b[1]:=1.2:
w[1]:=0.2:

eq := 1:

c := plot([x[eq],y[eq],s=a[eq]..b[eq]],color=blue):

l := plot([x[eq]+w[eq]*sin(rec(x[eq],y[eq])),y[eq]+w[eq]*cos(rec(x[eq],y[eq])),s=a[eq]..b[eq]]):

r := plot([x[eq]-w[eq]*sin(rec(x[eq],y[eq])),y[eq]-w[eq]*cos(rec(x[eq],y[eq])),s=a[eq]..b[eq]]):

m := diff(y[eq],s)/diff(x[eq],s)

3*s^2+4*s

B:=-1.5:
line := plot([x[eq], subs(s=B,m) * (x[eq]-subs(s=B,x[eq])) + subs(s=B,y[eq]), s=B..B+0.5], color=red):

[x[eq], subs(s=A,m) * (x[eq]-subs(s=A,x[eq])) + subs(s=A,y[eq]), s=A..A+0.5];

[s, (3*A^2+4*A)*(s-A)+A^3+2*A^2, s = A .. A+.5]

line := animate(plot,[[x[eq], subs(s=A,m) * (x[eq]-subs(s=A,x[eq])) + subs(s=A,y[eq]), s=A..A+0.5]],
                A=-2.3..1.2):

display(c,l,r,line);

 

Download agentpath_ac.mw

Is this what you mean, in your followup question? (I find your description somewhat unclear.)

xx := [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]:
yy := [3, 3, 3, 4, 5, 6, 8, 7, 6, 5]:
s := CurveFitting:-Spline(xx, yy, v):

plots:-animate(plot3d,[[s, th, v], v=0..y, th=0..2*Pi,
                       coords=cylindrical, style=surface],
               y=0..9);

Similarly, it seems to me as if your first example could be done in a straightforward manner.

restart;
s := exp(-0.1*x)*(2+sin(x)):

plots:-animate(plot3d,[[s, th, x], x=0..y, th=0..2*Pi,
                       coords=cylindrical, style=surface],
               y=0..6*Pi, frames=60);

Naturally, you can easily adjust the number of frames, surface style, title, or color, etc.

This produces a list of the 12 Matrices, ie. the 3d data for the vertices of each of the 5-sided polygons.

You can export those Matrices as you see fit.

export_ac.mw

If that's not the kind of data that you are after then you should explain properly what you want -- in a followup Reply/Comment here, not a duplicate Question. 

 

# Approach 1)  0.63 seconds real-time
#
# Utilize unapply to avoid invoking J0 and P1 for
# every numeric value of s.
# Also force NAG d01akc method for oscillatory integrand,
# which also has the effect of preventing evalf(Int(...))
# from trying to discern (symbolically) discontinuity points
# of the integrand.
#
# note: This passes Int calls to plot, which will evalf them.
#
restart;
P1:=(r,R)->(2/Pi)*(arccos(r/(2*R))-(r/(2*R))*sqrt(1-(r/(2*R))^2)):
J0:=(r,shk)-> BesselJ(0, 2*Pi*r*shk):
Jhk:=unapply((1/s)*Int(P1(r,R)*J0(r,shk)*sin(2*Pi*r*s), r=0..2*R,
                        method=_d01akc),[s,shk,R]);
CodeTools:-Usage(plot(s->Jhk(s,2.14,38), 0..5)):

proc (s, shk, R) options operator, arrow; (Int(2*(arccos((1/2)*r/R)-(1/4)*r*(4-r^2/R^2)^(1/2)/R)*BesselJ(0, 2*Pi*r*shk)*sin(2*Pi*r*s)/Pi, r = 0 .. 2*R, method = _d01akc))/s end proc

memory used=11.34MiB, alloc change=0 bytes, cpu time=623.00ms, real time=629.00ms, gc time=0ns

# Approach 2)   2.1 seconds real-time
#
# Utilize unapply to avoid invoking J0 and P1 for
# every numeric value of s.
# Utilize another unevaluated 'unapply' of the integrand
# within the Int call. This has the effect of preventing
# evalf(Int(...)) from trying to discern (symbolically)
# discontinuity points of the integrand. Notice that each
# call to Jhk will invoke unapply -- the overhead of that
# is greater than that of the Approach 1), and less
# than that of the Approach 3) below.
#
# note: This passes Int calls to plot, which will evalf them.
#
restart;
P1:=(r,R)->(2/Pi)*(arccos(r/(2*R))-(r/(2*R))*sqrt(1-(r/(2*R))^2)):
J0:=(r,shk)-> BesselJ(0, 2*Pi*r*shk):
Jhk:=unapply((1/s)*Int('unapply'(P1(r,R)*J0(r,shk)*sin(2*Pi*r*s),r), 0..2*R),[s,shk,R]);
CodeTools:-Usage(plot(s->Jhk(s,2.14,38), 0..5)):

proc (s, shk, R) options operator, arrow; (Int(unapply(2*(arccos((1/2)*r/R)-(1/4)*r*(4-r^2/R^2)^(1/2)/R)*BesselJ(0, 2*Pi*r*shk)*sin(2*Pi*r*s)/Pi, r), 0 .. 2*R))/s end proc

memory used=33.43MiB, alloc change=14.98MiB, cpu time=2.10s, real time=2.10s, gc time=0ns

# Approach 3)   4.8 seconds real-time
#
# Much like Kitonum's suggestion (given at end), in
# that the Jhk is the same.
#
# Utilize unapply to avoid invoking J0 and P1 for
# every numeric value of s.
#
# Here we retain the evalf within that unapply call. The
# effect is that the Pi within the sin call becomes an
# explicit float. This only mitigates the overhead of
# evalf(Int(...)) trying to discern (symbolically) discontinuity
# points of the integrand.
#
# note: The evalf call is not present within Jhk, and it is
# done here too by plot itself.
#
restart;
P1:=(r,R)->(2/Pi)*(arccos(r/(2*R))-(r/(2*R))*sqrt(1-(r/(2*R))^2)):
J0:=(r,shk)-> BesselJ(0, 2*Pi*r*shk):
Jhk:=unapply(evalf((1/s)*Int(P1(r,R)*J0(r,shk)*sin(2*Pi*r*s), r=0..2*R)),[s,shk,R]);
CodeTools:-Usage(plot(s->Jhk(s,2.14,38), 0..5)):

proc (s, shk, R) options operator, arrow; (Int(.6366197722*(arccos(.5000000000*r/R)-.2500000000*r*(4.-1.*r^2/R^2)^(1/2)/R)*BesselJ(0., 6.283185308*r*shk)*sin(6.283185308*r*s), r = 0. .. 2.*R))/s end proc

memory used=418.61MiB, alloc change=25.99MiB, cpu time=5.02s, real time=4.77s, gc time=485.94ms

# Approach 4)   3.9 seconds real-time
#
# Much like Kitonum's suggestion, and also much like the
# previous attempt, in that the Jhk is the same.
#
# note: I've added adaptive=false alongside numpoints=200, since
# the latter alone will not disable adaptive plotting and *might*
# consume more more time that with (default) numpoints=100.
# We could add   numpoints=200, adaptive=false  to the previous
# attempts, to make them a little faster, but they are fast
# enough and adaptive plotting can be useful for plots with
# these qualitative aspects.
#
# note: The evalf call is not present within Jhk, and it is
# done here too by plot itself.
#
restart;
P1:=(r,R)->(2/Pi)*(arccos(r/(2*R))-(r/(2*R))*sqrt(1-(r/(2*R))^2)):
J0:=(r,shk)->BesselJ(0, 2*Pi*r*shk):
Jhk:=unapply(evalf((1/s)*Int(P1(r,R)*J0(r,shk)*sin(2*Pi*r*s), r=0..2*R)),s,shk,R);
CodeTools:-Usage(plot(Jhk(s,2.14,38), s=0..5, numpoints=200, adaptive=false)):

proc (s, shk, R) options operator, arrow; (Int(.6366197722*(arccos(.5000000000*r/R)-.2500000000*r*(4.-1.*r^2/R^2)^(1/2)/R)*BesselJ(0., 6.283185308*r*shk)*sin(6.283185308*r*s), r = 0. .. 2.*R))/s end proc

memory used=329.58MiB, alloc change=25.99MiB, cpu time=4.09s, real time=3.90s, gc time=399.44ms

# Kitonum's suggestion)  4.8 seconds real-time
#
# See explanation for Approach 3)
#
restart;
P1:=(r,R)->(2/Pi)*(arccos(r/(2*R))-(r/(2*R))*sqrt(1-(r/(2*R))^2)):
J0:=(r,shk)->BesselJ(0, 2*Pi*r*shk):
Jhk:=unapply(evalf((1/s)*Int(P1(r,R)*J0(r,shk)*sin(2*Pi*r*s), r=0..2*R)),s,shk,R);
CodeTools:-Usage(plot(Jhk(s,2.14,38), s=0..5, numpoints=200, size=[1000,300])):

proc (s, shk, R) options operator, arrow; (Int(.6366197722*(arccos(.5000000000*r/R)-.2500000000*r*(4.-1.*r^2/R^2)^(1/2)/R)*BesselJ(0., 6.283185308*r*shk)*sin(6.283185308*r*s), r = 0. .. 2.*R))/s end proc

memory used=418.58MiB, alloc change=25.99MiB, cpu time=5.03s, real time=4.76s, gc time=526.53ms

 

Download unapply_quadrature_examples.mw

[edited] I forgot to provide the following simpler revision to the code originally posted.

The point of the unapply in this variant is not to avoid repeated calls to J0 and P1 (which will still occur for each separate numeric value of s from plot). The point of it here is that it calls evalf(Int(...)) with an operator for the integrand, rather than an expression. That makes the integrand more like a black-box, and prevents evalf(Int(...)) from incurring large overhead of examing the integrand for discontinuties in r.

Note that here I uneval-quote the call to Jhk in the plot call, to avoid premature evaluation. I could also have plotted  with plot(s->Jhk(s,2.14,38, 0..5) .

This is pretty fast, though not quite as fast as with the oscillatory NAG method forced as well.

restart;

P1:=(r,R)->(2/Pi)*(arccos(r/(2*R))-(r/(2*R))*sqrt(1-(r/(2*R))^2)):
J0:=(r,shk)-> BesselJ(0, 2*Pi*r*shk):
Jhk:=(s,shk,R)-> evalf((1/s)*Int(unapply(P1(r,R)*J0(r,shk)*sin(2*Pi*r*s), r), 0..2*R));
CodeTools:-Usage( plot('Jhk'(s,2.14,38), s=0..5) );

proc (s, shk, R) options operator, arrow; evalf((Int(unapply(P1(r, R)*J0(r, shk)*sin(2*Pi*r*s), r), 0 .. 2*R))/s) end proc

memory used=44.62MiB, alloc change=52.89MiB, cpu time=2.11s, real time=2.09s, gc time=98.09ms

 

Download _unapply_quadrature_followup.mw

It would be useful (and much simpler to explain) if evalf(Int(...)) simply got a new keyword option to forcibly disable discont checks. It's not always easy or even possible to get that effect if a NAG method is not suitable or if there is deep nesting of integrals or procedure calls.

You can use the view option to set the vertical range (or both horizontal and vertical ranges). You can use that with the plot and the implicitplot commands. For implicitplot there is also a convenient rangeasview option which forces the view to match the supplied ranges.

restart;
with(plots):

a,b:=7/10,0.2:

implicitplot(y=a*x+b, x=-10..10, y=-10..10,
             axis=[gridlines=[10,color=blue]],
             rangeasview);

implicitplot(y=a*x+b, x=-10..10, y=-10..10,
             axis=[gridlines=[10,color=blue]],
             view=[-10..10,-10..10]);

plot(a*x+b, x=-10..10, view=[-10..10,-10..10],
     axis=[gridlines=[10,color=blue]]);

plot(a*x+b, x=-10..10, view=[default,-10..10],
     axis=[gridlines=[10,color=blue]]);

plot(a*x+b, x=-10..10, view=-10..10,
     axis=[gridlines=[10,color=blue]]);

I suspect that the sky would not necessarily fall down if the Maple kernel developers allowed symbolic Matrix indexing to not throw an error.

restart

 

First, using a Matrix (capitalized), which we happen to
construct from the palette (but that aspect is not central).

 

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

Matrix(%id = 18446883814761169974)

b := m[a, 1]

Error, bad index into Matrix

b := 'm[a, 1]'

m[a, 1]

a := 1

1

b

1

a := 2

2

b

2

 

Note that attempting to evaluate m[a,1] will throw an error
if a is unassigned. The uneval-quoted reference to m[a,1] can

be tricky thing to manage, as you have to guard against

unintended evaluation.

 

a := 'a'

a

b

Error, bad index into Matrix

restart

 

Next, using a table-based matrix (uncapitalized).

Notice that the indexed reference m[a,1] does not produce an error
even when a is unassigned.

 

m := matrix(2, 2, [1, 3, 2, 4])

array( 1 .. 2, 1 .. 2, [( 1, 1 ) = (1), ( 1, 2 ) = (3), ( 2, 1 ) = (2), ( 2, 2 ) = (4)  ] )

b := m[a, 1]

m[a, 1]

a := 1

1

b

1

a := 2

2

b

2

a := 'a'

a

b

m[a, 1]

 

Download matrix_ref.mw

First lets see that PolynomialFit can indeed produce a better result, if its svdtolerance option is set fine enough.

restart;

kernelopts(version);

`Maple 2019.2, X86 64 LINUX, Oct 30 2019, Build ID 1430966`

Rok := Vector([2013, 2014, 2015, 2016, 2017, 2018], datatype = float):
TrzbyCelkemEmco := Vector([1028155, 1134120, 1004758, 929584, 995716, 1152042], datatype = float):

W := sort(Statistics:-PolynomialFit(3, Rok, TrzbyCelkemEmco, x, svdtolerance=1e-20));

HFloat(17490.364890056557)*x^3-HFloat(1.0573703094228968e8)*x^2+HFloat(2.1307569435843292e11)*x-HFloat(1.4312624246619506e14)

p1:=plot(W,x=min(Rok)..max(Rok),color=blue):
p2:=plot(Rok,TrzbyCelkemEmco,style=point,symbol=solidcircle,symbolsize=15):
plots:-display(p1,p2,view=min(TrzbyCelkemEmco)..max(TrzbyCelkemEmco),size=[500,400]);

 

Download svdtol.mw

Next, let's examine why the rank determination under PolynomialFit may have numeric difficulties using defaults. And we can do the least-squares solving using LinearAlgebra and a suitable raised working precision (Digits).

And then, at end, we'll see that CurveFitting:-LeastSquares can also find a similar result if Digits is set higher.

restart;

Rok := Vector([2013, 2014, 2015, 2016, 2017, 2018], datatype = float):
TrzbyCelkemEmco := Vector([1028155, 1134120, 1004758, 929584, 995716, 1152042], datatype = float):

pol3:=a*x^3+b*x^2+c*x+d;

a*x^3+b*x^2+c*x+d

resid:=eval~(pol3,x=~Rok)-TrzbyCelkemEmco;

Vector(6, {(1) = 8157016197.*a+4052169.000*b+2013.0*c+d-1028155.000, (2) = 8169178744.*a+4056196.000*b+2014.0*c+d-1134120.000, (3) = 8181353375.*a+4060225.000*b+2015.0*c+d-1004758.000, (4) = 8193540096.*a+4064256.000*b+2016.0*c+d-929584.0000, (5) = 8205738913.*a+4068289.000*b+2017.0*c+d-995716.0000, (6) = 8217949832.*a+4072324.000*b+2018.0*c+d-1152042.000})

A,B:=LinearAlgebra:-GenerateMatrix(convert(resid,list),[a,b,c,d]);

A, B := Matrix(6, 4, {(1, 1) = 8157016197., (1, 2) = 4052169.000, (1, 3) = 2013., (1, 4) = 1, (2, 1) = 8169178744., (2, 2) = 4056196.000, (2, 3) = 2014., (2, 4) = 1, (3, 1) = 8181353375., (3, 2) = 4060225.000, (3, 3) = 2015., (3, 4) = 1, (4, 1) = 8193540096., (4, 2) = 4064256.000, (4, 3) = 2016., (4, 4) = 1, (5, 1) = 8205738913., (5, 2) = 4068289.000, (5, 3) = 2017., (5, 4) = 1, (6, 1) = 8217949832., (6, 2) = 4072324.000, (6, 3) = 2018., (6, 4) = 1}), Vector(6, {(1) = 1028155.000, (2) = 1134120.000, (3) = 1004758.000, (4) = 929584., (5) = 995716., (6) = 1152042.000})

svd:=LinearAlgebra:-SingularValues(A);

Vector(6, {(1) = 0.2005517356e11, (2) = 8431.40030424569, (3) = 0.30315542273334386e-2, (4) = 0.9831935525e-9, (5) = 0., (6) = 0.})

Digits:=1+ceil(log[10](svd[1]/svd[4]));

21

V:=LinearAlgebra:-LeastSquares(A,B,method='SVD');
L:=Equate([a,b,c,d],convert(V,list)):
UU:=eval(pol3,L);

Vector(4, {(1) = 17490.3611109121787648, (2) = -105737008.1, (3) = 0.2130756483e12, (4) = -0.1431262115e15})

17490.3611109121787648*x^3-105737008.082130186122*x^2+213075648264.510954728*x-143126211485817.832367

p1UU:=plot(UU,x=min(Rok)..max(Rok)):
p2:=plot(Rok,TrzbyCelkemEmco,color=blue):
plots:-display(p1UU,p2,view=min(TrzbyCelkemEmco)..max(TrzbyCelkemEmco));

# Compare with CurveFitting:-LeastSquares,
# done using Digits=20 and Digits=21.
#
Digits:=20;
U:=CurveFitting:-LeastSquares(Rok,TrzbyCelkemEmco,x,curve=pol3);
p1U:=plot(U,x=min(Rok)..max(Rok)):
p2:=plot(Rok,TrzbyCelkemEmco,color=blue):
plots:-display(p1U,p2,view=min(TrzbyCelkemEmco)..max(TrzbyCelkemEmco));

20

55417.578996418964626+37231332.937900379590*x-36946.127657452668301*x^2+9.1658905419567079826*x^3

Digits:=21;
U:=CurveFitting:-LeastSquares(Rok,TrzbyCelkemEmco,x,curve=pol3);
p1U:=plot(U,x=min(Rok)..max(Rok)):
p2:=plot(Rok,TrzbyCelkemEmco,color=blue):
plots:-display(p1U,p2,view=min(TrzbyCelkemEmco)..max(TrzbyCelkemEmco));

21

-143126211486530.025986+213075648265.589385689*x-105737008.082662396656*x^2+17490.3611109998946169*x^3

sort(U);
sort(UU);

17490.3611109998946169*x^3-105737008.082662396656*x^2+213075648265.589385689*x-143126211486530.025986

17490.3611109121787648*x^3-105737008.082130186122*x^2+213075648264.510954728*x-143126211485817.832367

 

Download svd_fit.mw

The behavior which you describe with the "yyyyy" characters is its being rendered in italics, because the following space indicates to the parser that the preceding yyyyy is a name.

You can disable that behavior of rendering input variable names in italics in the CLI by issuing the command,
    interface(ansiedit=false):

But note that doing so will disable some colorizing of input, which is on by default.

See the ?interface Help page.

With a differently simplified form of the chareq expression (using Maple 18.02). The form of U2 below is quite a bit more compact, although I did not stress test any of the forms for the degree of susceptibility to numeric roundoff error in any range of En.

I had just finished writing this when I saw that dharr had submitted a second answer, also with a suggestion about NextZero (although, used with care, Calculus1:-Roots can also find more). I suppose a third Answer doesn't hurt.

restart

V0 := 1;

1

a := 2;

2

b := 1;

1

ks := sqrt(2*En);

2^(1/2)*En^(1/2)

qs := sqrt(2*(V0-En));

(2-2*En)^(1/2)

psi11 := proc (x) options operator, arrow; A1*exp(qs*x)+B1*exp(-qs*x) end proc

proc (x) options operator, arrow; A1*exp(qs*x)+B1*exp(-qs*x) end proc

psi21 := proc (x) options operator, arrow; A2*exp(I*ks*x)+B2*exp(-I*ks*x) end proc

proc (x) options operator, arrow; A2*exp(I*ks*x)+B2*exp(-I*ks*x) end proc

lp := a+b

3

psi12 := proc (x) options operator, arrow; exp(I*Kp*lp)*(A1*exp(qs*(x-lp))+B1*exp(-qs*(x-lp))) end proc

proc (x) options operator, arrow; exp(I*Kp*lp)*(A1*exp(qs*(x-lp))+B1*exp(-qs*(x-lp))) end proc

psi22 := proc (x) options operator, arrow; exp(I*Kp*lp)*(A2*exp(I*ks*(x-lp))+B2*exp(-I*ks*(x-lp))) end proc

proc (x) options operator, arrow; exp(I*Kp*lp)*(A2*exp(I*ks*(x-lp))+B2*exp(-I*ks*(x-lp))) end proc

eq1 := psi11(0) = psi21(0)

A1+B1 = A2+B2

eq2 := (D(psi11))(0) = (D(psi21))(0)

A1*(2-2*En)^(1/2)-B1*(2-2*En)^(1/2) = I*A2*2^(1/2)*En^(1/2)-I*B2*2^(1/2)*En^(1/2)

eq3 := psi21(a) = psi12(a)

A2*exp((2*I)*2^(1/2)*En^(1/2))+B2*exp(-(2*I)*2^(1/2)*En^(1/2)) = exp((3*I)*Kp)*(A1*exp(-(2-2*En)^(1/2))+B1*exp((2-2*En)^(1/2)))

eq4 := (D(psi21))(a) = (D(psi12))(a)

I*A2*2^(1/2)*En^(1/2)*exp((2*I)*2^(1/2)*En^(1/2))-I*B2*2^(1/2)*En^(1/2)*exp(-(2*I)*2^(1/2)*En^(1/2)) = exp((3*I)*Kp)*(A1*(2-2*En)^(1/2)*exp(-(2-2*En)^(1/2))-B1*(2-2*En)^(1/2)*exp((2-2*En)^(1/2)))

solve({eq1, eq2, eq3, eq4}, {A1, A2, B1, B2});

{A1 = 0, A2 = 0, B1 = 0, B2 = 0}

with(linalg):

row1 := vector([1, 1, -1, -1])

array( 1 .. 4, [( 1 ) = (1), ( 2 ) = (1), ( 3 ) = (-1), ( 4 ) = (-1)  ] )

row2 := vector([coeff((lhs-rhs)(eq2), A1), coeff((lhs-rhs)(eq2), B1), coeff((lhs-rhs)(eq2), A2), coeff((lhs-rhs)(eq2), B2)])

array( 1 .. 4, [( 1 ) = ((2-2*En)^(1/2)), ( 2 ) = (-(2-2*En)^(1/2)), ( 3 ) = (-I*2^(1/2)*En^(1/2)), ( 4 ) = (I*2^(1/2)*En^(1/2))  ] )

row3 := vector([coeff((lhs-rhs)(eq3), A1), coeff((lhs-rhs)(eq3), B1), coeff((lhs-rhs)(eq3), A2), coeff((lhs-rhs)(eq3), B2)])

array( 1 .. 4, [( 1 ) = (-exp((3*I)*Kp)*exp(-(2-2*En)^(1/2))), ( 2 ) = (-exp((2-2*En)^(1/2))*exp((3*I)*Kp)), ( 3 ) = (exp((2*I)*2^(1/2)*En^(1/2))), ( 4 ) = (exp(-(2*I)*2^(1/2)*En^(1/2)))  ] )

row4 := vector([coeff((lhs-rhs)(eq4), A1), coeff((lhs-rhs)(eq4), B1), coeff((lhs-rhs)(eq4), A2), coeff((lhs-rhs)(eq4), B2)])

array( 1 .. 4, [( 1 ) = (-exp((3*I)*Kp)*exp(-(2-2*En)^(1/2))*(2-2*En)^(1/2)), ( 2 ) = (exp((2-2*En)^(1/2))*exp((3*I)*Kp)*(2-2*En)^(1/2)), ( 3 ) = (I*2^(1/2)*En^(1/2)*exp((2*I)*2^(1/2)*En^(1/2))), ( 4 ) = (-I*2^(1/2)*En^(1/2)*exp(-(2*I)*2^(1/2)*En^(1/2)))  ] )

C := matrix(4, 4, 0)

array( 1 .. 4, 1 .. 4, [( 1, 2 ) = (0), ( 4, 2 ) = (0), ( 2, 2 ) = (0), ( 2, 1 ) = (0), ( 3, 4 ) = (0), ( 3, 1 ) = (0), ( 2, 3 ) = (0), ( 4, 4 ) = (0), ( 1, 4 ) = (0), ( 3, 2 ) = (0), ( 2, 4 ) = (0), ( 3, 3 ) = (0), ( 1, 1 ) = (0), ( 4, 1 ) = (0), ( 1, 3 ) = (0), ( 4, 3 ) = (0)  ] )

for i to 4 do C[1, i] := row1[i]; C[2, i] := row2[i]; C[3, i] := row3[i]; C[4, i] := row4[i] end do:

evalm(C)

array( 1 .. 4, 1 .. 4, [( 1, 2 ) = (1), ( 4, 2 ) = (exp((2-2*En)^(1/2))*exp((3*I)*Kp)*(2-2*En)^(1/2)), ( 2, 2 ) = (-(2-2*En)^(1/2)), ( 2, 1 ) = ((2-2*En)^(1/2)), ( 3, 4 ) = (exp(-(2*I)*2^(1/2)*En^(1/2))), ( 3, 1 ) = (-exp((3*I)*Kp)*exp(-(2-2*En)^(1/2))), ( 2, 3 ) = (-I*2^(1/2)*En^(1/2)), ( 4, 4 ) = (-I*2^(1/2)*En^(1/2)*exp(-(2*I)*2^(1/2)*En^(1/2))), ( 1, 4 ) = (-1), ( 3, 2 ) = (-exp((2-2*En)^(1/2))*exp((3*I)*Kp)), ( 2, 4 ) = (I*2^(1/2)*En^(1/2)), ( 3, 3 ) = (exp((2*I)*2^(1/2)*En^(1/2))), ( 1, 1 ) = (1), ( 4, 1 ) = (-exp((3*I)*Kp)*exp(-(2-2*En)^(1/2))*(2-2*En)^(1/2)), ( 1, 3 ) = (-1), ( 4, 3 ) = (I*2^(1/2)*En^(1/2)*exp((2*I)*2^(1/2)*En^(1/2)))  ] )

chareq := det(C);

(4*I)*2^(1/2)*exp((2-2*En)^(1/2))*(exp((3*I)*Kp))^2*exp(-(2-2*En)^(1/2))*(2-2*En)^(1/2)*En^(1/2)-(2*I)*2^(1/2)*exp(-(2*I)*2^(1/2)*En^(1/2))*exp((2-2*En)^(1/2))*exp((3*I)*Kp)*(2-2*En)^(1/2)*En^(1/2)-(2*I)*2^(1/2)*exp(-(2*I)*2^(1/2)*En^(1/2))*exp((3*I)*Kp)*exp(-(2-2*En)^(1/2))*(2-2*En)^(1/2)*En^(1/2)-(2*I)*2^(1/2)*exp((2*I)*2^(1/2)*En^(1/2))*exp((2-2*En)^(1/2))*exp((3*I)*Kp)*(2-2*En)^(1/2)*En^(1/2)-(2*I)*2^(1/2)*exp((2*I)*2^(1/2)*En^(1/2))*exp((3*I)*Kp)*exp(-(2-2*En)^(1/2))*(2-2*En)^(1/2)*En^(1/2)+(4*I)*(2-2*En)^(1/2)*exp((2*I)*2^(1/2)*En^(1/2))*2^(1/2)*En^(1/2)*exp(-(2*I)*2^(1/2)*En^(1/2))-4*En*exp(-(2*I)*2^(1/2)*En^(1/2))*exp((2-2*En)^(1/2))*exp((3*I)*Kp)+4*En*exp(-(2*I)*2^(1/2)*En^(1/2))*exp((3*I)*Kp)*exp(-(2-2*En)^(1/2))+4*En*exp((2*I)*2^(1/2)*En^(1/2))*exp((2-2*En)^(1/2))*exp((3*I)*Kp)-4*En*exp((2*I)*2^(1/2)*En^(1/2))*exp((3*I)*Kp)*exp(-(2-2*En)^(1/2))+2*exp(-(2*I)*2^(1/2)*En^(1/2))*exp((2-2*En)^(1/2))*exp((3*I)*Kp)-2*exp(-(2*I)*2^(1/2)*En^(1/2))*exp((3*I)*Kp)*exp(-(2-2*En)^(1/2))-2*exp((2*I)*2^(1/2)*En^(1/2))*exp((2-2*En)^(1/2))*exp((3*I)*Kp)+2*exp((2*I)*2^(1/2)*En^(1/2))*exp((3*I)*Kp)*exp(-(2-2*En)^(1/2))

U1 := `assuming`([simplify(combine(evalc(eval(chareq, Kp = 0))))], [En > 0, En < 1]);

(4*I)*(-2*cos(2*2^(1/2)*En^(1/2))*exp(2*(2-2*En)^(1/2))*(1-En)^(1/2)*En^(1/2)+2*En*sin(2*2^(1/2)*En^(1/2))*exp(2*(2-2*En)^(1/2))-2*cos(2*2^(1/2)*En^(1/2))*(1-En)^(1/2)*En^(1/2)+4*(1-En)^(1/2)*En^(1/2)*exp((2-2*En)^(1/2))-sin(2*2^(1/2)*En^(1/2))*exp(2*(2-2*En)^(1/2))-2*En*sin(2*2^(1/2)*En^(1/2))+sin(2*2^(1/2)*En^(1/2)))*exp(-(2-2*En)^(1/2))

U1 := `assuming`([collect(simplify(simplify(combine(evalc(eval(chareq, Kp = 0)))), size), [exp, sin, cos], proc (u) options operator, arrow; simplify(simplify(combine(u)), size) end proc)], [En > 0, En < 1]);

(16*I)*(1-En)^(1/2)*En^(1/2)*exp(-(2-2*En)^(1/2))*exp((2-2*En)^(1/2))+(((-4*I+(8*I)*En)*sin(2*2^(1/2)*En^(1/2))-(8*I)*(1-En)^(1/2)*cos(2*2^(1/2)*En^(1/2))*En^(1/2))*exp(2*(2-2*En)^(1/2))+(4*I-(8*I)*En)*sin(2*2^(1/2)*En^(1/2))-(8*I)*(1-En)^(1/2)*cos(2*2^(1/2)*En^(1/2))*En^(1/2))*exp(-(2-2*En)^(1/2))

Student:-Calculus1:-Roots(U1, En = 0 .. 1, numeric);

[0., .2669813621]

plot(([Re, Im])(U1), En = 0. .. 1.0, view = -20 .. 20, style = [line, point], thickness = 5, color = [red, blue]);

fImU1 := unapply(Im(U1), En):

.2669813620

U2 := `assuming`([simplify(combine(evalc(eval(chareq, Kp = 0))))], [En > 1]);

8*(-1+En)^(1/2)*En^(1/2)*cos(2*2^(1/2)*En^(1/2)-(-2+2*En)^(1/2))+8*(-1+En)^(1/2)*En^(1/2)*cos(2*2^(1/2)*En^(1/2)+(-2+2*En)^(1/2))-8*En*cos(2*2^(1/2)*En^(1/2)-(-2+2*En)^(1/2))+8*En*cos(2*2^(1/2)*En^(1/2)+(-2+2*En)^(1/2))+4*cos(2*2^(1/2)*En^(1/2)-(-2+2*En)^(1/2))-4*cos(2*2^(1/2)*En^(1/2)+(-2+2*En)^(1/2))-16*(-1+En)^(1/2)*En^(1/2)

U2 := `assuming`([collect(simplify(simplify(combine(evalc(eval(chareq, Kp = 0)))), size), [exp, sin, cos], proc (u) options operator, arrow; simplify(simplify(combine(u)), size) end proc)], [En > 1]);

(8*(-1+En)^(1/2)*En^(1/2)-8*En+4)*cos(2*2^(1/2)*En^(1/2)-(-2+2*En)^(1/2))+(8*(-1+En)^(1/2)*En^(1/2)+8*En-4)*cos(2*2^(1/2)*En^(1/2)+(-2+2*En)^(1/2))-16*(-1+En)^(1/2)*En^(1/2)

Student:-Calculus1:-Roots(U2, En = 1 .. 20, numeric);

[2.375133245, 2.712788471, 9.049238832, 9.174900984]

fU2 := unapply(U2, En):

2.375133244

2.712788471

9.049238831

9.174900983

Digits := 200:

plot(U2, En = 1 .. 1000, size = [600, 300]);

 

Download qm_maple_-_periodic_potentials_acc.mw

Two ways, using a symbolic option (instead of assumptions):

restart;

kernelopts(version);

`Maple 2019.2, X86 64 LINUX, Nov 26 2019, Build ID 1435526`

combine(exp(k*(ln(t) + ln(a))) - exp(ln(t) + ln(a))^k, symbolic);
simplify(%);

exp(k*ln(t*a))-(t*a)^k

0

expand(exp(k*(ln(t) + ln(a))) - exp(ln(t) + ln(a))^k);
simplify(%,symbolic);

t^k*a^k-(t*a)^k

0

 

NULL

Download symbolic_ac.mw

First 132 133 134 135 136 137 138 Last Page 134 of 336