Applications, Examples and Libraries

Sampling linearly correlated random variables

by: mmcdara 7831 Product: Maple

August 10 2020

4 0

I found in the Application Center a quite old work (2010) titled Generation of correlated random numbers (see here view.aspx).
This work contains a few errors that I thought it was worth correcting.

Basically the works I refer to concern the sampling of linearly correlated random variables (or correlation in the Pearson sense). Classical textbooks about the subject generally discuss this topic by considering only gaussian random variables and present two methods to generate linearlycorrelated samples: one base on the Cholesky decomposition of the correlation matrix, the other based on its SVD decomposition.

Now the question is: can we apply any of these two procedures to generate linearly correlated samples of arbitrary random variables?
The answer is NO and the reason why it is so is strongly related to a fundamental property of gaussian random variables (GRV) that is that any linear combination of GRVs is still a GRV.
But things are not that simple because even the multi gaussian case handmed with Cholesky's decomposition or SVD can lead to undesired results if no precautions are taken.

The aim of this post is to show what are those wrong results we obtain by thoughtlessly applying these decompositiond and, of course, to show how we must proceed to avoid them.

Let's start by a very simple point of natural good sense: suppose U1 and U2 are two independant identically distributed (iid) random variables and that we have some "function" F which, when applied to the couple (U1, U2) generate a couple (A1, A2) of linearly correlated random variables. Thus F(U1, U2) = (A1, A2).
Let's suppose this same relation holds if we replace U1 and U2 by "a sample of U1" and "a sample of U2", and thus (A1, A2) by "a sample of the bivariate (A1, A2) whose components are linearly correlated". Let's call S the cloud one could obtain by using for instance the ScatterPlot(A1, A2) procedure of Maple.

Let's suppose now that instead of computing F(U1, U2) I decide to compute (U2, U1). Let's call (A1' , A2') the corresponding joint sample and write S' := ScatterPlot(A2', A1').
It seems natural (and it is!) to think that S and S' will be the same up to sampling artifacts.

Any correct method to generate samples from (linearly or not) correlated random variables must verify this similarity of patterns between S and S' S and S'. But this is not the case in this work view.aspx.

The safer way to correlate, even in the Pearson's sense, random variables is to use the concept of COPULAS (there is a work on copulas in the Application Center, but for a quick overview see here Copula_(probability_theory)).
For this special case of linear correlation on can use copulas without knowing it, and this is very simple: as soon as our procedure F introduced above gives correct results if U1 and U2 are standard GRVs,

take any couple (R1, R2) of arbitrary random variables,
build a map M(R1, R2) --> (U1, U2),
generate (A1, A2) = F(U1, U2),
compute M^(-1)(A1, A2)

What is the point of correcting a work that is 10 years old?
A very simple answer is that the Cholesky's decomposition (or SVD) is still the emblematic method to use for linearly correlating random variable. This is the only one presented in scholar textbooks, the only one a lot of students have been taught about (unless they have they have had an extensive background in probability or statistics), and thus a systematic source of wrong results users are not even aware of.

Next point: it's well known that the Pearson's correlation cannot be lower than -1 or higher than +1, but this is common mistake to think any value between -1 and +1 can be reached.
This is guaranteed for GRVs, but not for some other random variables.
For a classical counter-example see 04_correlation_2016_cost_symposium_fkuo_tagged.pdf

The notation used in the attached file are mainly those used in the initial work view.aspx.

>

restart:

This work is aimed to correct the procedure used in https://fr.maplesoft.com/applications/view.aspx?SID=99806
to correlate arbitrary random variables in the (common) Pearson's sense.

>	with(LinearAlgebra): with(plots): with(Statistics):

GAUSSIAN RANDOM VARIABLES

>

# First example: both A1 and A2 are centered gaussian random variables
#                The order we use (A1 next A2 or A2 next A1) to define Ma doesn't matter

Y   := RandomVariable(Normal(0, .25)):
rho := .9:
Q   := 10^4:
A1 := Sample(Y, Q):
A2 := Sample(Y, Q):
Ma := `<,>`(`<,>`(A1), `<,>`(A2)):
MA := Transpose(Ma):
Cor := Matrix([[1, rho], [rho, 1]]):
Cd2 := LUDecomposition(Cor, method = 'Cholesky', output = ['U']):
Rs2 := MA . Cd2:

opts := titlefont = [TIMES, BOLD, 12], symbol=point, transparency=0.5:
A1A2 := ScatterPlot(Column(Rs2, 1), Column(Rs2, 2), title = "Correlated Normal RV", opts, color=blue):

Ma := `<,>`(`<,>`(A2), `<,>`(A1)):
MA := Transpose(Ma):
Cor := Matrix([[1, rho], [rho, 1]]):
Cd2 := LUDecomposition(Cor, method = 'Cholesky', output = ['U']):
Rs2 := MA . Cd2:

A2A1 := ScatterPlot(Column(Rs2, 2), Column(Rs2, 1), title = "Correlated Normal RV", opts, color=red):

display(A1A2, A2A1);

>

# Second example : both A1 and A2 are non-centered gaussian random variables with equal standard deviations.
#                  The order we use to define Ma does matter

Y   := RandomVariable(Normal(1, .25)):
rho := .9:
Q   := 10^4:
A1 := Sample(Y, Q):
A2 := Sample(Y, Q):
Ma := `<,>`(`<,>`(A1), `<,>`(A2)):
MA := Transpose(Ma):
Cor := Matrix([[1, rho], [rho, 1]]):
Cd2 := LUDecomposition(Cor, method = 'Cholesky', output = ['U']):
Rs2 := MA . Cd2:

opts := titlefont = [TIMES, BOLD, 12], symbol=point, transparency=0.5:

A1A2 := ScatterPlot(Column(Rs2, 1), Column(Rs2, 2), title = "Correlated Normal RV", opts, color=blue):

Ma := `<,>`(`<,>`(A2), `<,>`(A1)):
MA := Transpose(Ma):
Cor := Matrix([[1, rho], [rho, 1]]):
Cd2 := LUDecomposition(Cor, method = 'Cholesky', output = ['U']):
Rs2 := MA . Cd2:

A2A1 := ScatterPlot(Column(Rs2, 2), Column(Rs2, 1), title = "Correlated Normal RV", opts, color=red):

display(A1A2, A2A1);

>

# Second example corrected: to avoid order's dependency proceed this way
#    1/ center A1 and A2
#    2/ correlate the now centered rvs
#    3/ uncenter the couple of correlated rvs

C1 := convert(Scale(A1, scale=Mean), Vector[row]):
C2 := convert(Scale(A2, scale=Mean), Vector[row]):

Ma := `<,>`(`<,>`(C1), `<,>`(C2)):
MA := Transpose(Ma):
Cor := Matrix([[1, rho], [rho, 1]]):
Cd2 := LUDecomposition(Cor, method = 'Cholesky', output = ['U']):
Rs2 := MA . Cd2:

opts := titlefont = [TIMES, BOLD, 12], symbol=point, transparency=0.5:

A1A2 := ScatterPlot(Column(Rs2, 1)+~Mean(A1), Column(Rs2, 2)+~Mean(A2), title = "Correlated Normal RV", opts, color=blue):

Ma := `<,>`(`<,>`(C2), `<,>`(C1)):
MA := Transpose(Ma):
Cor := Matrix([[1, rho], [rho, 1]]):
Cd2 := LUDecomposition(Cor, method = 'Cholesky', output = ['U']):
Rs2 := MA . Cd2:

A2A1 := ScatterPlot(Column(Rs2, 2)+~Mean(A1), Column(Rs2, 1)+~Mean(A2), title = "Correlated Normal RV", opts, color=red):

display(A1A2, A2A1);

>

# Third example : both A1 and A2 are centered gaussian random variables with unequal standard deviations.
# The order we use to define Ma does matter

rho := .9:
Q := 10^4:
A1 := Sample(Normal(0, 1), Q):
A2 := Sample(Normal(0, 2), Q):
Ma := `<,>`(`<,>`(A1), `<,>`(A2)):
MA := Transpose(Ma):
Cor := Matrix([[1, rho], [rho, 1]]):
Cd2 := LUDecomposition(Cor, method = 'Cholesky', output = ['U']):
Rs2 := MA . Cd2:

opts := titlefont = [TIMES, BOLD, 12], symbol=point, transparency=0.5:

A1A2 := ScatterPlot(Column(Rs2, 1), Column(Rs2, 2), title = "Correlated Normal RV", opts, color=blue):

Ma := `<,>`(`<,>`(A2), `<,>`(A1)):
MA := Transpose(Ma):
Cor := Matrix([[1, rho], [rho, 1]]):
Cd2 := LUDecomposition(Cor, method = 'Cholesky', output = ['U']):
Rs2 := MA . Cd2:

A2A1 := ScatterPlot(Column(Rs2, 2), Column(Rs2, 1), title = "Correlated Normal RV", opts, color=red):

display(A1A2, A2A1);

>

# Third example corrected: to avoid order's dependency proceed this way
#    1/ scale A1 and A2
#    2/ correlate the now scaled rvs
#    3/ unscale the couple of correlated rvs

