Applications, Examples and Libraries

Share your work here

I have just posted an article with this title at Maplesoft Application Center here.
It was motivated by a question posed by  Markiyan Hirnyk  here and a test problem proposed there by Kitonum.

Now I just want to give the promissed complete solution to Kitonum's test:

Compute the plane area of the region defined by the inequalities:

R := [ (x-4)^2+y^2 <= 25, x^2+(y-3)^2 >= 9, (x+sqrt(7))^2+y^2 >= 16 ];

plots:-inequal(R, x=-7..10, y=-6..6, scaling=constrained);

The used procedures (for details see the mentioned article):

ranges:=proc(simpledom::list(relation), X::list(name))
local rez:=NULL, r,z,k,r1,r2;
if nops(simpledom)<>2*nops(X) then error "Domain not simple!" fi;
for k to nops(X) do    r1,r2:=simpledom[2*k-1..2*k][]; z:=X[k];
  if   rhs(r1)=z and lhs(r2)=z then rez:=z=lhs(r1)..rhs(r2),rez; #a<z,z<b
  elif lhs(r1)=z and rhs(r2)=z then rez:=z=lhs(r2)..rhs(r1),rez  #z<b,a<z
  else error "Strange order in a simple domain" fi
od;
rez
end proc:

MultiIntPoly:=proc(f, rels::list(relation(ratpoly)), X::list(name))
local r,rez,sol,irr,wirr, rels1, w;
irr:=[indets(rels,{function,realcons^realcons})[]];
wirr:=[seq(w[i],i=1..nops(irr))];
rels1:=eval(rels, irr=~wirr);
sol:=SolveTools:-SemiAlgebraic(rels1,X,parameters=wirr):
sol:=remove(hastype, eval(sol,wirr=~irr), `=`); 
add(Int(f,ranges(r,X)),r=sol)
end proc:

MeasApp:=proc(rel::{set,list}(relation), Q::list(name='range'(realcons)), N::posint)
local r, n:=0, X, t, frel:=evalf(rel)[];
if indets(rel,name) <> indets(Q,name)  then error "Non matching variables" fi;
r:=[seq(rand(evalf(rhs(t))), t=Q)];
X:=[seq(lhs(t),t=Q)];
to N do
  if evalb(eval(`and`(frel), X=~r())) then n:=n+1 fi;
od;
evalf( n/N*mul((rhs-lhs)(rhs(t)),t=Q) );
end proc:

Problem's solution:

MultiIntPoly(1, R, [x,y]):  # Unfortunately it's slow; patience needed!
radnormal(simplify(value(%)));

evalf(%) = MeasApp(R, [x=-7..10,y=-6..6], 10000); # A rough numerical check
           61.16217534 = 59.91480000

 

This presentation is about magnetic traps for neutral particles, first achieved for cold neutrons and nowadays widely used in cold-atom physics. The level is that of undergraduate electrodynamics and tensor calculus courses. Tackling this topic within a computer algebra worksheet as shown below illustrates well the kind of advanced computations that can be done today with the Physics package. A new feature minimizetensorcomponents and related functionality is used along the presentation, that requires the updated Physics library distributed at the Maplesoft R&D Physics webpage.
 

 

Magnetic traps in cold-atom physics

 

Pascal Szriftgiser1 and Edgardo S. Cheb-Terrab2 

(1) Laboratoire PhLAM, UMR CNRS 8523, Université Lille 1, F-59655, France

(2) Maplesoft

 

We consider a device constructed with a set of electrical wires fed with constant electrical currents. Those wires can have an arbitrary complex shape. The device is operated in a regime such that, in some region of interest, the moving particles experience a magnetic field that varies slowly compared to the Larmor spin precession frequency. In this region, the effective potential is proportional to the modulus of the field: LinearAlgebra[Norm](`#mover(mi("B"),mo("&rarr;"))`(x, y, z)), this potential has a minimum and, close to this minimum, the device behaves as a magnetic trap.

 

 

 

Figure 1: Schematic representation of a Ioffe-Pritchard magnetic trap. It is made of four infinite rods and two coils.

_________________________________________

 

Following [1], we show that:

 

  

a) For a time-independent magnetic field  `#mover(mi("B"),mo("&rarr;"))`(x, y, z) in vacuum, up to order two in the relative coordinates X__i = [x, y, z] around some point of interest, the coefficients of orders 1 and 2 in this expansion, `v__i,j` and `c__i,j,k` , respectively the gradient and curvature, contain only 5 and 7 independent components.

  

b) All stationary points of LinearAlgebra[Norm](`#mover(mi("B"),mo("&rarr;"))`(x, y, z))^2 (nonzero minima and saddle points) are confined to a curved surface defined by det(`&PartialD;`[j](B[i])) = 0.

  

c) The effective potential, proportional to LinearAlgebra[Norm](`#mover(mi("B"),mo("&rarr;"))`(x, y, z)), has no maximum, only a minimum.

 

Finally, we draw the stationary condition surface for the case of the widely used Ioffe-Pritchard magnetic trap.

  

 

  

Reference

  

[1] R. Gerritsma and R. J. C. Spreeuw, Topological constraints on magnetostatic traps,  Phys. Rev. A 74, 043405 (2006)

  

 

The independent components of `v__i,j` and `c__i,j,k` entering B[i] = u[i]+v[i, j]*X[j]+(1/2)*c[i, j, k]*X[j]*X[k]

   

The stationary points are within the surface det(`&PartialD;`[j](B[i])) = 0

   

U = LinearAlgebra[Norm](`#mover(mi("B",fontweight = "bold"),mo("&rarr;",fontweight = "bold"))`)^2 has only minima, no maxima

   

