## TechnicalSupport

Technical Support

## 645 Reputation

15 years, 132 days
Maplesoft

## Optimizing Mini Golf in MapleSim...

Playing mini-golf recently, I realized that my protractor can only help me so far since it can't calculate the speed of the swing needed.  I decided a more sophisticated tool was needed and modeled a trick-shot in MapleSim.

To start, I laid out the obstacles, the ball and club, the ground, and some additional visualizations in the MapleSim environment.

When running the simulation, my first result wasn't even close to the hole (similar to when I play in real life!).

The model clearly needed to be optimized. I went to the Optimization app in MapleSim (this can be found under Add Apps or Templates  on the left hand side).

Inside the app I clicked "Load System" then selected the parameters I wanted to optimize.

For this case, I'm optimizing 's' (the speed of the club) and 'theta' (the angle of the club). For the Objective Function I added a Relative Translation Sensor to the model and attached a probe to the Vector Norm of the output.

Inside the app, I switched to the Objective Function section.  Selecting Probes, I added the new probe as the Objective Function by giving it a weight of 1.

Scrolling down to "Execute Parameter Optimization", I checked the "Use Global Optimization Toolbox" checkbox, and clicked Run Parameter Optimization.

Following a run time of 120 seconds, the app returns the graph of the objective function.

Below the plot, optimal values for the parameters are given. Plugging these back into the parameter block for the simulation we see that the ball does in fact go into the hole. Success!

Mini_golf_Global_Optimization.msim

## Gingerbread House with 3D Geometry...

Maple 2019

We recently had a question about using some of the plotting commands in Maple to draw things. We were feeling creative and thought why not take it a step further and draw something in 3D.

Using the geom3d, plottools, and plots packages we decided to make a gingerbread house.

To make the base of the house we decided to use 2 cubes, as these would give us additional lines and segments for the icing on the house.

```point(p__1,[2,3,2]):
point(p__2,[3,3,2]):
cube(c1,p__1,2):
cube(c2,p__2,2):
base:=draw([c1,c2],color=tan);```

Using the same cubes but changing the style to be wireframe and point we made some icing lines and decorations for the gingerbread house.

```base_decor1:=draw([c1,c2],style=wireframe,thickness=3,color=red,transparency=0.2):
base_decor2:=draw([c1,c2],style=wireframe,thickness=10,color=green,linestyle=dot):
base_decor3:=draw([c1,c2],style=point,thickness=2,color="Silver",symbol=sphere):
base_decor:=display(base_decor1,base_decor2,base_decor3);```

To create the roof we found the vertices of the cubes and used those to find the top corners of the base.

```v1:=vertices(c1):
v2:=vertices(c2):
pc1:=seq(point(pc1||i,v1[i]),i=1..nops(v1)):
pc2:=seq(point(pc2||i,v2[i]),i=1..nops(v2)):
topCorners:=[pc1[5],pc1[6],pc2[1],pc2[2]]:
d1:=draw(topCorners):```

Using these top corners we found the midpoints (where the peak of the roof would be) and added the roof height to the coordinates.

```midpoint(lc1,topCorners[1],topCorners[2]):
detail(lc1);```

```point(cc1,[-(2*sqrt(3))/3 + 2, (2*sqrt(3))/3 + 3+1, 2]):
d3:=draw(cc1):

midpoint(lc2,topCorners[3],topCorners[4]):
detail(lc2);```

```point(cc2,[(2*sqrt(3))/3 + 3, (2*sqrt(3))/3 + 3+1, 2]):
d4:=draw(cc2):```

With the midpoints and vertices at the front and rear of the house we made two triangles for the attic of the gingerbread house.

```triangle(tf,[topCorners[1],topCorners[2],cc1]):
front:=draw(tf,color=brown):

triangle(tb,[topCorners[3],topCorners[4],cc2]):
back:=draw(tb,color=tan):```

Using these same points again we made more triangles to be the roof.