C1 := A1 /~ StandardDeviation(A1):
C2 := A2 /~ StandardDeviation(A2):

Ma := `<,>`(`<,>`(C1), `<,>`(C2)):
MA := Transpose(Ma):
Cor := Matrix([[1, rho], [rho, 1]]):
Cd2 := LUDecomposition(Cor, method = 'Cholesky', output = ['U']):
Rs2 := MA . Cd2:

opts := titlefont = [TIMES, BOLD, 12], symbol=point, transparency=0.5:

A1A2 := ScatterPlot(Column(Rs2, 1)*~StandardDeviation(A1), Column(Rs2, 2)*~StandardDeviation(A2), title = "Correlated Normal RV", opts, color=blue):

Ma := `<,>`(`<,>`(C2), `<,>`(C1)):
MA := Transpose(Ma):
Cor := Matrix([[1, rho], [rho, 1]]):
Cd2 := LUDecomposition(Cor, method = 'Cholesky', output = ['U']):
Rs2 := MA . Cd2:

A2A1 := ScatterPlot(Column(Rs2, 2)*~StandardDeviation(A1), Column(Rs2, 1)*~StandardDeviation(A2), title = "Correlated Normal RV", opts, color=red):

display(A1A2, A2A1);

>	# More generally: to avoid order's dependency proceed this way # 1/ transform A1 and A2 into standard gaussian random variables (mean and standard deviation scalings) # 2/ correlate the now scaled rvs # 3/ unscale the couple of correlated rvs

A MORE COMPLEX EXAMPLE:

NON GAUSSIAN RANDOM VARIABLES
(here two LogNormal rvs)

>

# Preliminary
#   the expectation (mean) of a LogNormal rv cannot be 0;
#   as a consequence it is expected that the order used to buid Ma will matter
#
# Proceed as Igor Hlivka did

Y   := RandomVariable(LogNormal(.5, .25)):
rho := .9:
Q   := 1000:
A1 := Sample(Y, Q):
A2 := Sample(Y, Q):
Ma := `<,>`(`<,>`(A1), `<,>`(A2)):
MA := Transpose(Ma):
Cor := Matrix([[1, rho], [rho, 1]]):
Cd2 := LUDecomposition(Cor, method = 'Cholesky', output = ['U']):
Rs2 := MA . Cd2:

ScatterPlot(A1, A2, color = red, title = ["Raw LogNormal RV", font = [TIMES, BOLD, 12]]):
A1A2 := ScatterPlot(Column(Rs2, 1), Column(Rs2, 2), title = "Correlated LogNormal RV", opts, color=blue):

# And now change, as usual, the order in Ma

Ma := `<,>`(`<,>`(A2), `<,>`(A1)):
MA := Transpose(Ma):
Cor := Matrix([[1, rho], [rho, 1]]):
Cd2 := LUDecomposition(Cor, method = 'Cholesky', output = ['U']):
Rs2 := MA . Cd2:

A2A1 := ScatterPlot(Column(Rs2, 2), Column(Rs2, 1), title = "Correlated LogNormal RV", opts, color=red):

display(A1A2, A2A1);

>

# How can we avoid that the order used to assemble Ma do matter?
#
# A close examination of what was done with gaussiann rvs show that in all the cases we
# went back to standard gaussian rvs before correlating them.
# So let's just do the same thing here.
#
# Of course it's not as immediate as previously...
# (please do not focus on the slowness of the code, it is written to clearly explain
# what is done, not to be fast)

#-------------------------------------- from Y to standard gaussian
G := RandomVariable(Normal(0, 1)):
G1 := Vector[row](Q, q -> Quantile(G, Probability(Y > A1[q], numeric), numeric)):
G2 := Vector[row](Q, q -> Quantile(G, Probability(Y > A2[q], numeric), numeric)):
# could be replaced by this faster code
#   cdf_Y := unapply(CDF(Y, z), z) assuming z > 0;
#   cdf_G := unapply(CDF(G, z), z);
#   S1    := sort(A1):
#   ini   := -10:
#   V     := Vector[row](Q):
#   for q from 1 to Q do
#     V[q] := fsolve(cdf_G(z)=cdf_Y(S1[q]), z=ini);
#     ini := V[q]:
#   end do:
#------------------------------------------------------------------

Ma := `<,>`(`<,>`(G1), `<,>`(G2)):
MA := Transpose(Ma):
Cor := Matrix([[1, rho], [rho, 1]]):
Cd2 := LUDecomposition(Cor, method = 'Cholesky', output = ['U']):
Rs2 := MA . Cd2:

opts := titlefont = [TIMES, BOLD, 12], symbol=point, transparency=0.5:

#-------------------------------------- from standard gaussian to Y
C1 := Vector[row](Q, q -> Quantile(Y, Probability(G > Rs2[q, 1], numeric), numeric)):
C2 := Vector[row](Q, q -> Quantile(Y, Probability(G > Rs2[q, 2], numeric), numeric)):
#------------------------------------------------------------------
A1A2 := ScatterPlot(C1, C2, title = "Correlated Normal RV", opts, color=blue):

Ma := `<,>`(`<,>`(G2), `<,>`(G1)):
MA := Transpose(Ma):
Cor := Matrix([[1, rho], [rho, 1]]):
Cd2 := LUDecomposition(Cor, method = 'Cholesky', output = ['U']):
Rs2 := MA . Cd2:

#-------------------------------------- from standard gaussian to Y
C1 := Vector[row](Q, q -> Quantile(Y, Probability(G > Rs2[q, 1], numeric), numeric)):
C2 := Vector[row](Q, q -> Quantile(Y, Probability(G > Rs2[q, 2], numeric), numeric)):
#------------------------------------------------------------------

A2A1 := ScatterPlot(C2, C1, title = "Correlated Normal RV", opts, color=red):

display(A1A2, A2A1);

CONCLUSION: Be extremely careful when correlating non standard gaussian random variables,
and more generally non gaussian random variables.

Correlating rvs the way Igor Hlivka did can be replaced in the more general framework of COPULA THEORY.

Mathematically a bidimensional copula C is a function from [0, 1] x [0, 1] to [0, 1] if C is joint CDF of a bivariate random variable
both with uniform marginals on [0, 1].
See for instanc here https://en.wikipedia.org/wiki/Copula_(probability_theory)

What I did here to "correlate" A1 and A2 was nothing but to apply in a step-by-step way a GAUSSIAN COPULA to the bivariate
(A1, A2) random variable.
In Quantile(G, Probability(Y > A1[q], numeric), numeric) the blue expression maps A1 onto [0, 1] (as it is needed
in the definition of a copula), while the brown sequence is the copula itself (when the same operation on A2 has been done).

>

Download LInear_Correlated_Random_Variables.mw

International Mathematics Competition for University...

by: vv 13750 Product: Maple

August 01 2020

5 4

This year, the International Mathematics Competition for University Students (IMC) took place online (due to Coronavirus), https://www.imc-math.org.uk/?year=2020

One of the sponsors was Maplesoft.

Here is a Maple solution for one of the most difficult problems.

Problem 4, Day 1.

A polynomial with real coeffcients satisfies the equation

, for all real .

Prove that for .

A Maple solution.

Obviously, the degree of the polynomial must be 101.

We shall find effectively p(x).

>

restart;

>

n:=100;

(1)

>	p:= x -> add(a[k]*x^k, k=0..n+1):

>	collect(expand( p(x+1) - p(x) - x^n ), x):

>	S:=solve([coeffs(%,x)]):

>	f:=unapply(expand(eval(p(1-x)-p(x), S)), x);

(2)

>	plot(f, 0..1); # Visual check: f(x)>0 for 0<x<1/2

>	f(0), f(1/4), f(1/2);

0, 2903528346661097497054603834764435875077553006646158945080492319146997643370625023889353447129967354174648294748510553528692457632980625125/3213876088517980551083924184682325205044405987565585670602752, 0

(3)

>	sturm(f(x), x, 0, 1/2);

(4)

So, the polynomial f has a unique zero in the interval (0, 1/2]. Since f(1/2) = 0 and f(1/4) > 0, it results that f > 0 in the interval (0, 1/2). Q.E.D.

Download imc2020-1-4.mw

Circumscribed ellipse

by: one man 1355 Product: Maple 17

July 31 2020

0 2

One way to find the equation of an ellipse circumscribed around a triangle. In this case, we solve a linear system of equations, which is obtained after fixing the values of two variables ( t1 and t2). These are five equations: three equations of the second-order curve at three vertices of the triangle and two equations of a linear combination of the coordinates of the gradient of the curve equation.
The solving of system takes place in the ELS procedure. When solving, hyperboles appear, so the program has a filter. The filter passes the equations of ellipses based on by checking the values of the invariants of the second-order curves.
FOR_ELL_ТR_OUT_PROCE_F.mw ( Fixed comments in the text 01, 08, 2020)

Doing plots "a la XKCD"

by: mmcdara 7831 Product: Maple

July 25 2020

6 0