Drawing the Ioffe-Pritchard Magnetic Trap

   


 

MagneticTraps.mw or in pdf format with the sections open: MagneticTraps.pdf

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

Much of this topic is developed using traditional techniques. Maple modernizes and optimizes solutions by displaying the necessary operators and simple commands to solve large problems. Using the conditions of equilibrium for both moment and force we find the forces and moments of reactions for any type of structure. In spanish.

Equlibrium.mw

https://www.youtube.com/watch?v=7zC8pGC4F2c

Lenin Araujo Castillo

Ambassador of Maple

Good book to start studying maple for engineering.

 


 

restart; with(plots)

Australopithecus := [[75, 25], [97, 30], [93, 40], [93, 45], [83, 50], [80, 55], [79, 60], [81, 73], [74, 76], [68, 81], [60, 82], [50, 83], [40, 80], [30, 71], [25, 60], [24, 50], [25, 37], [15, 33], [10, 30], [45, 10], [55, 16], [65, 10], [80, 8], [93, 14], [96, 24]]:

man := [[95, 39], [113, 40], [111, 47], [118, 53], [113, 62], [109, 72], [112, 88], [112, 95], [107, 112], [99, 117], [85, 122], [72, 122], [49, 117], [36, 104], [31, 78], [39, 52], [43, 43], [44, 34], [39, 16], [73, 3], [81, 17], [98, 14], [105, 17], [104, 26], [111, 33]]:

morph := proc (poly1, poly2, t) if nops(poly1) <> nops(poly2) then ERROR("mensaje.") end if; [seq([(1-t)*op(1, op(k, poly1))+t*op(1, op(k, poly2)), (1-t)*op(2, op(k, poly1))+t*op(2, op(k, poly2))], k = 1 .. nops(poly1))] end proc:

display([seq(polygonplot(morph(Australopithecus, man, (1/20)*k), scaling = constrained), k = 0 .. 19)], insequence = true, axes = none);

 

NULL


 

Download Australopithecus_updated.mw

http://www.gatewaycoalition.org/includes/display_project.aspx?ID=279&maincatid=105&subcatid=1019&thirdcatid=0

Lenin Araujo Castillo

Ambassador of Maple

In this application you can visualize the impulse generated by a constant and variable force for the interaction of a particle with an object in a state of rest or movement. It is also the calculation of the momentum-momentum equation by entering the mass of the particle to solve initial and final velocities respectively according to the case study. Engineering students can quickly display the calculations and then their interpretation. In spanish.

Plot_of_equation_impulse-momentum.mw

Exercises_of_Momentum-Impulse_Linear.mw

Lenin Araujo Castillo

Ambassador of Maple

So I have recently finished up a project that took different sounds found in nature, and through the Spectrogram command, plotted the frequency of each sound over time with some really cool results!

https://www.maplesoft.com/applications/view.aspx?SID=154346 

The contrast between sounds produced by the weather such as tornadoes, thunder, and hail versus something as innocuous as a buzzing bee, a chorus of crickets, or a purring cat really shows the variance in the different sounds we hear in our day to day life, while also creating some very interesting imagery.

My personal favourite was the cricket chorus, producing a very ordered image with some really cool spikes through many different frequencies as the crickets chirped, as shown here:

Using this plot, we can do some interesting things, like count the number of chirps in 8 seconds, which turns out to be 18-18.5. Why Is this important? Well, there’s a law known as Dolbear’s Law(shown here: https://en.wikipedia.org/wiki/Dolbear%27s_law)  which uses the number of chirps in a minute for Fahrenheit, or 8 seconds for Celsius to calculate the temperature. Celsius is very simple, and just requires adding 5 to the number of chirps in 8 seconds, which gets us a temperature of 23C.

Tc= 5 + N8

For Fahrenheit, it’s a bit more complicated, as we need the chirps in a minute. This is around 132 chirps in our case. Putting that into the formula:

TF= 50 +((N60 – 40)/4)

Which gets us 73F, or 22.7C, so you can see that it works pretty well! Pretty cool, huh?

 

There was also some really cool images that were produced, like the thunder plot:

Which I personally really like due to the contrasting black and yellow spike that occurs. Overall this was a very fun project to do, getting to tweak the different colours and scales of each spectrogram, creating a story out of a sound. Hope you all enjoy it!

 I accidentally stumbled on this problem in the list of tasks for mathematical olympiads. I quote its text in Russian-English translation:

"The floor in the drawing room of Baron Munchausen is paved with the identical square stone plates.
 Baron claims that his new carpet (made of one piece of a material ) covers exactly 24 plates and
 at the same time each vertical and each horizontal row of plates in the living room contains 
exactly 4 plates covered with carpet. Is not the Baron deceiving?"

At first glance this seems impossible, but in fact the Baron is right. Several examples can be obtained simply by hand, for example

                                        or        

 

The problem is to find all solutions. This post is dedicated to this problem.

We put in correspondence to each such carpet a matrix of zeros and ones, such that in each row and in each column there are exactly 2 zeros and 4 ones. The problem to generate all such the matrices was already discussed here and Carl found a very effective solution. I propose another solution (based on the method of branches and boundaries), it is less effective, but more universal. I've used this method several times, for example here and here.
There will be a lot of such matrices (total 67950), so we will impose natural limitations. We require that the carpet be a simply connected set that has as its boundary a simple polygon (non-self-intersecting).

Below we give a complete solution to the problem.


restart;
R:=combinat:-permute([0,0,1,1,1,1]);
# All lists of two zeros and four units

# In the procedure OneStep, the matrices are presented as lists of lists. The procedure adds one row to each matrix so that in each column there are no more than 2 zeros and not more than 4 ones

OneStep:=proc(L::listlist)
local m, k, l, r, a, L1;
m:=nops(L[1]); k:=0;
for l in L do
for r in R do
a:=[op(l),r];
if `and`(seq(add(a[..,j])<=4, j=1..6)) and `and`(seq(m-add(a[..,j])<=2, j=1..6)) then k:=k+1; L1[k]:=a fi;
od; od;
convert(L1, list);
end proc:

# M is a list of all matrices, each of which has exactly 2 zeros and 4 units in each row and column

L:=map(t->[t], R):
M:=(OneStep@@5)(L):
nops(M);

                                            67950

M1:=map(Matrix, M):

# From the list of M1 we delete those matrices that contain <1,0;0,1> and <0,1;1,0> submatrices. This means that the boundaries of the corresponding carpets will be simple non-self-intersecting curves

k:=0:
for m in M1 do
s:=1;
for i from 2 to 6 do
for j from 2 to 6 do
if (m[i,j]=0 and m[i-1,j-1]=0 and m[i,j-1]=1 and m[i-1,j]=1) or (m[i,j]=1 and m[i-1,j-1]=1 and m[i,j-1]=0 and m[i-1,j]=0) then s:=0; break fi;
od: if s=0 then break fi; od:
if s=1 then k:=k+1; M2[k]:=m fi;
od:
M2:=convert(M2, list):
nops(M2);

                                             394

# We find the list T of all segments from which the boundary consists

T:='T':
n:=0:
for m in M2 do
k:=0: S:='S':
for i from 1 to 6 do
for j from 1 to 6 do
if m[i,j]=1 then
if j=1 or (j>1 and m[i,j-1]=0) then k:=k+1; S[k]:={[j-1/2,7-i-1/2],[j-1/2,7-i+1/2]} fi;
if i=1 or (i>1 and m[i-1,j]=0) then k:=k+1; S[k]:={[j-1/2,7-i+1/2],[j+1/2,7-i+1/2]} fi;
if j=6 or (j<6 and m[i,j+1]=0) then k:=k+1; S[k]:={[j+1/2,7-i+1/2],[j+1/2,7-i-1/2]} fi;
if i=6 or (i<6 and m[i+1,j]=0) then k:=k+1; S[k]:={[j+1/2,7-i-1/2],[j-1/2,7-i-1/2]} fi; 
fi;
od: od:
n:=n+1; T[n]:=[m,convert(S,set)];
od:
T:=convert(T, list):

# Choose carpets with a connected border

C:='C': k:=0:
for t in T do
a:=t[2]; v:=op~(a);
G:=GraphTheory:-Graph([$1..nops(v)], subs([seq(v[i]=i,i=1..nops(v))],a));
if GraphTheory:-IsConnected(G) then k:=k+1; C[k]:=t fi;
od:
C:=convert(C,list):
nops(C);
                                             
 208

# Sort the list of border segments so that they go one by one and form a polygon

k:=0: P:='P':
for c in C do
a:=c[2]: v:=op~(a);
G1:=GraphTheory:-Graph([$1..nops(v)], subs([seq(v[i]=i,i=1..nops(v))],a));
GraphTheory:-IsEulerian(G1,'U');
U; s:=[op(U)];
k:=k+1; P[k]:=[seq(v[i],i=s[1..-2])];
od:
P:=convert(P, list):

# We apply AreIsometric procedure from here to remove solutions that coincide under a rotation or reflection

P1:=[ListTools:-Categorize( AreIsometric, P)]:
nops(P1);

                                                 28


We get 28 unique solutions to this problem.

Visualization of all these solutions:

interface(rtablesize=100):
E1:=seq(plottools:-line([1/2,i],[13/2,i], color=red),i=1/2..13/2,1):
E2:=seq(plottools:-line([i,1/2],[i,13/2], color=red),i=1/2..13/2,1):
F:=plottools:-polygon([[1/2,1/2],[1/2,13/2],[13/2,13/2],[13/2,1/2]], color=yellow):
plots:-display(Matrix(4,7,[seq(plots:-display(plottools:-polygon(p,color=red),F, E1,E2), p=[seq(i[1],i=P1)])]), scaling=constrained, axes=none, size=[800,700]);

 

 

Carpet1.mw

The code was edited.

 

 

Using the syntax in Maple we develop the energy with conservation equations here we are applying the commands int, factor, solve among others. We also integrate vector functions through the scalar product and finally we calculate conservative fields applying the rotational to a field of force. Exclusive for engineering students. In spanish.

Work_of_a_Force.mw

Lenin Araujo castillo

Ambassafor of Maple


 

Maple 2017.1 and 2017.2 introduced several improvements in the solution of PDE & Boundary conditions problems (exact solutions).  Maple 2017.3 includes more improvements in this same area.

 

The following is a set of 25 examples of different PDE & Boundary Conditions problems that are solvable in Maple 2017.3 but not in Maple 2017.2 or previous releases. In the examples that follow, in some cases the PDE is different, in other cases the boundary conditions are of a different kind or solving the problem involves different computational strategies.

As usual, at the end there is a link pointing to the corresponding worksheet.

 

pde[1] := diff(u(x, t), t) = k*(diff(u(x, t), x, x)); bc[1] := u(x, 0) = 6+4*cos(3*Pi*x/L), (D[1](u))(0, t) = 0, (D[1](u))(L, t) = 0

diff(u(x, t), t) = k*(diff(diff(u(x, t), x), x))

 

u(x, 0) = 6+4*cos(3*Pi*x/L), (D[1](u))(0, t) = 0, (D[1](u))(L, t) = 0

(1)

`assuming`([pdsolve([pde[1], bc[1]], u(x, t))], [L > 0, k > 0])

u(x, t) = 6+4*cos(3*Pi*x/L)*exp(-9*k*Pi^2*t/L^2)

(2)

pde[2] := diff(g(t, x), t) = diff(g(t, x), x, x)+a*g(t, x); bc[2] := g(0, x) = 1

diff(g(t, x), t) = diff(diff(g(t, x), x), x)+a*g(t, x)

 

g(0, x) = 1

(3)

pdsolve([pde[2], bc[2]])

g(t, x) = exp(a*t)

(4)

pde[3] := diff(u(x, t), t) = diff(u(x, t), x, x)

bc[3] := u(x, 0) = f(x), u(-1, t) = 0, u(1, t) = 0

diff(u(x, t), t) = diff(diff(u(x, t), x), x)

 

u(x, 0) = f(x), u(-1, t) = 0, u(1, t) = 0

(5)

pdsolve([pde[3], bc[3]], u(x, t))

u(x, t) = Sum((Int(f(x)*sin(n*Pi*x), x = -1 .. 1))*sin(n*Pi*x)*exp(-Pi^2*n^2*t), n = 1 .. infinity)

(6)

pde[4] := diff(u(x, t), t) = k*(diff(u(x, t), x, x)); bc[4] := u(0, t) = 0, u(L, t) = 0, u(x, 0) = piecewise(0 < x and x <= (1/2)*L, 1, (1/2)*L < x and x < L, 2)

pde[4] := diff(u(x, t), t) = k*(diff(u(x, t), x, x))

 

u(0, t) = 0, u(L, t) = 0, u(x, 0) = piecewise(0 < x and x <= (1/2)*L, 1, (1/2)*L < x and x < L, 2)

(7)

`assuming`([pdsolve([pde[4], bc[4]], u(x, t))], [L > 0])

u(x, t) = Sum((2*cos((1/2)*Pi*n)+2+4*(-1)^(1+n))*sin(n*Pi*x/L)*exp(-k*Pi^2*n^2*t/L^2)/(Pi*n), n = 1 .. infinity)

(8)

pde[5] := diff(u(x, t), t) = k*(diff(u(x, t), x, x))

bc[5] := (D[1](u))(0, t) = 0, (D[1](u))(L, t) = 0, u(x, 0) = -3*cos(8*Pi*x/L)

diff(u(x, t), t) = k*(diff(diff(u(x, t), x), x))

 

(D[1](u))(0, t) = 0, (D[1](u))(L, t) = 0, u(x, 0) = -3*cos(8*Pi*x/L)

(9)

`assuming`([pdsolve([pde[5], bc[5]], u(x, t))], [0 < L, 0 < k])

u(x, t) = -3*cos(8*Pi*x/L)*exp(-64*k*Pi^2*t/L^2)

(10)

pde[6] := diff(u(x, t), t) = k*(diff(u(x, t), x, x))+f(x, t); bc[6] := u(0, t) = 0, u(l, t) = 0, u(x, 0) = g(x)

diff(u(x, t), t) = k*(diff(diff(u(x, t), x), x))+f(x, t)

 

u(0, t) = 0, u(l, t) = 0, u(x, 0) = g(x)

(11)

pdsolve([pde[6], bc[6]], u(x, t))

u(x, t) = Sum(2*(Int(g(tau1)*sin(Pi*n1*tau1/l), tau1 = 0 .. l))*sin(Pi*n1*x/l)*exp(-k*Pi^2*n1^2*t/l^2)/l, n1 = 1 .. infinity)+Int(Sum(2*(Int(f(x, tau1)*sin(Pi*n*x/l), x = 0 .. l))*sin(Pi*n*x/l)*exp(-k*Pi^2*n^2*(t-tau1)/l^2)/l, n = 1 .. infinity), tau1 = 0 .. t)

(12)

pde[7] := diff(u(x, t), t) = diff(u(x, t), x, x); bc[7] := u(x, 0) = f(x), u(-1, t) = 0, u(1, t) = 0

diff(u(x, t), t) = diff(diff(u(x, t), x), x)

 

u(x, 0) = f(x), u(-1, t) = 0, u(1, t) = 0

(13)

pdsolve([pde[7], bc[7]], u(x, t))

u(x, t) = Sum((Int(f(x)*sin(n*Pi*x), x = -1 .. 1))*sin(n*Pi*x)*exp(-Pi^2*n^2*t), n = 1 .. infinity)

(14)

pde[8] := diff(u(x, t), t) = k*(diff(u(x, t), x, x))-h*u(x, t); bc[8] := u(x, 0) = sin(x), u(-Pi, t) = u(Pi, t), (D[1](u))(-Pi, t) = (D[1](u))(Pi, t)

diff(u(x, t), t) = k*(diff(diff(u(x, t), x), x))-h*u(x, t)

 

u(x, 0) = sin(x), u(-Pi, t) = u(Pi, t), (D[1](u))(-Pi, t) = (D[1](u))(Pi, t)

(15)

pdsolve([pde[8], bc[8]], u(x, t))

u(x, t) = sin(x)*exp(-t*(k+h))

(16)

pde[9] := diff(u(x, t), t) = diff(u(x, t), x, x)

bc[9] := u(0, t) = 20, u(1, t) = 50, u(x, 0) = 0

diff(u(x, t), t) = diff(diff(u(x, t), x), x)

 

u(0, t) = 20, u(1, t) = 50, u(x, 0) = 0

(17)

pdsolve([pde[9], bc[9]], u(x, t))

u(x, t) = 20+Sum((-40+100*(-1)^n)*sin(n*Pi*x)*exp(-Pi^2*n^2*t)/(Pi*n), n = 1 .. infinity)+30*x

(18)

pde[10] := diff(u(x, y), x, x)+diff(u(x, y), y, y) = 0; bc[10] := u(x, 0) = 0, u(x, 1) = 0, u(0, y) = y^2-y, u(1, y) = 0

diff(diff(u(x, y), x), x)+diff(diff(u(x, y), y), y) = 0

 

u(x, 0) = 0, u(x, 1) = 0, u(0, y) = y^2-y, u(1, y) = 0

(19)

pdsolve([pde[10], bc[10]], u(x, y))

u(x, y) = Sum(-4*((-1)^n-1)*sin(Pi*y*n)*(exp(Pi*n*(2*x-1))-exp(Pi*n))*exp(-Pi*n*(x-1))/((exp(2*Pi*n)-1)*Pi^3*n^3), n = 1 .. infinity)

(20)

pde[11] := diff(u(x, y), x, x)+diff(u(x, y), y, y) = 0

bc[11] := (D[1](u))(0, y) = 0, (D[1](u))(L, y) = 0, u(x, H) = f(x), u(x, 0) = 0

diff(diff(u(x, y), x), x)+diff(diff(u(x, y), y), y) = 0

 

(D[1](u))(0, y) = 0, (D[1](u))(L, y) = 0, u(x, H) = f(x), u(x, 0) = 0

(21)

`assuming`([pdsolve([pde[11], bc[11]], u(x, y))], [0 < L, 0 < H])

u(x, y) = Sum(2*(Int(cos(Pi*x*n/L)*f(x), x = 0 .. L))*cos(Pi*x*n/L)*exp(Pi*n*(H-y)/L)*(exp(2*Pi*y*n/L)-1)/(L*(exp(2*Pi*H*n/L)-1)), n = 1 .. infinity)

(22)

pde[12] := diff(u(x, y), x, x)+diff(u(x, y), y, y) = 0

bc[12] := (D[1](u))(L, y) = 0, u(x, H) = 0, u(x, 0) = 0, (D[1](u))(0, y) = g(y)

diff(diff(u(x, y), x), x)+diff(diff(u(x, y), y), y) = 0

 

(D[1](u))(L, y) = 0, u(x, H) = 0, u(x, 0) = 0, (D[1](u))(0, y) = g(y)

(23)

`assuming`([pdsolve([pde[12], bc[12]], u(x, y))], [0 < x, x <= L, 0 < y, y <= H])

u(x, y) = Sum(-2*(Int(sin(Pi*y*n/H)*g(y), y = 0 .. H))*sin(Pi*y*n/H)*(exp(-Pi*n*(L-2*x)/H)+exp(Pi*L*n/H))*exp(Pi*n*(L-x)/H)/(Pi*n*(exp(2*Pi*L*n/H)-1)), n = 1 .. infinity)

(24)

pde[13] := diff(u(x, y), x, x)+diff(u(x, y), y, y) = 0

bc[13] := (D[1](u))(0, y) = 0, u(x, 0) = 0, u(x, H) = 0, u(L, y) = g(y)

diff(diff(u(x, y), x), x)+diff(diff(u(x, y), y), y) = 0

 

(D[1](u))(0, y) = 0, u(x, 0) = 0, u(x, H) = 0, u(L, y) = g(y)

(25)

`assuming`([pdsolve([pde[13], bc[13]], u(x, y))], [0 < L, 0 < H])

u(x, y) = Sum(2*(Int(sin(Pi*y*n/H)*g(y), y = 0 .. H))*sin(Pi*y*n/H)*exp(Pi*n*(L-x)/H)*(exp(2*Pi*x*n/H)+1)/(H*(exp(2*Pi*L*n/H)+1)), n = 1 .. infinity)

(26)

pde[14] := diff(u(x, y), x, x)+diff(u(x, y), y, y) = 0

bc[14] := u(0, y) = g(y), u(L, y) = 0, (D[2](u))(x, 0) = 0, u(x, H) = 0

diff(diff(u(x, y), x), x)+diff(diff(u(x, y), y), y) = 0

 

u(0, y) = g(y), u(L, y) = 0, (D[2](u))(x, 0) = 0, u(x, H) = 0

(27)

`assuming`([pdsolve([pde[14], bc[14]], u(x, y))], [0 < x, x <= L, 0 < y, y <= H])

u(x, y) = Sum(-2*(exp(-(L-2*x)*(1/2+n)*Pi/H)-exp((1/2)*Pi*(1+2*n)*L/H))*cos((1/2)*Pi*(1+2*n)*y/H)*(Int(cos((1/2)*Pi*(1+2*n)*y/H)*g(y), y = 0 .. H))*exp((1/2)*Pi*(1+2*n)*(L-x)/H)/(H*(exp(Pi*(1+2*n)*L/H)-1)), n = 0 .. infinity)

(28)

pde[15] := diff(u(x, y), x, x)+diff(u(x, y), y, y) = 0

bc[15] := u(0, y) = 0, u(L, y) = 0, u(x, 0) = (D[2](u))(x, 0), u(x, H) = f(x)

diff(diff(u(x, y), x), x)+diff(diff(u(x, y), y), y) = 0

 

u(0, y) = 0, u(L, y) = 0, u(x, 0) = (D[2](u))(x, 0), u(x, H) = f(x)

(29)

`assuming`([pdsolve([pde[15], bc[15]], u(x, y))], [0 < x, x <= L, 0 < y, y <= H])

u(x, y) = Sum(2*exp(Pi*n*(H-y)/L)*sin(Pi*x*n/L)*((Pi*n+L)*exp(2*Pi*y*n/L)+Pi*n-L)*(Int(sin(Pi*x*n/L)*f(x), x = 0 .. L))/(L*((Pi*n+L)*exp(2*Pi*H*n/L)+Pi*n-L)), n = 1 .. infinity)

(30)

pde[16] := diff(u(x, t), t, t) = c^2*(diff(u(x, t), x, x))

bc[16] := u(0, t) = 0, (D[1](u))(L, t) = 0, (D[2](u))(x, 0) = 0, u(x, 0) = f(x)

diff(diff(u(x, t), t), t) = c^2*(diff(diff(u(x, t), x), x))

 

u(0, t) = 0, (D[1](u))(L, t) = 0, (D[2](u))(x, 0) = 0, u(x, 0) = f(x)

(31)

`assuming`([pdsolve([pde[16], bc[16]], u(x, t))], [0 < x, x <= L])

u(x, t) = Sum(2*cos((1/2)*c*Pi*(1+2*n)*t/L)*(Int(sin((1/2)*Pi*(1+2*n)*x/L)*f(x), x = 0 .. L))*sin((1/2)*Pi*(1+2*n)*x/L)/L, n = 0 .. infinity)

(32)

pde[17] := diff(w(x1, x2, x3, t), t) = diff(w(x1, x2, x3, t), x1, x1)+diff(w(x1, x2, x3, t), x2, x2)+diff(w(x1, x2, x3, t), x3, x3)

bc[17] := w(x1, x2, x3, 0) = exp(x1)+x2*x3^5

diff(w(x1, x2, x3, t), t) = diff(diff(w(x1, x2, x3, t), x1), x1)+diff(diff(w(x1, x2, x3, t), x2), x2)+diff(diff(w(x1, x2, x3, t), x3), x3)

 

w(x1, x2, x3, 0) = exp(x1)+x2*x3^5

(33)

pdsolve([pde[17], bc[17]])

w(x1, x2, x3, t) = Sum(t^n*((proc (U) options operator, arrow; diff(diff(U, x1), x1)+diff(diff(U, x2), x2)+diff(diff(U, x3), x3) end proc)@@n)(exp(x1)+x2*x3^5)/factorial(n), n = 0 .. infinity)

(34)

pde[18] := diff(u(x, t), t)+u(x, t)*(diff(u(x, t), x)) = -x

bc[18] := u(x, 0) = x

diff(u(x, t), t)+u(x, t)*(diff(u(x, t), x)) = -x

 

u(x, 0) = x

(35)

pdsolve([pde[18], bc[18]], u(x, t))

u(x, t) = x/tan(t+(1/4)*Pi)

(36)

pde[19] := diff(u(x, t), t)-u(x, t)^2*(diff(u(x, t), x)) = 3*u(x, t)

bc[19] := u(x, 0) = x

diff(u(x, t), t)-u(x, t)^2*(diff(u(x, t), x)) = 3*u(x, t)

 

u(x, 0) = x

(37)

pdsolve([pde[19], bc[19]], u(x, t))

u(x, t) = -exp(3*t)*((-6*exp(6*t)*x+6*x+9)^(1/2)-3)/(exp(6*t)-1)

(38)

pde[20] := diff(u(x, t), t)-x*u(x, t)*(diff(u(x, t), x)) = 0

bc[20] := u(x, 0) = x

diff(u(x, t), t)-x*u(x, t)*(diff(u(x, t), x)) = 0

 

u(x, 0) = x

(39)

pdsolve([pde[20], bc[20]], u(x, t))

u(x, t) = -LambertW(-t*x)/t

(40)

pde[21] := diff(u(x, t), t)+u(x, t)*(diff(u(x, t), x)) = 0

bc[21] := u(x, 0) = x

diff(u(x, t), t)+u(x, t)*(diff(u(x, t), x)) = 0

 

u(x, 0) = x

(41)

pdsolve([pde[21], bc[21]], u(x, t))

u(x, t) = x/(t+1)

(42)

pde[22] := (t+u(x, t))*(diff(u(x, t), x))+t*(diff(u(x, t), t)) = 0; bc[22] := u(x, 1) = x

(t+u(x, t))*(diff(u(x, t), x))+t*(diff(u(x, t), t)) = 0

 

u(x, 1) = x

(43)

pdsolve([pde[22], bc[22]], u(x, t))

u(x, t) = (t-x-1)/(ln(1/t)-1)

(44)

pde[23] := (diff(r*(diff(u(r, theta), r)), r))/r+(diff(u(r, theta), theta, theta))/r^2 = 0; bc[23] := u(a, theta) = f(theta), u(r, -Pi) = u(r, Pi), (D[2](u))(r, -Pi) = (D[2](u))(r, Pi)

(diff(u(r, theta), r)+r*(diff(diff(u(r, theta), r), r)))/r+(diff(diff(u(r, theta), theta), theta))/r^2 = 0

 

u(a, theta) = f(theta), u(r, -Pi) = u(r, Pi), (D[2](u))(r, -Pi) = (D[2](u))(r, Pi)

(45)

`assuming`([pdsolve([pde[23], bc[23]], u(r, theta), HINT = boundedseries)], [a > 0])

u(r, theta) = (1/2)*(2*(Sum(r^n*(sin(n*theta)*(Int(f(theta)*sin(n*theta), theta = -Pi .. Pi))+cos(n*theta)*(Int(f(theta)*cos(n*theta), theta = -Pi .. Pi)))*a^(-n)/Pi, n = 1 .. infinity))*Pi+Int(f(theta), theta = -Pi .. Pi))/Pi

(46)

pde[24] := diff(g(t, x), t) = diff(g(t, x), x, x)+a*g(t, x); bc[24] := g(0, x) = f(x)

diff(g(t, x), t) = diff(diff(g(t, x), x), x)+a*g(t, x)

 

g(0, x) = f(x)

(47)

pdsolve([pde[24], bc[24]])

g(t, x) = exp(a*t)*invfourier(fourier(f(x), x, s1)*exp(-t*s1^2), s1, x)

(48)

pde[25] := diff(u(x, y), x, x)+diff(u(x, y), y, y) = 0; bc[25] := u(0, y) = y*(-y+1), u(1, y) = 0, (D[2](u))(x, 0) = 0, (D[2](u))(x, 1) = 0

diff(diff(u(x, y), x), x)+diff(diff(u(x, y), y), y) = 0

 

u(0, y) = y*(-y+1), u(1, y) = 0, (D[2](u))(x, 0) = 0, (D[2](u))(x, 1) = 0

(49)

pdsolve([pde[25], bc[25]], u(x, y))

u(x, y) = Sum(2*((-1)^n+1)*cos(Pi*y*n)*(exp(Pi*n*(2*x-1))-exp(Pi*n))*exp(-Pi*n*(x-1))/((exp(2*Pi*n)-1)*Pi^2*n^2), n = 1 .. infinity)

(50)

``


 

Download ImprovementsInPdsolve.mw

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

And the Nobel prize in physics 2017 went for work in General Relativity ! Actually, experimental work involving sophisticated detectors and Numerical Relativity, one of the branches of GR. The prize was awarded to Rainer Weiss (85 years old, 1/2 of the prize), Barry Barish (81 years old, 1/4 of the prize) and Kip Thorne (77 years old, 1/4 of the prize) who have "shaken the world again" with their work on Ligo experiment, which was able to detect ripples in the fabric of spacetime.

General Relativity continues to be at the center of work in theoretical and experimental physics. I take this opportunity to note that, in Maple 2017, among the several improvements in the Physics package regarding General Relativity, there is a new package, Physics:-ThreePlusOne, all dedicated to the symbolic manipulations necessary to formulate problems in Numerical Relativity.

The GR functionality implemented in Physics, Physics:-Tetrads and Physics:-ThreePlusOne is unique in computer algebra systems and reflects the Maplesoft intention, for several years now, to provide the very best possible computer algebra environment for Physics, regarding current research activity as well as related education in advanced mathematical-physics methods.

For what is going on in theoretical physics nowadays and its connection with General Relativity, in very short: the unification of gravity with the other forces, check for instance this map from Aug/2015 (by the way a very nice summary for whoever is interested):

It is also interesting the article behind this map of topics as well as this brilliant and accessible presentation by Nima Arkani-Hamed (Princeton):  Quantum Mechanics and Spacetime in the 21st Century, given in the Perimeter Institute for Theoretical Physics (Waterloo), November 2014. 

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

Group of exercises solved using Maple scientific software, with the necessary considerations of some basic commands: evalf and convert that will show the solutions with the user-defined digits and the angular measurement in sexagesimal degrees. Important use of the law of the triangle through of vector position applied to vectors in vector spaces, vector force and vector moment for engineering students. In spanish.

Exercises_of_vectors_forces_and_moment_with_Maple.mw

Videotutorial:

https://www.youtube.com/watch?v=DxpO0gc5GCA

Lenin Araujo Castillo

Ambassador of Maple

The development of the calculation of moments using force vectors is clearly observed by taking a point and also a line. Different exercises are solved with the help of Maple syntax. We can also visualize the vector behavior in the different configurations of the position vector. Applications designed exclusively for engineering students. In Spanish.

Moment_of_a_force_using_vectors_updated.mw

Lenin Araujo Castillo

Ambassador of Maple

A project that I have been working on is adding some functionality for Cluster Analysis to Maple (a small part of a much bigger project to increase Maple’s toolkit for exploratory data mining and data analysis). The launch of the MapleCloud package manager gave me a way to share my code for the project as it evolves, providing others with some useful new tools and hopefully gathering feedback (and collaborators) along the way.

At this point, there aren’t a lot of commands in the ClusterAnalysis package, but I have already hit upon several interesting applications. For example, while working on a command for plotting clusters of points, one problem I encountered was how to draw the minimal volume enclosing ellipsoid around a group (or cluster) of points. After doing some research, I stumbled upon Khachiyan’s Algorithm, which related to solving linear programming problems with rational data. The math behind this is definitely interesting, but I’m not going to spend any time on it here. For further reading, you can explore the following:

Khachiyan’s Algorithm had previously been applied in some other languages, but to the best of my knowledge, did not have any Maple implementations. As such, the following code is an implementation of Khachiyan’s Algorithm in 2-D, which could be extended to N-dimensional space rather easily.

This routine accepts an Nx2 dataset and outputs either a plot of the minimum volume enclosing ellipsoid (MVEE) or a list of results as described in the details for the ‘output’ option below.

MVEE( X :: DataSet, optional arguments, additional arguments passed to the plotting command );

The optional arguments are as follows:

  • tolerance : realcons;  specifies the convergence criterion
  • maxiterations : posint; specifies the maximum number of iterations
  • output : {identical(data,plot),list(identical(data,plot))}; specifies the output. If output includes plot, then a plot of the enclosing ellipsoid is returned. If output includes data, then the return includes is a list containing the matrix A, which defines the ellipsoid, the center of the ellipse, and the eigenvalues and eigenvectors that can be used to find the semi-axis coordinates and the angle of rotation, alpha, for the ellipse.
  • filled : truefalse; specifies if the returned plot should be filled or not

Code:

#Minimum Volume Enclosing Ellipsoid
MVEE := proc(XY, 
              {tolerance::positive:= 1e-4}, #Convergence Criterion
              {maxiterations::posint := 100},
              {output::{identical(data,plot),list(identical(data,plot))} := data},
              {filled::truefalse := false} 
            )

    local alpha, evalues, evectors, i, l_error, ldata, ldataext, M, maxvalindex, n, ncols, nrows, p1, semiaxes, stepsize, U, U1, x, X, y;
    local A, center, l_output; #Output

    if hastype(output, 'list') then
        l_output := output;
    else
        l_output := [output];
    end if;

    kernelopts(opaquemodules=false):

    ldata := Statistics:-PreProcessData(XY, 2, 'copy');

    nrows, ncols := upperbound(ldata);
    ldataext := Matrix([ldata, Vector[column](nrows, ':-fill' = 1)], 'datatype = float');

    if ncols <> 2 then
        error "expected 2 columns of data, got %1", ncols;
    end if;

    l_error := 1;

    U := Vector[column](1..nrows, 'fill' = 1/nrows);

    ##Khachiyan Algorithm##
    for n to maxiterations while l_error >= tolerance do

        X := LinearAlgebra:-Transpose(ldataext) . LinearAlgebra:-DiagonalMatrix(U) . ldataext;
        M := LinearAlgebra:-Diagonal(ldataext . LinearAlgebra:-MatrixInverse(X) . LinearAlgebra:-Transpose(ldataext));
        maxvalindex := max[index](map['evalhf', 'inplace'](abs, M));
        stepsize := (M[maxvalindex] - ncols - 1)/((ncols + 1) * (M[maxvalindex] - 1));
        U1 := (1 - stepsize) * U;
        U1[maxvalindex] := U1[maxvalindex] + stepsize;
        l_error := LinearAlgebra:-Norm(LinearAlgebra:-DiagonalMatrix(U1 - U));
        U := U1;

    end do;

    A := (1/ncols) * LinearAlgebra:-MatrixInverse(LinearAlgebra:-Transpose(ldata) . LinearAlgebra:-DiagonalMatrix(U) . ldata - (LinearAlgebra:-Transpose(ldata) . U) . LinearAlgebra:-Transpose((LinearAlgebra:-Transpose(ldata) . U)));
    center := LinearAlgebra:-Transpose(ldata) . U;
    evalues, evectors := LinearAlgebra:-Eigenvectors(A);
    evectors := evectors(.., sort[index](1 /~ (sqrt~(Re~(evalues))), `>`, ':-output' = ':-permutation'));
    semiaxes := sort(1 /~ (sqrt~(Re~(evalues))), `>`);
    alpha := arctan(Re(evectors[2,1]) / Re(evectors[1,1]));

    if l_output = [':-data'] then
        return A, center, evectors, evalues;
    elif has( l_output, ':-plot' ) then
            x := t -> center[1] + semiaxes[1] * cos(t) * cos(alpha) - semiaxes[2] * sin(t) * sin(alpha);
            y := t -> center[2] + semiaxes[1] * cos(t) * sin(alpha) + semiaxes[2] * sin(t) * cos(alpha);
            if filled then
                p1 := plots:-display(subs(CURVES=POLYGONS, plot([x(t), y(t), t = 0..2*Pi], ':-transparency' = 0.95, _rest)));
            else
                p1 := plot([x(t), y(t), t = 0..2*Pi], _rest);
            end if;
        return p1, `if`( has(l_output, ':-data'), op([A, center, evectors, evalues]), NULL );
    end if;

end proc:

 

You can run this as follows:

M:=Matrix(10,2,rand(0..3)):

plots:-display([MVEE(M,output=plot,filled,transparency=.3),
                plots:-pointplot(M, symbol=solidcircle,symbolsize=15)],
size=[0.5,"golden"]);

 

 

As it stands, this is not an export from the “work in progress” ClusterAnalysis package – it’s actually just a local procedure used by the ClusterPlot command. However, it seemed like an interesting enough application that it deserved its own post (and potentially even some consideration for inclusion in some future more geometry-specific package). Here’s an example of how this routine is used from ClusterAnalysis:

with(ClusterAnalysis);

X := Import(FileTools:-JoinPath(["datasets/iris.csv"], base = datadir));

kmeans_results := KMeans(X[[`Sepal Length`, `Sepal Width`]],
    clusters = 3, epsilon = 1.*10^(-7), initializationmethod = Forgy);

ClusterPlot(kmeans_results, style = ellipse);

 

 

The source code for this is stored on GitHub, here:

https://github.com/dskoog/Maple-ClusterAnalysis/blob/master/src/MVEE.mm

Comments and suggestions are welcomed.

 

If you don’t have a copy of the ClusterAnalysis package, you can install it from the MapleCloud window, or by running:

PackageTools:-Install(5629844458045440);

 

These worksheets provide the volume calculations  of a small causal diamond near the tip of the past light cone, using dimensional analysis and particular test metrics.

I recommend them for anyone working in causet theory on the problem of finding higher order corrections.

2D.mw

4D.mw

4Dflat.mw

 

 

 

 

With this app you will be able to interpret the curvatures generated by two position vectors, either in the plane or in space. Just enter the position vectors and drag the slider to calculate the curvature at different times and you will of course be able to observe its respective graph. At first I show you how it is developed using the natural syntax of Maple and then optimize our
 app with the use of buttons. App made in Maple for engineering students. In spanish.

Plot_of_Curvature.mw

Videotutorial:

https://www.youtube.com/watch?v=SbXFgr_5JDE

Lenin Araujo Castillo

Ambassador of Maple

First 29 30 31 32 33 34 35 Last Page 31 of 77