```triangle(trl,[cc1,cc2,pc1[5]]):
triangle(trh,[pc2[2],pc1[6],cc1]):
triangle(tll,[cc1,cc2,pc2[2]]):
triangle(tlh,[pc2[1],pc1[5],cc2]):
roof:=draw([trl,trh,tll,tlh],color="Chocolate");```

Our gingerbread house now had four walls, a roof, and icing, but no door. Creating the door was as easy as making a parallelepiped, but what is a door without more icing?

```door:=display(plottools:-parallelepiped([1,0,0],[0,1.2,0],[0,0,0.8],[0.8,1.9,1.6]),color="DarkRed"):
door_decor1:=display(plottools:-parallelepiped([1,0,0],[0,1.2,0],[0,0,0.8],[0.8,1.9,1.6]),color="Gold",style=line):
door_decor2:=display(plottools:-parallelepiped([1,0,0],[0,1.2,0],[0,0,0.8],[0.8,1.9,1.6]),color="Silver", style=line,linestyle=dot,thickness=5):
door_decor:=display(door_decor1,door_decor2):```

Now having a door we could have left it like this, but what better way to decorate a gingerbread house than with candy canes? Naturally, if we were going to have one candy cane we were going to have lots of candy canes. To facilitate this we made a candy cane procedure.

```candy_pole:=proc(c:=[0,0,0], {segR:=0.1}, {segH:=0.1}, {segn:=7}, {tilt_theta:=0}, {theta:=0}, {arch:=true}, {flip:=false})
local cane1,cane2,cane_s,cane_c,cane0,cane,i,cl,cd,ch, cane_a,tmp,cane_ac,cane_a1,cane00,cane01,cane02,cane_a1s,tmp2,cane_a2s:
uses plots,geom3d:
cl:=c[1]:
cd:=c[2]:
ch:=c[3]:
cane1:=plottools:-cylinder([cd, ch, cl], segR, segH,style=surface):

cane2:=display(plottools:-rotate(cane1,Pi/2,[[cd,ch,cl],[cd+1,ch,cl]])):
cane_s:=[cane2,seq(display(plottools:-translate(cane2,0,segH*i,0)),i=1..segn-1)]:
cane_c:=seq(ifelse(type(i,odd),red,white),i=1..segn):

cane0:=display(cane_s,color=[cane_c]):

if arch then

cane_a:=plottools:-translate(cane2,0,segH*segn-segH/2,0):
tmp:=i->plottools:-rotate(cane_a,i*Pi/24, [ [cd,ch+segH*segn-segH/2,cl+segR*2] , [cd+1,ch+segH*segn-segH/2,cl+segR*2] ] ):

cane_ac:=seq(ifelse(type(i,odd),red,white),i=1..24):

cane_a1s:=seq(plottools:-translate(tmp(i),0,segH*i/12,segR*i/4),i=1..12):

tmp2:=i->plottools:-rotate(cane_a1s[12],i*Pi/24,[[cd,ch+segH*segn-0.05,cl+segR*2],[cd+1,ch+segH*segn-0.05,cl+segR*2]]):

cane_a2s:=seq(plottools:-translate(tmp2(i),0,-segH*i/500,segR*i/4),i=1..12):
cane_a1:=display(cane_a1s,cane_a2s,color=[cane_ac]):
cane00:=display(cane0,cane_a1);

if flip then

cane01:=plottools:-rotate(cane00,tilt_theta,[[cd,ch,cl],[cd+1,ch,cl]]):
cane02:=plottools:-rotate(cane01,theta,[[cd,ch,cl],[cd,ch+1,cl]]):
cane:=plottools:-reflect(cane01,[[cd,ch,cl],[cd,ch+1,cl]]):

else

cane01:=plottools:-rotate(cane00,tilt_theta,[[cd,ch,cl],[cd+1,ch,cl]]):
cane:=plottools:-rotate(cane01,theta,[[cd,ch,cl],[cd,ch+1,cl]]):

end if:

return cane:

else

cane:=plottools:-rotate(cane0,tilt_theta,[[cd,ch,cl],[cd+1,ch,cl]]):

return cane:

end if:

end proc:```