A lot of scientific software propose packages enabling drawing figures in XKCD style/
Up to now I thought this was restricted to open products (R, Python, ...) but I recently discovered Matlab and even Mathematica were doing same.

Layton S (2012). “XKCDIFY! Adding flair to boring Matlab Axes one plot at a time.” Last accessed on December 08, 2014, URL https://github.com/slayton/matlab-xkcdify.

Woods S (2012). “xkcd-style graphs.” Last accessed on December 08, 2014, URL http://mathematica.stackexchange.com/questions/11350/xkcd-style-graphs/ 11355#11355.

So why not Maple?

As a regular user of R, I could have visualize the body of the corresponding procedures to see how these drawings were made and just translate theminto Maple.
But copying for the sake of copying is not of much interest.
So I started to develop some primitives for "XKCD-drawing" lines, polygons, circles and even histograms.
My goal is not to write an XKCD package (I don't have the skills for that) but just to arouse the interest of (maybe) a few people here who could continue this preliminary work

A main problem is the one of the XKCD fonts: no question to redefine them in Maple and I guess using them in a commercial code is not legal (?). So no XKCD font in this first work, nor even the funny guy who appears recurently on the drawings (but it could be easily constructed in Maple).

In a recent post (Plot styling - experimenting with Maple's plotting...) Samir Khan proposed a few styles made of several plotting options, some of which he named "Excel style" or "Oscilloscope style"... maybe a future "XKCD xtyle" in Maple ?

This work has been done with Maple 2015 and reuses an old version of a 1D-Kriging procedure

>

restart:

>	with(LinearAlgebra): with(plots): with(Statistics):

The principle is always the same:
    1/   Let L a straight line which is either defined by its two ending points (xkvd_hline) or taken as the default [0, 0], [1, 0] line.
          For xkvd_hline the given line L is firstly rotate to be aligned with the horizontal axis.

    2/   Let P1, ..., PN N points on L. Each Pn writes [xn, yn]

    3/   A random perturbation rn is added yo the values y1, ..., yN

    4/   A stationnary random process RP, with gaussian correlation function is used to build a smooth curve passing through the points
          (x1, y1+r1), ..., (xN, yN+rN) (procedure KG where "KG" stands for "Kriging")

    5/   The result is drawn or mapped to some predefined shape :
                  xkcd_hist,
                  xkcd_polyline,
                  xkcd_circle

    6/   A procedure xkcd_func is also provided to draw functions defined by an explicit relation.

>

KG := proc(X, Y, psi, sigma)
  local NX, DX, K, mu, k, y:
  NX := numelems(X);
  DX := < seq(Vector[row](NX, X), k=1..NX) >^+:
  K := (sigma, psi) -> evalf( sigma^2 *~ exp~(-((DX - DX^+) /~ psi)^~2) ):
  mu := add(Y) / NX;
  k := (x, sigma, psi) -> evalf( convert(sigma^2 *~ exp~(-((x -~ X ) /~ psi)^~2), Vector[row]) ):
  y := mu + k(x, sigma, psi) . (K(sigma, psi))^(-1) . convert((Y -~ mu), Vector[column]):
  return y
end proc:

xkcd_hline := proc(p1::list, p2::list, a::nonnegative, lc::positive, col)
  # p1 : first ending point
  # p2 : second ending point
  # a : amplitude of the random perturbations
  # lc : correlation length
  # col: color
  local roll, NX, LX, X, Z:
  roll := rand(-1.0 .. 1.0):
  NX   := 10:
  LX   := p2[1]-p1[1]:
  X    := [seq(p1[1]..p2[1], LX/(NX-1))]:
  Z    := [p1[2], seq(p1[2]+a*roll(), k=1..NX-1)]:
  return plot(KG(X, Z, lc*LX, 1), x=min(X)..max(X), color=col, scaling=constrained):
end proc:

xkcd_line := proc(L::list, a::nonnegative, lc::positive, col, {lsty::integer:=1})
  # L : list which contains the two ending point
  # a : amplitude of the random perturbations
  # lc : correlation length
  # col: color
  local T, roll, NX, DX, DY, LX, A, m, M, X, Z, P:
  T    := (a, x0, y0, l) ->
             plottools:-transform(
               (x,y) -> [ x0 + l * (x*cos(a)-y*sin(a)), y0 + l * (x*sin(a)+y*cos(a)) ]
             ):
  roll := rand(-1.0 .. 1.0):
  NX   := 5:
  DX   := L[2][1]-L[1][1]:
  DY   := L[2][2]-L[1][2]:
  LX := sqrt(DX^2+DY^2):
  if DX <> 0 then
     A := arcsin(DY/LX):
  else
     A:= Pi/2:
  end if:
  X := [seq(0..1, 1/(NX-1))]:
  Z := [ seq(a*roll(), k=1..NX)]:
  P := plot(KG(X, Z, lc, 1), x=0..1, color=col, scaling=constrained, linestyle=lsty):
  return T(A, op(L[1]), LX)(P)
end proc:

xkcd_func := proc(f, r::list, NX::posint, a::positive, lc::positive, col)
  # f : function to draw
  # r : plot range
  # NX : number of equidistant "nodes" in the range r (boundaries included)
  # a : amplitude of the random perturbations
  # lc : correlation length
  # col: color
  local roll, F, LX, Pf, Xf, Zf:
  roll := rand(-1.0 .. 1.0):
  F    := unapply(f, indets(f, name)[1]);
  LX   := r[2]-r[1]:
  Pf   := [seq(r[1]..r[2], LX/(NX-1))]:
  Xf   := Pf +~ [seq(a*roll(), k=1..numelems(Pf))]:
  Zf   := F~(Pf) +~ [seq(a*roll(), k=1..numelems(Pf))]:
  return plot(KG(Xf, Zf, lc*LX, 1), x=min(Xf)..max(Xf), color=col):
end proc:

xkcd_hist := proc(H, ah, av, ax, ay, lch, lcv, lcx, lcy, colh, colxy)
  # H   : Histogram
  # ah : amplitude of the random perturbations on the horizontal boundaries of the bins
  # av : amplitude of the random perturbations on the vertical boundaries of the bins
  # ax : amplitude of the random perturbations on the horizontal axis
  # ay : amplitude of the random perturbations on the vertical axis
  # lch : correlation length on the horizontal boundaries of the bins
  # lcv : correlation length on the vertical boundaries of the bins
  # lcx : correlation length on the horizontal axis
  # lcy : correlation length on the vertical axis
  # colh: color of the histogram
  # col : color of the axes
  local data, horiz, verti, horizontal_lines, vertical_lines, po, rpo, p1, p2:
  data := op(1..-2, op(1, H)):
  verti := sort( [seq(data[n][3..4][], n=1..numelems([data]))] , key=(x->x[1]) );
  verti := verti[1],
           map(
                n -> if verti[n][2] > verti[n+1][2] then
                        verti[n]
                      else
                        verti[n+1]
                      end if,
                [seq(2..numelems(verti)-2,2)]
           )[],
           verti[-1];
  horiz := seq(data[n][[4, 3]], n=1..numelems([data])):

  horizontal_lines := NULL:
  for po in horiz do
    horizontal_lines := horizontal_lines, xkcd_hline(po[1], po[2], ah, lch, colh):
  end do:

  vertical_lines := NULL:
  for po in [verti] do
    rpo := po[[2, 1]]:
    vertical_lines := vertical_lines, xkcd_hline([0, rpo[2]], rpo, av, lcv, colh):
  end do:

  p1 := [2*verti[1][1]-verti[2][1], 0]:
  p2 := [2*verti[-1][1]-verti[-2][1], 0]:

  return
    display(
      horizontal_lines,
      T~([vertical_lines]),
      xkcd_hline(p1, p2, ax, lcx, colxy),
      T(xkcd_hline([0, 0], [1.2*max(op~(2, [verti])), 0], ay, lcy, colxy)),
      axes=none,
      scaling=unconstrained
    );
end proc:

xkcd_polyline := proc(L::list, a::nonnegative, lc::positive, col)
  # xkcd_polyline reduces to xkcd_line whebn L has 2 elements
  # L : List of points
  # a : amplitude of the random perturbations
  # lc : correlation length
  # col: color
  local T, roll, NX, n, DX, DY, LX, A, m, M, X, Z, P:
  T    := (a, x0, y0, l) ->
             plottools:-transform(
               (x,y) -> [ x0 + l * (x*cos(a)-y*sin(a)), y0 + l * (x*sin(a)+y*cos(a)) ]
             ):
  roll := rand(-1.0 .. 1.0):
  NX   := 5:
  for n from 1 to numelems(L)-1 do
    DX   := L[n+1][1]-L[n][1]:
    DY   := L[n+1][2]-L[n][2]:
    LX := sqrt(DX^2+DY^2):
    if DX <> 0 then
      A := evalf(arcsin(abs(DY)/LX)):
      if DX >= 0 and DY <= 0 then A := -A end if:
      if DX <= 0 and DY > 0 then A := Pi-A end if:
      if DX <= 0 and DY <= 0 then A := Pi+A end if:
    else
      A:= Pi/2:
      if DY < 0 then A := 3*Pi/2 end if:
    end if:
    if n=1 then
      X := [seq(0..1, 1/(NX-1))]:
      Z := [seq(a*roll(), k=1..NX)]:
    else
      X := [0    , seq(1/(NX-1)..1, 1/(NX-1))]:
      Z := [Z[NX], seq(a*roll(), k=1..NX-1)]:
    end if:
    P    := plot(KG(X, Z, lc, 1), x=0..1, color=col, scaling=constrained):
    P||n := T(A, op(L[n]), LX)(P):
  end do;
  return seq(P||n, n=1..numelems(L)-1)
end proc:

xkcd_circle := proc(a::nonnegative, lc::positive, r::positive, cent::list, col)
  # a   : amplitude of the random perturbations
  # lc : correlation length
  # r   : redius of the circle
  # cent: center of the circle
  # col : color
  local roll, NX, LX, X, Z, xkg, A:
  roll := rand(-1.0 .. 1.0):
  NX   := 10:
  X    := [seq(0..1, 1/(NX-1))]:
  Z    := [0, seq(a*roll(), k=1..NX-1)]:
  xkg := KG(X, Z, lc, 1):
  A    := Pi*roll():
  return plot([cent[1]+r*(1+xkg)*cos(2*Pi*x+A), cent[2]+r*(1+xkg)*sin(2*Pi*x+A), x=0..1], color=col)
end proc:

T := plottools:-transform((x,y) -> [y, x]):

>	# Axes plot x_axis := xkcd_hline([0, 0], [10, 0], 0.03, 0.5, black): y_axis := xkcd_hline([0, 0], [10, 0], 0.03, 0.5, black): display( x_axis, T(y_axis), axes=none, scaling=constrained )

>	# A simple function f := 1+10*(x/5-1)^2: F := xkcd_func(f, [0.5, 9.5], 6, 0.3, 0.4, red): display( x_axis, T(y_axis), F, axes=none, scaling=constrained )

>	# An histogram S := Sample(Normal(0,1),100): H := Histogram(S, maxbins=6): xkcd_hist(H, 0, 0.02, 0.001, 0.01, 1, 0.1, 0.01, 1, red, black)

>

# Axes plus grid with two red straight lines

r := rand(-0.1 .. 0.1):

x_axis := xkcd_line([[-2, 0], [10, 0]], 0.01, 0.2, black):
y_axis := xkcd_line([[0, -2], [0, 10]], 0.01, 0.2, black):
d1     := xkcd_line([[-1, 1], [9, 9]] , 0.01, 0.2, red):
d2     := xkcd_line([[-1, 9], [9, -1]], 0.01, 0.2, red):
display(
  x_axis, y_axis,
  seq( xkcd_line([[-2+0.3*r(), u+0.3*r()], [10+0.3*r(), u+0.3*r()]], 0.005, 0.5, gray), u in [seq(-1..9, 2)]),
  seq( xkcd_line([[u+0.3*r(), -2+0.3*r()], [u+0.3*r(), 10+0.3*r()]], 0.005, 0.5, gray), u in [seq(-1..9, 2)]),
  d1, d2,
  axes=none,
  scaling=constrained
)

>	# Axes and a couple of polygonal lines d1 := xkcd_polyline([[0, 0], [1, 3], [3, 5], [7, 1], [9, 7]], 0.01, 1, red): d2 := xkcd_polyline([[0, 9], [2, 8], [5, 2], [8, 3], [10, -1]], 0.01, 1, blue): display( x_axis, y_axis, d1, d2, axes=none, scaling=constrained )

>

# A few polygonal shapes

display(
  xkcd_polyline([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]], 0.01, 1, red),
  xkcd_polyline([[1/3, 1/3], [2/3, 1/3], [2/3, 4/3], [-1, 4/3], [1/3, 1/3]], 0.01, 1, blue),
  xkcd_polyline([[-1/3, -1/3], [4/3, 1/2], [1/2, 1/2], [1/2,-1], [-1/3, -1/3]], 0.01, 1, green),
  axes=none,
  scaling=constrained
)

>

# A few circles

cols := [red, green, blue, gold, black]:                                # colors
cents := convert( Statistics:-Sample(Uniform(-1, 3), [5, 2]), listlist): # centers
radii := Statistics:-Sample(Uniform(1/2, 2), 5):                         # radii
lcs   := Statistics:-Sample(Uniform(0.2, 0.7), 5):                       # correlations lengths

display(
  seq(
    xkcd_circle(0.02, lcs[n], radii[n], cents[n], cols[n]),
    n=1..5
  ),
  axes=none,
  scaling=constrained
)

>

# A 3D drawing

x_axis := xkcd_line([[0, 0], [5, 0]], 0.01, 0.2, black):
y_axis := xkcd_line([[0, 0], [4, 2]], 0.01, 0.2, black):
z_axis := xkcd_line([[0, 0], [0, 5]], 0.01, 0.2, black):

f1 := 4*cos(x/6)-1:
F1 := xkcd_func(f1, [0.5, 5], 6, 0.001, 0.8, red):
F2 := xkcd_line([[0.5, eval(f1, x=0.5)], [3, 4]], 0.01, 0.2, red):
f3 := 4*cos((x-2)/6):
F3 := xkcd_func(f3, [3, 7], 6, 0.001, 0.8, red):
F4 := xkcd_line([[5, eval(f1, x=5)], [7, eval(f3, x=7)]], 0.01, 0.2, red):

dx := xkcd_line([[2, 1], [4, 1]], 0.01, 0.2, gray, lsty=3):
dy := xkcd_line([[2, 0], [4, 1]], 0.01, 0.2, gray, lsty=3):
dz := xkcd_line([[4, 1], [4, 3]], 0.01, 0.2, gray, lsty=3):

po := xkcd_circle(0.02, 0.3, 0.1, [4, 3], blue):

# Numerical value come from "probe info + copy/paste"

nvect   := xkcd_polyline([[4, 3], [4.57, 4.26], [4.35, 4.1], [4.57, 4.26], [4.58, 4.02]], 0.01, 1, blue):
tg1vect := xkcd_polyline([[4, 3], [4.78, 2.59], [4.49, 2.87], [4.78, 2.59], [4.46, 2.57]], 0.01, 1, blue):
tg2vect := xkcd_polyline([[4, 3], [4.79, 3.35], [4.70, 3.13], [4.79, 3.35], [4.46, 3.35]], 0.01, 1, blue):
rec1    := xkcd_polyline([[4.118, 3.286], [4.365, 3.396], [4.257, 3.108]], 0.01, 1, blue):
rec2    := xkcd_polyline([[4.257, 3.108], [4.476, 2.985], [4.259, 2.876]], 0.01, 1, blue):

display(
  x_axis, y_axis, z_axis,
  F1, F2, F3, F4,
  dx, dy, dz,
  po,
  nvect, tg1vect, tg2vect, rec1, rec2,
  axes=none,
  scaling=constrained
)

>

# Arrow

d1 := xkcd_polyline([[0, 0], [1, 0], [0.9, 0.05], [1, 0], [0.9, -0.05]], 0.01, 1, red):

T := (a, x0, y0, l) ->
             plottools:-transform(
               (x,y) -> [ x0 + l * (x*cos(a)-y*sin(a)), y0 + l * (x*sin(a)+y*cos(a)) ]
             ):

display(
  seq( T(2*Pi*n/10, 0.5, 0, 1/2)(
           display(
              xkcd_polyline(
                  [[0, 0], [1, 0], [0.9, 0.05], [1, 0], [0.9, -0.05]],
                  0.01,
                  1,
                  ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12])
               )
           )
        ),
       n=1..10
  ),
  axes=none,
  scaling=constrained
)

