Applications, Examples and Libraries

Share your work here

Let us consider the linear integer programming problem:

A := Matrix([[1, 7, 1, 3], [1, 6, 4, 6], [17, 1, 5, 1], [1, 6, 10, 4]]):
 n := 4; z := add(add(A[i, j]*x[i, j], j = 1 .. n), i = 1 .. n):
restr := {seq(add(x[i, j], i = 1 .. n) = 1, j = 1 .. n), seq(add(x[i, j], j = 1 .. n) = 1, i = 1 .. n)}:
 sol := Optimization[LPSolve](z, restr, assume = binary);

Error, (in Optimization:-LPSolve) no feasible integer point found; 
use feasibilitytolerance option to adjust tolerance

sol1 := Optimization[LPSolve](z, restr, assume = binary, feasibilitytolerance = 100, integertolerance = 1);

Error, (in Optimization:-LPSolve) no feasible integer point found;
 use feasibilitytolerance option to adjust tolerance

That was OK in Maple 16, outputting

.

The bug in one of the principal Maple commands lasts since Maple 2015, where the above code causes "Kernel connection has been lost". The SCRs about it were submitted three times (see http://www.mapleprimes.com/questions/204750-Bug-In-LPSolve-In-Maple-20151).

The Möbius strip  Mobius_strip_rolling.mw

Variants :


The line and the curve on the surface.

 

Recently, I came across an addendum to a problem that appears in many calculus texts, an addendum I had never explored. It intrigued me, and I hope it will capture your attention too.

The problem is that of girding the equator of the earth with a belt, then extending by one unit (here, taken as the foot) the radius of the circle so formed. The question is by how much does the circumference of the belt increase. This problem usually appears in the section of the calculus text dealing with linear approximations by the differential. It turns out that the circumference of the enlarged band is 2*Pi ft greater than the original band.

(An alternate version of this has the circumference of the band increased by one foot, with the radius then being increased by 0.16 ft.)

The addendum to the problem then asked how high would the enlarged band be over the surface of the earth if it were lifted at one point and drawn as tight as possible around the equator. At first, I didn't know what to think. Would the height be some surprisingly large number? And how would one go about calculating this height.

It turns out that the enlarged and lifted band would be some 616.67 feet above the surface of the earth! This is significantly larger than the increase in the diameter of the original band. So, the result is a surprise, at least to me.

This is the kind of amusement that retirement affords. I heartily recommend both the amusement and the retirement. The supporting calculations can be found in the attached worksheet: Girding.mw

Let us consider 

restart; 
MultiSeries:-limit(sin(n)/n, n = infinity, complex);
0

The answer is wrong: in view of the Casorati-Weierstrass theorem the limit does not exist. Let us try another limit command of Maple

limit(sin(n)/n, n = infinity, complex);


(lim) (sin(n))/(n)

which fails. Therefore, Maple user does not obtain the correct answer. 

Suppose we have some simple animations. Our goal - to build a more complex animation, combining the original animations in different ways.
We show how to do it on the example of the three animations. The technique is general and can be applied to any number of animations.

Here are the three simple animations:

restart;
with(plots):
A:=animate(plot, [sin(x), x=-Pi..a, color=red, thickness=3], a=-Pi..Pi):
B:=animate(plot, [x^2-1, x=-2..a, thickness=3, color=green], a=-2..2): 
C:=animate(plot, [[4*cos(t),4*sin(t), t=0..a], color=blue, thickness=3], a=0..2*Pi):

 

In Example 1 all three animation executed simultaneously:

display([A, B, C], view=[-4..4,-4..4]);

                                

 

In Example 2, the same animation performed sequentially. Note that the previous animation disappears completely when the next one begins to execute:

display([A, B, C], insequence);

                                 

 

Below we show how to save the last frame of every previous animation into subsequent animations:

display([A, display(op([1,-1,1],A),B), display(op([1,-1,1],A),op([1,-1,1],B),C)], insequence);

                                 

 

Using this technique, we can anyhow combine the original animations. For example, in the following example at firstly animations   and  B  are executed simultaneously, afterwards C is executed:

display([display(A, B), display(op([1,-1,1],A),op([1,-1,1],B),C)], insequence);

                                     

 

The last example in 3D I have taken from here:

restart;
with(plots):
A:=animate(plot3d,[[2*cos(phi),2*sin(phi),z], z =0..a, phi=0..2*Pi, style=surface, color=red], a=0..5):
B:=animate(plot3d,[[(2+6/5*(z-5))*cos(phi), (2+6/5*(z-5))*sin(phi),z], z=5..a, phi=0..2*Pi, style=surface, color=blue], a=5..10):
C:=animate(plot3d,[[8*cos(phi),8*sin(phi),z], z =10..a, phi=0..2*Pi, style=surface, color=green], a=10..20):
display([A, display(op([1,-1,1],A),B), display(op([1,-1,1],A),op([1,-1,1],B),C)], insequence, scaling=constrained, axes=normal);

                        


 

AA.mw

I'd like to pay attention of Maple community to the recent work by Alex Degtyarev in algebraic geometry done with Maple.

Bertini.zip

Let us consider 

with(Statistics);
U := RandomVariable(DiscreteUniform(-10, 10)):
V := RandomVariable(DiscreteUniform(-10, 10)):
Probability(U^2-V^2 <= 1/9, numeric);
  0.

, whereas a positive number greater than 1/21 is expected. 

 

Let us consider the example from Maple help to ?ProbabilityFunction (also see ?Geometric)

with(Statistics):
ProbabilityFunction(Geometric(1/3), 5);
                              32 /729
                             

Let us continue the investigation

ProbabilityFunction(Geometric(1/3), 5.1);
0.4215152817e-1
ProbabilityFunction(Geometric(1/3), 5.12);
0.4181109090e-1
ProbabilityFunction(Geometric(1/3), 51/10)
(32/2187)*2^(1/10)*3^(9/10)

whereas the result 0 is expected in all the three cases up to Wiki. I am aware of the line

"t-algebraic; point (assumed to be an integer)"

in the help. However, 

ProbabilityFunction(Geometric(1/3), -.5);
                               0

The same issue with the DiscreteUniform distribution. This bug lasts from  at least Maple 16. The question arises: may we trust Maple?

Everything is simple, until you go underwater – This is what the University of Waterloo Submarine Racing team, or in short ‘WatSub’ coined as their motto. Never mind learning to scuba dive, and dealing with such things as rust, this newly formed team would have to compete against university teams with a decade or more of experience.

But that did not deter the team, and they started work on Ontario’s first submarine racing project. The team approached Maplesoft to be a sponsor and we are proud to have supported this ingenious venture. The team has used Maplesoft technology in the design and testing of the submarine.

“Maple has been our go-to calculations and analysis tool throughout the development of Amy (2015-2016 season), and we will continue using it throughout the development of Bolt (2016-2017 season),” said Gonzalo Espinoza Graham, President of the WatSub Team. “Its familiar interface and computing environment allowed us to set design benchmark targets from early on the design process and follow through with them on the later stage.”

What started as an engineering project in December 2014, becoming officially the first submarine racing team in Ontario. The team soon grew to over 130 general members and a tight core-team, who were eager to tackle new challenges.  The team resides inside the Sedra Student Design Centre, University of Waterloo’s state of the art facility that houses over 25 student teams, the largest of its kind in North America.  

WatSub made its first appearance on the European International Submarine Races (eISR) back in July 2016, with its 1st submarine ‘Amy’, where a single scuba diver piloted the submarine and propelled it through an unforgiving winding course marked by obstacles and turns 10 meters underwater. The team has since then participated in other competitions and is constantly improving the design and performance of the submarine, learning from each competition they participate in.  Next year Amy will participate in the 14th edition of the eISR international competition. “I think the greatest thing we learned is never to give up,” said Ana Krstanovic, a third-year political science student who manages communications for the team. “We’re more motivated now than ever.”

 

Ojaswi Tagore, Gonzalo Espinoza Graham, and Janna Henzl represented WatSub at the European International Submarine Race in Gosport, UK.

 

Another example of an innovative project that Maplesoft supported in 2016 is Waterloop: The Canadian SpaceX Hyperloop Competition Team, Canada's only SpaceX Hyperloop Pod Competition team. This project, which could change the way we travel in the future, is driven by a group of dedicated University of Waterloo students who have taken on the challenge to design and build a functional prototype Hyperloop pod. They will test it on a one-mile test track in Hawthorne, California in January 2017, pitting it against 22 of the 1200+ teams who originally entered the competition.

The Hyperloop is a conceptual next generation high-speed transit system that will take commuters between cities at speeds over 1,000 km/h. The technology will differ from previous rail transit by having pods ride on a cushion of air in a reduced pressure tube in order to reach greater speeds with a smoother ride, and is powered entirely by renewable energy.

 The Hyperloop Pod Competition was launched by Elon Musk, the billionaire engineer and founder of SpaceX and Tesla Motors.  The competition is separated into 3 rounds. The first one was held in late December, where selected teams sent in their initial designs to be reviewed. From there, 180 teams were chosen to compete at Texas A&M University. Each team set up a booth and a panel of judges critiqued them and chose 31 teams to move onto the final, build and test stage.

Waterloop Goose I

Waterloop Goose X

The GOOSE I is Waterloop’s half-scale, functional prototype vehicle pod, which will be the one in the competition.  The GOOSE X pod is a conceptual full size Hyperloop vehicle inspired by the prototype they are building. The full size pod will have a capacity of 26 passengers per pod.

"Our prototype has been designed to be as simple and economical as possible, while still performing all necessary functions for the full size Hyperloop. If it is successful, it has the potential to revolutionize the transit industry in the same manner the train and airplane has before it," said Montgomery de Luna, architectural design lead for Waterloop. “We would like to thank Maplesoft for their generous support.  Without sponsors like Maplesoft supporting our vision and encouraging innovative student projects, we wouldn’t be able to achieve our goal.”

Revolutionizing the transportation industry isn’t easy and is at times frustrating and time consuming for these teams, but having the best tools and resources will ensure that the teams have a good chance at excelling in competitions and creating innovative models that could change our future.

   

 

The code for the animation:

L:=[[-0.12,2],[-0.14,0],[0.14,0],[0.12,2]]:
L1:=[[0.05,2],[4,1],[2,4],[3.5,3.5],[1,7],[2,6.5],[0,10]]:
A:=plot(L, color=brown, thickness=10):
B:=plot([op(L1),op(map(t->[-t[1],t[2]],ListTools:-Reverse(L1)))], color="Green", thickness=10):
C:=plottools:-polygon([op(L1),op(map(t->[-t[1],t[2]],ListTools:-Reverse(L1)))], color=green):
Tree:=plots:-display([A, B, C], scaling=constrained, axes=none):
T:=[[-3.2,-2, Happy, color=blue, font=[times,bold,30]], [0,-2,New, color=blue, font=[times,bold,30]], [2.5,-2,Year, color=blue, font=[times,bold,30]], [-5,-3.5, "&", color=yellow, font=[times,bold,30]],[-2.5,-3.5, Merry, color=red, font=[times,bold,30]], [2.3,-3.5, Christmas!, color=red, font=[times,bold,30]], [0,-5, "2017", color=cyan, font=[times,bold,36]]$5]:
F:=k->plottools:-homothety(Tree, k, [0,5]):
A:=plots:-animate(plots:-display, ['F'(k)], k=0..1, frames=60, paraminfo=false):
B:=plots:-animate(plots:-textplot,[T[1..round(i)]], i=0..nops(T), frames=60, paraminfo=false):
plots:-display(A, B, size=[500,550], scaling=constrained);


Christmas_Tree.mw

 Edit.

 

Parametric equation of second-order curve in 3d. Draghilev method.
PLAN_CURVE_3d_1.mw
Examples:
x1^2+x1*x3+13*x2^2+x3-1=0;
x1+x2+x3=0;


 x1^2+0.1*x2^2+x3^2-9=0;
 x1+3*x3+1=0;


 x1^2-0.1*x2^2+x3^2-9=0;
 x1+3*x3+1=0;

Parametric equation of a circle in 3d by three points. Draghilev method.

CIRCLE_3_POINTS_geom3d_2.mw

In this post I want to present an easy method to obtain a discrete parametrization of a surface S defined implicitly (f(x,y,z)=0).
This problem was discussed here several times, the most recent post is
http://www.mapleprimes.com/posts/207661-Isolation-Of-Sides-Of-The-Surface-On-The-Graph

S is supposed to be the boundary of a convex body having (x0,y0,z0) an interior point and contained in a ball of radius R centered at (x0,y0,z0).
Actually, the procedure also works if the body is only star-shaped with respect to the interior point, and it is also possible to plot only a part of the surface
inside a solid angle centered at (x0,y0,z0).

Usage:
Par3d(f, x=x0, y=y0, z=z0, R, m, n,  theta1 .. theta2,  phi1 .. phi2)

f           is an expression depending on the variables x, y, z
x0, y0, z0  are the coordinates of the interior point
R           is the radius of the ball which contains the surface,
m, n        are the numbers of the grid lines which will be generated
The last two parameters are optional and are used when only a part of S will be parametrized.

The procedure Par3d returns a MESH structure M, which can be plotted with PLOT3D(M).

Par3d :=proc(f,x::`=`,y::`=`,z::`=`,R,m,n,th:=0..2*Pi,ph:=0..Pi)
    local A,i,j, rij,fij,Cth,Sth,Cph,Sph, theta,phi, r;
    A:=Array(1..m+1,1..n+1,1..3,datatype=float[8]);
    for i from 0 to m do for j from 0 to n do
      theta:=op(1,th)+i/m*(op(2,th)-op(1,th));
      phi:=op(1,ph)+j/n*(op(2,ph)-op(1,ph));
      Cth:=evalf(cos(theta)); Sth:=evalf(sin(theta));
      Cph:=evalf(cos(phi));   Sph:=evalf(sin(phi));
      fij:= eval(f, [lhs(x)=rhs(x)+r*Sph*Cth, lhs(y)=rhs(y)+r*Sph*Sth, lhs(z)=rhs(z)+r*Cph]);
      rij:=fsolve(fij,r=0..R);  if [rij]::list(numeric) then rij:=min(rij) fi; 
      if [rij]=[] or not(type(rij,numeric)) then print(['i'=i,'j'=j], fij); rij:=undefined fi; 
      A[i+1,j+1,1]:=evalf(rhs(x)+rij*Sph*Cth);
      A[i+1,j+1,2]:=evalf(rhs(y)+rij*Sph*Sth);
      A[i+1,j+1,3]:=evalf(rhs(z)+rij*Cph);
    od;od:
    MESH(A);
end:

The procedure is not optimized, e.g.
- Cth, etc could be Vectors computed outside the loops
- Some small changes to use evalhf.

###### EXAMPLES ######

f1 := x^2+3*y^2+4*z^2 - x*y - 2*y*z - 10:
plots:-implicitplot3d(f1, x=-5..5, y=-5..5, z=-2..2);

M:=Par3d(f1, x=0,y=0,z=0,5,40,40):
PLOT3D(M);

f2 := x^4+y^4+z^4-1:
M:=Par3d(f2, x=0,y=0,z=0,5,40,40):
PLOT3D(M);

M:=Par3d(f2, x=0,y=0,z=0, 5,40,40, 0..Pi, 0 .. Pi/3): #Plot half of the top only
plots:-display(PLOT3D(M), scaling=constrained);

M:=Par3d(f2,      x=0,y=0,z=0, 5,30,30, 0..Pi, 0 .. Pi):
N:=Par3d(f2+0.01, x=0,y=0,z=0, 5,30,30, 0..Pi, 0 .. Pi):
plots:-display(PLOT3D(M), color=red):
plots:-display(PLOT3D(N), color=green):
plots:-display(%,%%, orientation=[-40,65,10]);

 

f3 := (x^2+y^2-1)^2+(z+sin(x*y+z))^4-120:
plots:-implicitplot3d(f3, x=-4..4,y=-4..4,z=-5..5, numpoints=10000);

Par3d(f3, x=0,y=0,z=0,5, 30,30):
PLOT3D(%);

Note.
The procedure could be used to plot locally around a point (x0,y0,z0)
One may use the spherical coordinates (theta0,phi0) and then call the procedure taking theta0-a .. theta0+a,  phi0-b, .. phi0+b  for the trailing parameters
The spherical coordonates can be computed using:

ThetaPhi :=proc(x,y,z, X,Y,Z)
    local r:=sqrt((X-x)^2+(Y-y)^2+(Z-z)^2);
    ['theta'=arctan(Y-y,X-x), 'phi'=arccos((Z-z)/r)]
end:

ThetaPhi(10,20,30, 11,21,28);evalf(%);

 

 

One way is coloring a surface on both sides. We build equidistant surface with very small radius and stain the equidistant surface in color different from the color of the original surface.
Examples coloring of surfaces on both sides.  Radius equal to abs (0.0001).
x3-0.5*exp(sin(x1+2.5*x2+x3))=0;
(x1^2+x2^2-0.4)^2+(x3+sin(x1*x2+x3))^4-0.1=0;

2_COLORS.mw


Let us consider

restart; Digits := 20; evalf(Int(abs(cos(1/t)), t = 0 .. 0.1e-1), 3);
   -0.639e-2

Pay your attention to the minus sign. Simply no words. Mma produces 0.006377.

evalf@Int.mw

First 28 29 30 31 32 33 34 Last Page 30 of 71