With this procedure we decided to add candy canes to the front, back, and sides of the gingerbread house. In addition we added two candy poles.

Candy Canes in front of the house:

```cane1:=candy_pole([1.2,0,2],segn=9,arch=false):
cane2:=candy_pole([2.8,0,2],segn=9,arch=false):
cane3:=candy_pole([2.7,0.8,3.3],segn=9,segR=0.05,tilt_theta=-Pi/7):
cane4:=candy_pole([1.3,0.8,3.3],segn=9,segR=0.05,tilt_theta=-Pi/7,flip=true):
front_canes:=display(cane1,cane2,cane3,cane4):```

Candy Canes at the back of the house:

```caneb3:=candy_pole([1.5,4.2,2.5],segn=15,segR=0.05,tilt_theta=-Pi/3,flip=true):
caneb4:=candy_pole([2.5,4.2,2.5],segn=15,segR=0.05,tilt_theta=-Pi/3):}
back_canes:=display(caneb3,caneb4):```

Candy Canes at the left of the house:

```canel1:=candy_pole([0.8,1.5,2.5],segn=15,segR=0.05,tilt_theta=-Pi/7,theta=Pi/2):
canel2:=candy_pole([0.8,2.5,2.5],segn=15,segR=0.05,tilt_theta=-Pi/7,theta=Pi/2):
canel3:=candy_pole([0.8,4,2.5],segn=15,segR=0.05,tilt_theta=-Pi/7,theta=Pi/2):
left_canes:=display(canel1,canel2,canel3):```

Candy Canes at the right of the house:

```caner1:=candy_pole([3.2,1.5,2.5],segn=15,segR=0.05,tilt_theta=-Pi/7,theta=Pi/2):
caner2:=candy_pole([3.2,2.5,2.5],segn=15,segR=0.05,tilt_theta=-Pi/7,theta=Pi/2):
caner3:=candy_pole([3.2,4,2.5],segn=15,segR=0.05,tilt_theta=-Pi/7,theta=Pi/2):
right_canes:=display(caner1,caner2,caner3):

canes:=display(front_canes,back_canes,right_canes,left_canes):```

With these canes in place all that was left was to create the ground and display our Gingerbread House.

```ground:=display(plottools:-parallelepiped([5,0,0],[0,0.5,0],[0,0,4],[0,1.35,0]),color="DarkGreen"):

display([door,door_decor,d1,base,base_decor,d3,d4,front,back,roof,ground,canes],orientation=[-100,0,95]);```

We are currently in the process of updating the support FAQs at https://faq.maplesoft.com. We’ve been working on updating the existing content for clarity, and have added several new articles already.

The majority of our FAQs are from questions people ask us in Technical Support at support@maplesoft.com, but we’d also like to like to add content from other sources.

Since we have such a great community here at MaplePrimes, we wanted to reach out and ask if there are any articles or questions that you'd like to see added to our FAQ.

We look forward to hearing your feedback!

## Using the Lindstedt-Poincaré Method to A...

Maple

We occasionally get asked questions about methods of Perturbation Theory in Maple, including the Lindstedt-Poincaré Method. Presented here is the most famous application of this method.

Introduction

During the dawn of the 20th Century, one problem that bothered astronomers and astrophysicists was the precession of the perihelion of Mercury. Even when considering the gravity from other planets and objects in the solar system, the equations from Newtonian Mechanics could not account for the discrepancy between the observed and predicted precession.

One of the early successes of Einstein's General Theory of Relativity was that the new model was able to capture the precession of Mercury, in addition to the orbits of all the other planets. The Einsteinian model, when applied to the orbit of Mercury, was in fact a non-negligible perturbation of the old model. In this post, we show how to use Maple to compute the perturbation, and derive the formula for calculating the precession.

In polar coordinates, the Einsteinian model can be written in the following form, where `u(theta)=a(1-e^2)/r(theta)`, with `theta` being the polar angle, `r(theta)` being the radial distance, `a` being the semi-major axis length, and `e` being the eccentricity of the orbit:

``````# Original system.
f := (u,epsilon) -> -1 - epsilon * u^2;
omega := 1;
u0, du0 := 1 + e, 0;
de1 := diff( u(theta), theta, theta ) + omega^2 * u(theta) + f( u(theta), epsilon );
ic1 := u(0) = u0, D(u)(0) = du0;
``````

The small parameter `epsilon` (along with the amount of precession) can be found in terms of the physical constants, but for now we leave it arbitrary:

```# Parameters.
P := [
a = 5.7909050e10 * Unit(m),
c = 2.99792458e8 * Unit(m/s),
e = 0.205630,
G = 6.6740831e-11 * Unit(N*m^2/kg^2),
M = 1.9885e30 * Unit(kg),
alpha = 206264.8062,
beta = 415.2030758
];
epsilon = simplify( eval( 3 * G * M / a / ( 1 - e^2 ) / c^2, P ) ); # approximately 7.987552635e-8```

Note that `c` is the speed of light, `G` is the gravitational constant, `M` is the mass of the sun, `alpha` is the number of arcseconds per radian, and `beta` is the number of revolutions per century for Mercury.

We will show that the radial distance, predicted by Einstein's model, is close to that for an ellipse, as predicted by Newton's model, but the perturbation accounts for the precession (42.98 arcseconds/century). During one revolution, the precession can be determined to be approximately

``sigma = simplify( eval( 6 * Pi * G * M / a / ( 1 - e^2 ) / c^2, P ) ); # approximately 5.018727337e-7``

and so, per century, it is `alpha*beta*sigma`, which is approximately 42.98 arcseconds/century.
It is worth checking out this question on Stack Exchange, which includes an animation generated numerically using Maple for a similar problem that has a more pronounced precession.

Lindstedt-Poincaré Method in Maple

In order to obtain a perturbation solution to the perturbed differential equation `u'+omega^2*u=1+epsilon*u^2` which is periodic, we need to write both `u` and `omega` as a series in the small parameter `epsilon`. This is because otherwise, the solution would have unbounded oscillatory terms ("secular terms"). Using this Lindstedt-Poincaré Method, we substitute arbitrary series in `epsilon` for `u` and `omega` into the initial value problem, and then choose the coefficient constants/functions so that both the initial value problem is satisfied and there are no secular terms. Note that a first-order approximation provides plenty of agreement with the measured precession, but higher-order approximations can be obtained.

To perform this in Maple, we can do the following:

``````# Transformed system, with the new independent variable being the original times a series in epsilon.
de2 := op( PDEtools:-dchange( { theta = phi/b }, { de1 }, { phi }, params = { b, epsilon }, simplify = true ) );
ic2 := ic1;```

```# Order and series for the perturbation solutions of u(phi) and b. Here, n = 1 is sufficient.
n := 1;
U := unapply( add( p[k](phi) * epsilon^k, k = 0 .. n ), phi );
B := omega + add( q[k] * epsilon^k, k = 1 .. n );```

```# DE in terms of the series.
de3 := series( eval( de2, [ u = U, b = B ] ), epsilon = 0, n + 1 );```

```# Successively determine the coefficients p[k](phi) and q[k].
for k from 0 to n do```

```    # Specify the initial conditions for the kth DE, which involves p[k](phi).
# The original initial conditions appear only in the coefficient functions with index k = 0,
# and those for k > 1 are all zero.
if k = 0 then
ic3 := op( expand( eval[recurse]( [ ic2 ], [ u = U, epsilon = 0 ] ) ) );
else
ic3 := p[k](0), D(p[k])(0);
end if:```

```    # Solve kth DE, which can be found from the coefficients of the powers of epsilon in de3, for p[k](phi).
# Then, update de3 with the new information.
soln := dsolve( { simplify( coeff( de3, epsilon, k ) ), ic3 } );
p[k] := unapply( rhs( soln ), phi );
de3 := eval( de3 );```

```    # Choose q[k] to eliminate secular terms. To do this, use the frontend() command to keep only the terms in p[k](phi)
# which have powers of t, and then solve for the value of q[k] which makes the expression zero.
# Note that frontend() masks the t-dependence within the sine and cosine terms.
# Note also that this method may need to be amended, based on the form of the terms in p[k](phi).
if k > 0 then
q[1] := solve( frontend( select, [ has, p[k](phi), phi ] ) = 0, q[1] );
de3 := eval( de3 );
end if;```

`end do:`

```# Final perturbation solution.
'u(theta)' = eval( eval( U(phi), phi = B * theta ) ) + O( epsilon^(n+1) );```

```# Angular precession in one revolution.
sigma := convert( series( 2 * Pi * (1/B-1), epsilon = 0, n + 1 ), polynom ):
epsilon := 3 * G * M / a / ( 1 - e^2 ) / c^2;
'sigma' = sigma;```

```# Precession per century.
xi := simplify( eval( sigma * alpha * beta, P ) ); # returns approximately 42.98``````