>

Download XKCD.mw

An ellipse inscribed in a given triangle

by: one man 1355 Product: Maple 17

July 24 2020

3 11

An attempt to find the equation of an ellipse inscribed in a given triangle.
The program works on the basis of the ELS procedure. After the procedure works, the solutions are filtered.
ELS procedure solves the system of equations f1, f2, f3, f4, f5 for the coefficients of the second-order curve.
The equation f1 corresponds to the condition that the side of the triangle intersects t a curve of the second order at one point.
The equation f2 corresponds to the condition that the point x1,x2 belongs to a curve of the second order.
Equation f3 corresponds to the condition that the side of the triangle is tangent to the second order curve at the point x1,x2.
The equation f4 is similar to the equation f2, and the equation f5 is similar to the equation f3.
FOR_ELL_ТR_PROCE.mw
For example

Plot styling - experimenting with Maple's plotting...

by: Samir Khan 2091 Product: Maple

July 13 2020

11 4

I like tweaking plots to get the look and feel I want, and luckily Maple has many plotting options that I often play with. Here, I visualize the same data several times, but each time with different styling.

First, some data.

restart:
data_1 := [[0,0],[1,2],[2,1.3],[3,6]]:
data_2 := [[0.5,3],[1,1],[2,5],[3,2]]:
data_3 := [[-0.5,3],[1.3,1],[2.5,5],[4.5,2]]:

This is the default look.

plot([data_1, data_2, data_3])

I think the darker background on this plot makes it easier to look at.