Maple Worksheet: Lindstedt-Poincare_Method.mw

## Writing Derivatives at a Point Using Lei...

Maple

Maple users often want to write a derivative evaluated at a point using Leibniz notation, as a matter of presentation, with appropriate variables and coordinates. For instance:

Now, Maple uses the D operator for evaluating derivatives at a point, but this can be a little clunky:

```p := D[1,2,2,3](f)(a,b,c);

q := convert( p, Diff );```

```u := D[1,2,2,3](f)(5,10,15);

v := convert( u, Diff );```

How can we tell Maple, programmatically, to print this in a nicer way? We amended the print command (see below) to do this. For example:

```print( D[1,2,2,3](f)(a,b,c), [x,y,z] );

print( D[1,2,2,3](f)(5,10,15), [x,y,z] );```

`print( 'D(sin)(Pi/6)', theta );`

Here's the definition of the custom version of print:

```# Type to check if an expression is a derivative using 'D', e.g. D(f)(a) and D[1,2](f)(a,b).

'Dexpr',

proc( f )

if op( [0,0], f ) <> D and op( [0,0,0], f ) <> D then

return false;

end if;

if not type( op( [0,1], f ), 'name' ) or not type( { op( f ) }, 'set(algebraic)' ) then

return false;

end if;

if op( [0,0,0], f ) = D and not type( { op( [0,0,..], f ) }, 'set(posint)' ) then

return false;

end if;

return true;

end proc

):

# Create a local version of 'print', which will print expressions like D[1,2](f)(a,b) in a custom way,

# but otherwise print in the usual fashion.

local print := proc()

local A, B, f, g, L, X, Y, Z;

# Check that a valid expression involving 'D' is passed, along with a variable name or list of variable names.

if ( _npassed < 2 ) or ( not _passed[1] :: 'Dexpr' ) or ( not passed[2] :: 'Or'('name','list'('name')) ) then

return :-print( _passed );

end if;

# Extract important variables from the input.

g := _passed[1]; # expression

X := _passed[2]; # variable name(s)

f := op( [0,1], g ); # function name in expression

A := op( g ); # point(s) of evaluation

# Check that the number of variables is the same as the number of evaluation points.

if nops( X ) <> nops( [A] ) then

return :-print( _passed );

end if;

# The differential operator.

L := op( [0,0], g );

# Find the variable (univariate) or indices (multivariate) for the derivative(s).

B := `if`( L = D, X, [ op( L ) ] );

# Variable name(s) as expression sequence.

Y := op( X );

# Check that the point(s) of evaluation is/are distinct from the variable name(s).

if numelems( {Y} intersect {A} ) > 0 then

return :-print( _passed );

end if;

# Find the expression sequence of the variable names.

Z := `if`( L = D, X, X[B] );

return print( Eval( Diff( f(Y), Z ), (Y) = (A) ) );

end proc:

```

Do you use Leibniz Notation often? Or do you have an alternate method? We’d love to hear from you!

 1 2 3 4 5 Page 4 of 5
﻿