gray_grid :=
 background      = "LightGrey"
,color           = [ ColorTools:-Color("RGB",[150/255, 40 /255, 27 /255])
                    ,ColorTools:-Color("RGB",[0  /255, 0  /255, 0  /255])
                    ,ColorTools:-Color("RGB",[68 /255, 108/255, 179/255]) ]
,axes            = frame
,axis[2]         = [color = black, gridlines = [10, thickness = 1, color = ColorTools:-Color("RGB", [1, 1, 1])]]
,axis[1]         = [color = black, gridlines = [10, thickness = 1, color = ColorTools:-Color("RGB", [1, 1, 1])]]
,axesfont        = [Arial]
,labelfont       = [Arial]
,size            = [400*1.78, 400]
,labeldirections = [horizontal, vertical]
,filled          = false
,transparency    = 0
,thickness       = 5
,style           = line:

plot([data_1, data_2, data_3], gray_grid);

I call the next style Excel, for obvious reasons.

excel :=
 background      = white
,color           = [ ColorTools:-Color("RGB",[79/255,  129/255, 189/255])
                    ,ColorTools:-Color("RGB",[192/255, 80/255,   77/255])
                    ,ColorTools:-Color("RGB",[155/255, 187/255,  89/255])]
,axes            = frame
,axis[2]         = [gridlines = [10, thickness = 0, color = ColorTools:-Color("RGB",[134/255,134/255,134/255])]]
,font            = [Calibri]
,labelfont       = [Calibri]
,size            = [400*1.78, 400]
,labeldirections = [horizontal, vertical]
,transparency    = 0
,thickness       = 3
,style           = point
,symbol          = [soliddiamond, solidbox, solidcircle]
,symbolsize      = 15:

plot([data_1, data_2, data_3], excel)

This style makes the plot look a bit like the oscilloscope I have in my garage.

dark_gridlines :=
 background      = ColorTools:-Color("RGB",[0,0,0])
,color           = white
,axes            = frame
,linestyle       = [solid, dash, dashdot]
,axis            = [gridlines = [10, linestyle = dot, color = ColorTools:-Color("RGB",[0.5, 0.5, 0.5])]]
,font            = [Arial]
,size            = [400*1.78, 400]:

plot([data_1, data_2, data_3], dark_gridlines);

The colors in the next style remind me of an Autumn morning.

autumnal :=
 background      =  ColorTools:-Color("RGB",[236/255, 240/255, 241/255])
,color           = [  ColorTools:-Color("RGB",[144/255, 54/255, 24/255])
                     ,ColorTools:-Color("RGB",[105/255, 108/255, 51/255])
                     ,ColorTools:-Color("RGB",[131/255, 112/255, 82/255]) ]
,axes            = frame
,font            = [Arial]
,size            = [400*1.78, 400]
,filled          = true
,axis[2]         = [gridlines = [10, thickness = 1, color = white]]
,axis[1]         = [gridlines = [10, thickness = 1, color = white]]
,symbol          = solidcircle
,style           = point
,transparency    = [0.6, 0.4, 0.2]:

plot([data_1, data_2, data_3], autumnal);

In honor of a friend and ex-colleague, I call this style "The Swedish".

swedish :=
 background      = ColorTools:-Color("RGB", [0/255, 107/255, 168/255])
,color           = [ ColorTools:-Color("RGB",[169/255, 158/255, 112/255])
                    ,ColorTools:-Color("RGB",[126/255,  24/255,   9/255])
                    ,ColorTools:-Color("RGB",[254/255, 205/255,   0/255])]
,axes            = frame
,axis            = [gridlines = [10, color = ColorTools:-Color("RGB",[134/255,134/255,134/255])]]
,font            = [Arial]
,size            = [400*1.78, 400]
,labeldirections = [horizontal, vertical]
,filled          = false
,thickness       = 10:

plot([data_1, data_2, data_3], swedish);

This looks like a plot from a journal article.

experimental_data_mono :=

background       = white
,color           = black
,axes            = box
,axis            = [gridlines = [linestyle = dot, color = ColorTools:-Color("RGB",[0.5, 0.5, 0.5])]]
,font            = [Arial, 11]
,legendstyle     = [font = [Arial, 11]]
,size            = [400, 400]
,labeldirections = [horizontal, vertical]
,style           = point
,symbol          = [solidcircle, solidbox, soliddiamond]
,symbolsize      = [15,15,20]:

plot([data_1, data_2, data_3], experimental_data_mono, legend = ["Annihilation", "Authority", "Acceptance"]);

If you're willing to tinker a little bit, you can add some real character and personality to your visualizations. Try it!

I'd also be very interested to learn what you find attractive in a plot - please do let me know.

Realistic 3D rendering of a complex system

by: mmcdara 7831 Product: Maple 2015

July 02 2020

5 1

Hi,

I would like to share this work I've done.
No big math here, just a demonstrator of Maple's capabilities in 3D visualization.

All the plots in the file have been discarded to reduce the size of this post. Here is a screen capture to give you an idea of what is inside the file.

Download 3D_Visualization.mw

Numerical evaluation of integrals over polygons...

by: mmcdara 7831

June 28 2020

2 2

Hi,

In a recent post (Monte Carlo Integration) Radaar shared its work about the numerical integration, with the Monte Carlo method, of a function defined in polar coordinates.
Radaar used a raw strategy based on a sampling in cartesian coordinates plus an ad hoc transformation.
Radaar obtained reasonably good results, but I posted a comment to show how Monte Carlo summation in polar coordinates can be done in a much simpler way. Behind this is the choice of a "good" sampling distribution which makes the integration problem as simple as Monte Carlo integration over a 2D rectangle with sides parallel to the co-ordinate axis.

This comment I sent pushed me to share the present work on Monte Carlo integration over simple polygons ("simple" means that two sides do not intersect).
Here again one can use raw Monte Carlo integration on the rectangle this polygon is inscribed in. But as in Radaar's post, a specific sampling distribution can be used that makes the summation method more elegant.

This work relies on three main ingredients:

The Dirichlet distribution, whose one form enables sampling the 2D simplex in a uniform way.
The construction of a 1-to-1 mapping from this simplex into any non degenerated triangle (a mapping whose jacobian is a constant equal to the ratio of the areas of the two triangles).
A tesselation into triangles of the polygon to integrate over.

This work has been carried out in Maple 2015, which required the development of a module to do the tesselation. Maybe more recent Maple's versions contain internal procedures to do that.

Monte_Carlo_Integration.mw

Why Staying at Home is Good to Avoid the Spread...

by: Roanes-Lozano 70 Product: Maple

June 16 2020

13 0

Hi. My name is Eugenio and I’m a Professor at the Departamento de Didáctica de las Ciencias Experimentales, Sociales y Matemáticas at the Facultad de Educación of the Universidad Complutense de Madrid (UCM) and a member of the Instituto de Matemática Interdisciplinar (IMI) of the UCM.

I have a 14-year-old son. In the beginning of the pandemic, a confinement was ordered in Spain. It is not easy to make a kid understand that we shouldn't meet our friends and relatives for some time and that we should all stay at home in the first stage. So, I developed a simplified explanation of virus propagation for kids, firstly in Scratch and later in Maple, the latter using an implementation of turtle geometry that we developed long ago and has a much better graphic resolution (E. Roanes-Lozano and E. Roanes-Macías. An Implementation of “Turtle Graphics” in Maple V. MapleTech. Special Issue, 1994, 82-85). A video (in Spanish) of the Scratch version is available from the Instituto de Matemática Interdisciplinar (IMI) web page: https://www.ucm.es/imi/other-activities

Introduction

Surely you are uncomfortable being locked up at home, so I will try to justify that, although we are all looking forward going out, it is good not to meet your friends and family with whom you do not live.

I firstly need to mention a fractal is. A fractal is a geometric object whose structure is repeated at any scale. An example in nature is Romanesco broccoli, that you perhaps have eaten (you can search for images on the Internet). You can find a simple fractal in the following image (drawn with Maple):

Notice that each branch is divided into two branches, always forming the same angle and decreasing in size in the same proportion.

We can say that the tree in the previous image is of “depth 7” because there are 7 levels of branches.

It is quite easy to create this kind of drawing with the so called “turtle geometry” (with a recursive procedure, that is, a procedure that calls itself). Perhaps you have used Scratch programming language at school (or Logo, if you are older), which graphics are based in turtle geometry.

All drawings along these pages have been created with Maple. We can easily reform the code that generated the previous tree so that it has three, four, five,… branches at each level, instead of two.

But let’s begin with a tale that explains in a much simplified way how the spread of a disease works.

- o O o -

Let's suppose that a cat returns sick to Catland suffering from a very contagious disease and he meets his friends and family, since he has missed them so much.

We do not know very well how many cats each sick cat infects in average (before the order to STAY AT HOME arrives, as cats in Catland are very obedient and obey right away). Therefore, we’ll analyze different scenarios:

Each sick cat infects two other cats.
Each sick cat infects three other cats.
Each sick cat infects five other cats

1. Each Sick Cat Infects Two Cats

In all the figures that follow, the cat initially sick is in the center of the image. The infected cats are represented by a red square.

· Before everyone gets confined at home, it only takes time for that first sick cat to see his friends, but then confinement is ordered (depth 1)

As you can see, with the cat meeting his friends and family, we already have 3 sick cats.

· Before all cats confine themselves at home, the first cat meets his friends, and these in turn have time to meet their friends (depth 2)

In this case, the number of sick cats is 7.

· Before every cat is confined at home, there is time for the initially sick cat to meet his friends, for these to meet their friends, and for the latter (friends of the friends of the first sick cat) to meet their friends (depth 3).

There are already 15 sick cats...

· Depth 4: 31 sick cats.

· Depth 5: 63 sick cats.

Next we’ll see what would happen if each sick cat infected three cats, instead of two.

2. Every Sick Cat Infects Three Cats

· Now we speed up, as you’ve got the idea.

The first sick cat has infected three friends or family before confining himself at home. So there are 4 infected cats.

· If each of the recently infected cats in the previous figure have in turn contact with their friends and family, we move on to the following situation, with 13 sick cats:

· And if each of these 13 infected cats lives a normal life, the disease spreads even more, and we already have 40!

· At the next step we have 121 sick cats:

· And, if they keep seeing friends and family, there will be 364 sick cats (the image reminds of what is called a Sierpinski triangle):

4. Every Sick Cat Infects Five Cats

· In this case already at depth 2 we already have 31 sick cats.

5. Conclusion

This is an example of exponential growth. And the higher the number of cats infected by each sick cat, the worse the situation is.

Therefore, avoiding meeting friends and relatives that do not live with you is hard, but good for stopping the infection. So, it is hard, but I stay at home at the first stage too!

Monte Carlo Integration

by: radaar 202 Product: Maple

June 15 2020

1 1

Monte Carlo integration uses random sampling unlike classical techniques like the trapezoidal or Simpson's rule in evaluating the integration numerically.

Download MONTE_CARLO_INTEGRATION1.mw

Minimal Road Radius - solving a physics problem...

by: Scot Gould 1039 Product: Maple 2020

June 14 2020

3 2

The purpose of this document is:

a) to correct the physics that was used in the document "Minimal Road Radius for Highway Superelevation" recently submitted to the Maple Applications Center;

b) to confirm the values found in the manual for the American Association of State Highway and Transportation Officials (AASHTO) that engineers use to design and build these banked curves are physically sound.

c) to highlight the pedagogical value inherent in the Maple language to distinguish between assignment ( := ) and equivalence ( = );

d) but most importantly, to demonstrate the pedagogical value Maple has in thinking about solving a problem involving a physical process. Given Maple's symbolic mathematics capabilities, one can implement a top-down approach to the physics and the mathematics, working from the general principle to the specific example. This allows one to avoid the types of errors that occur when translating the problem into a bottom up approach, from specific values of the example to the general principle, an approach that is required by most other computational systems.

I hope that others are willing to continue to engage in discussions related to the pedagogical value of Maple beyond mathematics.

I was asked to post this document to both here and the Maple Applications Center

[Document edited for typos.]

Minimum_Road_Radius.mw

Solving multi-span Euler beams with the Method...

by: Rouben Rostamian 8456 Product: Maple 2020

June 08 2020

14 15

Maple's pdsolve() is quite capable of solving the PDE that describes the motion of a single-span Euler beam. As far as I have been able to ascertain, there is no obvious way of applying pdsolve() to solve multi-span beams. The worksheet attached to this post provides tools for solving multi-span Euler beams. Shown below are a few demos. The worksheet contains more demos.

A module for solving the Euler beam with the method of lines

beamsolve

The beamsolve proc solves a (possibly multi-span) Euler beam equation:

subject to initial and boundary conditions. The solution is the

transverse deflection of the beam at point at time , subject to the load
density (i.e., load per unit length) given by . The coefficient

is the beam's mass density (mass per unit length), is the Young's modulus of

the beam's material, and is the beam's cross-sectional moment of inertia

about the neutral axis. The figure below illustrates a 3-span beam (drawn in green)
supported on four supports, and loaded by a variable density load (drawn in gray)
which may vary in time. The objective is to determine the deformed shape of the
beam as a function of time.

The number of spans, their lengths, and the nature of the supports may be specified

as arguments to beamsolve.

In this worksheet we assume that are constants for simplicity. Since only
the product of the coefficients and enters the calculations, we lump the two

together into a single variable which we indicate with the two-letter symbol .

Commonly, is referred to as the beam's rigidity.

The PDE needs to be supplied with boundary conditions, two at each end, each

condition prescribing a value (possibly time-dependent) for one of
(that's 36 possible combinations!) where I have used subscripts to indicate

derivatives. Thus, for a single-span beam of length , the following is an admissible

set of boundary conditions:
. (Oops, coorection, that last
condition was meant to be u_xxx(L,t) = sin t.)
Additionally, the PDE needs to be supplied with initial conditions which express

the initial displacement and the initial velocity:

The PDE is solved through the Method of Lines. Thus, each span is subdivided into

subintervals and the PDE's spatial derivatives are approximated through finite differences.

The time derivatives, however, are not discretized. This reduces the PDE into a set of

ODEs which are solved with Maple's dsolve().

Calling sequence:

beamsolve(, options)

Parameters:

: List of span lengths, in order from left to right, as in .

The number of subintervals in the shortest span (for the finite difference approximation)

Notes:

•	It is assumed that the spans are laid back-to-back along the axis, with the left end of the overall beam at .

•

The interior supports, that is, those supports where any two spans meet, are assumed
to be of the so-called simple type. A simple support is immobile and it doesn't exert
a bending moment on the beam. Supports at the far left and far right of the beam can
be of general type; see the BC_left and BC_right options below.

•	If the beam consists of a single span, then the argument may be entered as a number rather than as a list. That is, is equivalent to .

Options:

All options are of the form option_name=value, and have built-in default values.

Only options that are other than the defaults need be specified.

rho: the beam's (constant) mass density per unit length (default: rho = 1)

EI: the beam's (constant) rigidity (default: EI = 1)

T: solve the PDE over the time interval (default: T = 1)

        F: an expression in and that describes the applied force   (default: F = 0)
        IC: the list of the initial displacement and velocity, as
                expressions in (default: IC = [0,0])

        BC_left: a list consisting of a pair of boundary conditions at the left end of
                the overall (possibly multi-span beam. These can be any two of
                . The right-hand sides of these equations

can be any expression in . The left-hand sides should be entered literally as indicated.

If a right-hand side is omitted, it is taken to be zero. (default: BC_left = which

corresponds to a simple support).

BC_right: like BC_left, but for the right end of the overall beam (default: BC_right =

The returned module:

A call to beamsolve returns a module which presents three methods. The methods are:

plot (t, refine=val, options)

plots the solution at time . If the discretization in the direction

is too coarse and the graph looks non-smooth, the refine option

(default: refine=1) may be specified to smooth out the graph by introducing

val number of intermediate points within each discretized subinterval.

All other options are assumed to be plot options and are passed to plots:-display.

plot3d (m=val, options)

plots the surface . The optional specification requests

a grid consisting of subintervals in the time direction (default:

Note that this grid is for plotting purposes only; the solution is computed

as a continuous (not discrete) function of time. All other options are assumed

to be plot3d options and are passed to plots:-display.

animate (frames=val, refine=val, options)

produces an animation of the beam's motion. The frames option (default = 50)

                specifies the number of animation frames. The refine option is passed to plot
                (see the description above. All other options are assumed to be plot options and
                are passed to plots:-display.

Note:

In specifying the boundary conditions, the following reminder can be helpful. If the beam

is considered to be horizontal, then is the vertical displacement, is the slope,

is the bending moment, and is the transverse shear force.

>

A single-span simply-supported beam with initial velocity

The function is an exact solution of a simply supported beam with

The solution is periodic in time with period .

>	sol := beamsolve(1, 25, 'T'=2/Pi, 'IC'=[0, Pi^2sin(Pix)]): sol:-animate(size=[600,250]);

The initial condition does not lead to a separable form, and

therefore the motion is more complex.

>	sol := beamsolve(1, 25, 'T'=2/Pi, 'IC'=[0, 1]): sol:-animate(frames=200, size=[600,250]);

>

A single-span cantilever beam

A cantilever beam with initial condition where is the
first eigenmode of its free vibration (calculated in another spreadsheet). The motion is
periodic in time, with period

>	g := 0.5cos(1.875104069x) - 0.5cosh(1.875104069x) - 0.3670477570sin(1.875104069x) + 0.3670477570sinh(1.875104069x): sol := beamsolve(1, 25, 'T'=1.787018777, 'BC_left'=[u,u_x], 'BC_right'=[u_xx,u_xxx], 'IC'=[g, 0]): sol:-animate(size=[600,250]);

If the initial condition is not an eigenmode, then the solution is rather chaotic.

>	sol := beamsolve(1, 25, 'T'=3.57, 'BC_left'=[u,u_x], 'BC_right'=[u_xx,u_xxx], 'IC'=[-x^2, 0]): sol:-animate(size=[600,250], frames=100);

>

A single-span cantilever beam with a weight hanging from its free end

>	sol := beamsolve(1, 25, 'T'=3.57, 'BC_left'=[u,u_x], 'BC_right'=[u_xx,u_xxx=1]): sol:-animate(size=[600,250], frames=100);

>

A single-span cantilever beam with oscillating support

>	sol := beamsolve(1, 25, 'T'=Pi, 'BC_left'=[u=0.1sin(10t),u_x], 'BC_right'=[u_xx,u_xxx]): sol:-animate(size=[600,250], frames=100);

>

A dual-span simply-supported beam with moving load

Load moves across a dual-span beam.

The beam continues oscillating after the load leaves.

>

d := 0.4: T := 4: nframes := 100:
myload := - max(0, -6*(x - t)*(d + x - t)/d^3):
sol := beamsolve([1,1], 20, 'T'=T, 'F'=myload):
sol:-animate(frames=nframes):
plots:-animate(plot, [2e-3*myload(x,t), x=0..2, thickness=1, filled=[color="Green"]], t=0..T, frames=nframes):
plots:-display([%%,%], size=[600,250]);

>

A triple-span simply-supported beam with moving load

Load moves across a triple-span beam.

The beam continues oscillating after the load leaves.

>

d := 0.4: T := 6: nframes := 100:
myload := - max(0, -6*(x - t)*(d + x - t)/d^3):
sol := beamsolve([1,1,1], 20, 'T'=T, 'F'=myload):
sol:-plot3d(m=50);
sol:-animate(frames=nframes):
plots:-animate(plot, [2e-3*myload(x,t), x=0..3, thickness=1, filled=[color="Green"]], t=0..T, frames=nframes):
plots:-display([%%,%], size=[600,250]);

>

z3d;

>

A triple-span beam, moving load falling off the cantilever end

In this demo the load move across a multi-span beam with a cantilever section at the right.

As it skips past the cantilever end, the beam snaps back violently.

>

d := 0.4: T := 8: nframes := 200:
myload := - max(0, -6*(x - t/2)*(d + x - t/2)/d^3):
sol := beamsolve([1,1,1/2], 10, 'T'=T, 'F'=myload, BC_right=[u_xx, u_xxx]):
sol:-animate(frames=nframes):
plots:-animate(plot, [1e-2*myload(x,t), x=0..3, thickness=1, filled=[color="Green"]], t=0..T, frames=nframes):
plots:-display([%%,%], size=[600,250]);

>

Download worksheet: euler-beam-with-method-of-lines.mw

Scalable Subplots within Figure

by: Bland3 30 Product: Maple 2019

May 29 2020

0 3

This is my attempt to produce a subplot within a larger plot for magnifying data in a small region, and putting that subplot into the white space of the figure.
Based on the questions: How to insert a plot into another plot? and Inset figure in Maple, I wrote a couple of procedures that create sub-plots and allow the user to place the subplot window as he/she chooses. This avoids the graininess issues mentioned by acer in the second link (and experienced by me).

So far, I only have this completed for point plots, but using acer's method of piecewise functions posted in the plotin2b.mw of the second article, with the subplot function being defined only if it satisfies your conditions, would allow the subplot generating procedure to be generalized easily enough. But the data I'm working with all point plots, so that's the example here.

The basic idea is to use one procedure to create boxes, make tickmarks on the expanded region, and make tickmark labels, combine all of those into one graph. Then create scaled and shifted versions of the data series, then make graphs of those. Lastly, combine them all into one picture.

Hope this helps someone who has to do the same.

Mapleprimes isn't inserting the contents, but here is the download of the file: SubPlotBoxesandVectorDataSeries.mw

Partitioning a cube into pyramids: An animation

by: Rouben Rostamian 8456 Product: Maple

May 21 2020

7 3

Here is a little cute demo that shows how a cube may be paritioned into three congruent pyramids. This was inspired by a Mathematica demo that I found in the web but I think this one's better :-)

A Cube as a union of three right pyramids

Here is an animated demo of the well-known fact that a cube may be partitioned

into three congruent right pyramids.

2020-05-21

>

restart;

>	with(plots):

>	with(plottools):

A proc to plot a general polyhedron.
list of vertices
list of faces

An entry in describes a face made of the vertices etc.

>	polyhedron := proc(V::list, F::list) seq(plottools:-polygon([seq( V[F[i][j]], j=1..nops(F[i]))]), i=1..nops(F)); plots:-display(%); end proc:

Define the vertices and faces of a pyramid:

>	v := [[0,0,0],[1,0,0],[1,1,0],[0,1,0],[0,0,1]]; f := [ [1,2,3,4], [5,2,3], [5,3,4], [1,5,4], [1,2,5] ];

Build three such pyramids:

>	P1 := polyhedron(v, f): P2 := reflect(P1, [[1,0,0],[1,1,0],[1,0,1]]): P3 := reflect(P1, [[0,1,0],[1,1,0],[0,1,1]]):

This is what we have so far:

>	display(P1,P2,P3, scaling=constrained);

Define an animation frame. The parameter goes from 0 to 1.

Any extra options are assumed to be plot3d options and are

passed to plots:-display.

>	frame := proc(t) plots:-display( P1, rotate(P2, Pi/2t, [[1,1,0],[1,0,0]]), rotate(P3, Pi/2t, [[0,1,0],[1,1,0]]), color=["Red", "Green", "Blue"], _rest); end proc:

Animate:

>	display(frame(0) $40, seq(frame(t), t=0..1, 0.01), frame(1) $40, insequence, scaling=constrained, axes=none, orientation=[45,0,120], viewpoint=circleleft);

Download square-partitioned-into-pyramids.mw

The equations of motion in curvilinear coordinates,...

by: ecterrab 14535 Product: Maple

May 21 2020

5 0

The equations of motion in curvilinear coordinates, tensor notation and Coriolis force

The formulation of the equations of motion of a particle is simple in Cartesian coordinates using vector notation. However, depending on the problem, for example when describing the motion of a particle as seen from a non-inertial system of references (e.g. a rotating planet, like earth), there is advantage in using curvilinear coordinates and also tensor notation. When the particle's movement is observed from such a rotating referential, we also see "acceleration" that is not due to any force but to the fact that the referential itself is accelerated. The material below discusses and formulates these topics, and derives the expression for the Coriolis and centripetal force in cylindrical coordinates as seen from a rotating system of references.

The computation below is reproducible in Maple 2020 using the Maplesoft Physics Updates v.681 or newer.

Vector notation

Generally speaking the equations of motion of a particle are easy to formulate: the position vector is a function of time, the velocity is its first derivative and the acceleration is its second derivative. For instance, in Cartesian coordinates

>

(1)

>

(2)

>

(3)

Newton's 2^nd law, that in an inertial system of references when there is force there is acceleration and viceversa, is then given by

>

(4)

where represents the total force acting on the particle. This vectorial equation represents three second order differential equations which, for given initial conditions, can be integrated to arrive at a closed form expression for as a function of .

Tensor notation

In Cartesian coordinates, the tensorial form of the equations (4) is also straightforward. In a flat spacetime - Galilean system of references - the three space coordinates form a 3D tensor, and so does its first derivate and the second one. Set the spacetime to be 3-dimensional and Euclidean and use lowercaselatin indices for the corresponding tensors

>

Physics:-g_[mu, nu] = Matrix(%id = 18446744078329083054)

(5)

The position, velocity and acceleration vectors are expressed in tensor notation as done in (1), (2) and (3)

>

(6)

>

(7)

>

(8)

Setting a tensor to represent the three Cartesian components of the force

>

${Physics:-Dgamma[a], F[j], Physics:-Psigma[a], Physics:-d_[a], Physics:-g_[a, b], Physics:-LeviCivita[a, b, c], Physics:-SpaceTimeVector[a](X)}$

(9)

Newton's 2^nd law (4), now expressed in tensorial notation, is given by

>

(10)

The three differential equations behind this tensorial form of (4) are as expected

>

${F__x(t) = m*(diff(diff(x(t), t), t)), F__y(t) = m*(diff(diff(y(t), t), t)), F__z(t) = m*(diff(diff(z(t), t), t))}$

(11)

Things are straightforward in Cartesian coordinates because the components of the line element are exactly the components of the tensor

>

(12)

and so, the form factors (see related Mapleprimes post) are all equal to 1.

In the case of no external forces, and the equations of motion, whose solution are the trajectory, can be formulated as the path of minimal length between two points, a geodesic. In the case under consideration, because the spacetime is flat (Galilean) those two points lie on a plane, we are talking about Euclidean geometry, that information is encoded in the metric (the 3x3 identity matrix (5)), and the geodesic is a straight line. The differential equations of this geodesic are thus the equations of motion (11) with , and can be computed using Geodesics

>

(13)

Curvilinear coordinates

Vector notation

The form of these equations in the case of curvilinear coordinates, for example in cylindrical or spherical variables, is obtained performing a change of coordinates.

>

(14)

This change keeps the z axis unchanged, so the corresponding unit vector remains unchanged.

Changing the basis and coordinates used to represent the position vector , it becomes

>

(15)

where since in (1) the coordinates () depend on t, in (15), not just and but also the unit vector depends on t. The velocity is computed as usual, differentiating

>

(16)

The second term is due to the dependency of on the coordinate together with the chain rule diff(`#mover(mi("ρ",fontstyle = "normal"),mo("&and;"))`(t), t) = (diff(`#mover(mi("ρ",fontstyle = "normal"),mo("&and;"))`, phi))*(diff(phi(t), t)) and (diff(`#mover(mi("ρ",fontstyle = "normal"),mo("&and;"))`, phi))*(diff(phi(t), t)) = `#mover(mi("φ",fontstyle = "normal"),mo("&and;"))`(t)*(diff(phi(t), t)) . The dependency of curvilinear unit vectors on the coordinates is automatically taken into account when differentiating due to the Setup setting geometricdifferentiation = true.

For the acceleration,

>

(17)

where the term involving comes from differentiating in (16) taking into account the dependency of on the coordinate This result can also be obtained by directly changing variables in the acceleration , in equation (3)

>

(18)

Newton's 2^nd law becomes

>

(19)

In the absence of external forces, equating to 0 the vector components (coefficients of the unit vectors) of the acceleration we get the system of differential equations in cylindrical coordinates whose solution is the trajectory of the particle expressed in the (

>

$`~`[`=`]({coeffs(rhs(diff(diff(r_(t), t), t) = _rho(t)*(diff(diff(rho(t), t), t)-rho(t)*(diff(phi(t), t))^2)+_phi(t)*(2*(diff(rho(t), t))*(diff(phi(t), t))+rho(t)*(diff(diff(phi(t), t), t)))+(diff(diff(z(t), t), t))*_k), [`#mover(mi("ρ",fontstyle = "normal"),mo("&and;"))`(t), `#mover(mi("φ",fontstyle = "normal"),mo("&and;"))`(t), `#mover(mi("k"),mo("&and;"))`])}, 0)$

${2*(diff(rho(t), t))*(diff(phi(t), t))+rho(t)*(diff(diff(phi(t), t), t)) = 0, diff(diff(rho(t), t), t)-rho(t)*(diff(phi(t), t))^2 = 0, diff(diff(z(t), t), t) = 0}$

(20)

>	$solve({2(diff(rho(t), t))(diff(phi(t), t))+rho(t)(diff(diff(phi(t), t), t)) = 0, diff(diff(rho(t), t), t)-rho(t)(diff(phi(t), t))^2 = 0, diff(diff(z(t), t), t) = 0}, {diff(phi(t), t, t), diff(rho(t), t, t), diff(z(t), t, t)})$

${diff(diff(phi(t), t), t) = -2*(diff(rho(t), t))*(diff(phi(t), t))/rho(t), diff(diff(rho(t), t), t) = rho(t)*(diff(phi(t), t))^2, diff(diff(z(t), t), t) = 0}$

(21)

In this formulation (21) with , although , no acceleration in the direction, is naturally expected, the same cannot be said about the other two equations for and . Those two equations are discussed below under Coriolis and Centripetal forces. The key observation at this point, however, is that the right-hand sides of both unexpected equations involve , rotation around the z axis.

Tensor notation

The same equations (19) and (21) result when using tensor notation. For that purpose, one can transform the position, velocity and acceleration tensors (6), (7), (8), but since they are expressed as functions of a parameter (the time), it is simpler to transform only the underlying metric using TransformCoordinates. That has the advantage that all the geometrical subtleties of curvilinear coordinates, like scale factors and dependency of unit vectors on curvilinear coordinates, get automatically, very succinctly, encoded in the metric:

>

Physics:-g_[a, b] = Matrix(%id = 18446744078263848958)

(22)

The computation of geodesics assumes that the coordinates () depend on a parameter. That parameter is passed as the first argument to Geodesics

>

(23)

These equations of motion (23) are the same as the equations (21) computed using standard vector notation, differentiating and taking into account the dependency of curvilinear unit vectors on the curvilinear coordinates in (16) and (17). One of the interesting features of computing with tensors is, as said, that all those geometrical algebraic subtleties of curvilinear coordinates are automatically encoded in the metric (22).

To understand how are the geodesic equations computed in one go in (23), one can perform the calculation in steps:

1.	Make be a function of directly in the metric

2.	Compute - not the final form of the equations (23) - but the intermediate form expressing the geodesic equation using tensor notation, in terms of Christoffel symbols

3.	Compute the components of that tensorial equation for the geodesic (using TensorArray)

For step 1, we have

>

Physics:-g_[a, b] = Matrix(%id = 18446744078354237910)

(24)

Set this metric where

>

Physics:-g_[a, b] = Matrix(%id = 18446744078342604430)

(25)

Step 2, the geodesic equations in tensor notation with the coordinates depending on the time t are computed passing the optional argument tensornotation

>

(26)

Step 3: compute the components of this tensorial equation

>

(27)

These are the same equations (23).

Having the tensorial equation (26) is also useful to formulate the equations of motion in tensorial form in the presence of force. For that purpose, redefine the contravariant tensor to represent the force in the cylindrical basis

>

${Physics:-D_[a], Physics:-Dgamma[a], F[j], Physics:-Psigma[a], Physics:-Ricci[a, b], Physics:-Riemann[a, b, c, d], Physics:-Weyl[a, b, c, d], Physics:-d_[a], Physics:-g_[a, b], Physics:-Christoffel[a, b, c], Physics:-Einstein[a, b], Physics:-LeviCivita[a, b, c], Physics:-SpaceTimeVector[a](X)}$

(28)

Newton's 2^nd law (19)

>

(29)

now using tensorial notation, becomes

>

(30)

>

Array(%id = 18446744078329063774)

(31)

where we recall (see related Mapleprimes post) that to obtain the vector components entering from these tensor components we need to multiply the latter by the scale factors (), the component of in the direction of is given by .

Coriolis force and centripetal force

After changing variables the position vector of the particle got expressed in (15) as

A distinction needs to be made here, according to whether the unit vector depends or not on the time , the former being the general case. When is a constant, the value of the coordinate - the angle between and the x axis - does not change, there is no rotation around the z axis. On the other hand, when , the orientation of and so the coordinate changes with time, so either the force acting on the particle has a component in the direction that produces rotation around the z axis, or the system of references - itself - is rotating around the z axis.

Likewise, the expression (15) can represent the position vector measured in the original Galilean (inertial) system of references, where a force is producing rotation around the z axis, or it can represent the position of the particle measured in a rotating, non-inertial system references. Hence the transformation (14) can also be interpreted in two different ways, as representing a choice of different functions (generalized coordinates) to represent the position of the particle in the original inertial system of references, or it can represent a transformation from an inertial to another rotating, non-inertial, system of references.

This equivalence between the trajectory of a particle subject to an external force, as observed in an inertial system of references, and the same trajectory observed from a non-inertial accelerated system of references where there is no external force, actually at the root of the formulation of general relativity, is also well known in classical mechanics. The (apparent) forces observed in the rotating non-inertial system of references, due only to its acceleration, are called Coriolis and centripetal forces.

To see that the equations

diff(rho(t), t, t) = (diff(phi(t), t))^2*rho(t), diff(phi(t), t, t) = -2*(diff(phi(t), t))*(diff(rho(t), t))/rho(t)

that appeared in (27) when in the inertial system of references , are related to the Coriolis and centripetal forces in the non-inertial referencial, following [1] introduce a vector representing the rotation of that referencial around the z axis (when, in the inertial system of references, the non-inertial system rotates clockwise, in the non-inertial system increases in value in the anti-clockwise direction)

>

(32)

According to [1], (39.7), the acceleration in the inertial system is expressed in terms of the quantities in the non-inertial rotating system by the sum of the following three vectorial terms.

First, the components of the acceleration measured in the non-inertial system are given by the second derivatives of the coordinates () multiplied by the scale factors, which in this case are given by () (see this post in Mapleprimes)

>

(33)

Second, the Coriolis force divided by the mass, by definition given by

>

(34)

Third the centripetal force divided by the mass, defined by

>

(35)

Adding these three terms,

>

a_(t)+2*Physics:-Vectors:-`&x`(diff(r_(t), t), omega_)+Physics:-Vectors:-`&x`(omega_, Physics:-Vectors:-`&x`(r_(t), omega_)) = _rho(t)*(diff(diff(rho(t), t), t)-rho(t)*(diff(phi(t), t))^2)+_phi(t)*(2*(diff(rho(t), t))*(diff(phi(t), t))+rho(t)*(diff(diff(phi(t), t), t)))+(diff(diff(z(t), t), t))*_k

(36)

So that

>

(37)

and where the right-hand side of (36) is, actually, the result computed lines above in (18)

>

(38)

>

(39)

From (37), in the absence of external forces and so the acceleration measured in the rotating system is given by the sum of the Coriolis and centripetal accelerations

>

(40)

In other words: in the absence of external forces, the acceleration of a particle observed in a rotating (non-inertial) system of references is not zero.

Expressing this equation (40) in terms of () we get