Posts - MaplePrimes

Eigenfunctions of a multi-span Euler-Bernoulli...

by: Rouben Rostamian 8596 Product: Maple 2020

August 19 2020

5 6

Question about deflection and vibration of beams occur with some regularity in this forum. Search for "beam" to see several pages of hits.

In this post I present a general approach to calculating the vibrational modes of a beam that applies to both single-span and multi-span beams. The code is not perfectly polished, but it is sufficiently documented to enable the interested user to modify/extend it as needed.

Vibrational modes of multi-span Euler-Bernoulli beams

through Krylov-Dunction functions

Rouben Rostamian
2020-07-19

>

restart;

Note: Maple defines the imaginary unit . We want to use the
symbol as the beam's cross-sectional moment of inertia.
Therefore we redefine the imaginary unit (for which we have no

use) as and free up the symbol for our use.

>	interface(imaginaryunit=II):

>	with(LinearAlgebra):

>

The Euler-Bernoulli beam equation
"rho*A*((&PartialD;)^2u)/((&PartialD;)^( )t^2)+E*I*((&PartialD;)^(4)u)/((&PartialD;)^( )x^(4))=0" .

We wish to determine the natural modes of vibration of

a possibly multi-span Euler-Bernoulli beam.

Separate the variables by setting . We get
-
"(rho*A)/(E*I)*(T ' ')/(T)=(X^((4)))/(X)=mu^(4) "
whence
.

Let omega = sqrt(I*E/(rho*A))*mu^2 . Then

and

The idea behind the Krylov-Duncan technique is to express

in terms an alternative (and equivalent) set of basis
functions through ,, as
,

where the functions through are defined in the next section.

In some literature the symbols , are used for these

functions but I find it more sensible to use the indexed function

notation.

The Krylov-Duncan approach is particularly effective in formulating
and finding a multi-span beam's natural modes of vibration.

>

The Krylov-Duncan functions

The K[i](x) defined by this proc evaluates to the th

Krylov-Duncan function.

Normally the index will be in the set, however the proc is

set up to accept any integer index (positive or negative). The proc evaluates

the index modulo 4 to bring the index into the set . For

instance, K[5](x) and K[-3](x)i are equivalent to K[1](x) .

>

K := proc(x)
        local n := op(procname);

        if not type(n, integer) then
                return 'procname'(args);
        else
                n := 1 + modp(n-1,4); # reduce n modulo 4
        end if;

        if n=1 then
                (cosh(x) + cos(x))/2;
        elif n=2 then
                (sinh(x) + sin(x))/2;
        elif n=3 then
                 (cosh(x) - cos(x))/2;
        elif n=4 then
                (sinh(x) - sin(x))/2;
        else
                error "shouldn't be here!";
        end if;

end proc:

Here are the Krylov-Duncan basis functions:

>	seq(print(cat(`K__`,i)(x) = K[i](x)), i=1..4);

and here is what they look like. All grow exponentially for large
but are significantly different near the origin.

>	plot([K[i](x) $i=1..4], x=-Pi..Pi, color=["red","Green","blue","cyan"], thickness=2, legend=['K[1](x)', 'K[2](x)', 'K[3](x)', 'K[4](x)']);

The cyclic property of the derivatives:
We have . Let's verify that:

>	diff(K[i](x),x) - K[i-1](x) $i=1..4;

The fourth derivative of each function equals itself. This is a consequence of the cyclic property:

>	diff(K[i](x), x$4) - K[i](x) $ i=1..4;

The essential property of the Krylov-Duncan basis function is that their

zeroth through third derivatives at form a basis for :

>	seq((D@@n)(K[1])(0), n=0..3); seq((D@@n)(K[2])(0), n=0..3); seq((D@@n)(K[3])(0), n=0..3); seq((D@@n)(K[4])(0), n=0..3);

As noted earlier, in the case of a single-span beam, the modal shapes

are expressed as
.

Then, due to the cyclic property of the derivatives of the Krylov-Duncan

functions, we see that:
.
.
.
Let us note, in particular, that
,
,
,
.

>

A general approach for solving multi-span beams

In a multi-span beam, we write for the deflection of the th span, where

and where is the span's length. The coordinate indicates the

location within the span, with corresponding to the span's left endpoint.

Thus, each span has its own coordinate system.

We assume that the interface of the two adjoining spans is supported on springs

which (a) resist transverse displacement proportional to the displacement (constant of

proportionality of (d for displacement), and (b) resist rotation proportional to the
slope (constant of proportionality of (t for torsion or twist). The spans are numbered

from left to right. The interface conditions between spans and +1 are

1.	The displacements at the interface match: .

2.	The slopes at the interface match (0).

3.	The difference of the moments just to the left and just to the right of the support is due to the torque exerted by the torsional spring:

4.	The difference of the shear forces just to the left and just to the right of the support is due to the force exerted by the linear spring:

The special case of a pinned support corresponds to and .
In that case, condition 3 above implies that ,
and condition 4 implies that

Let us write the displacements and in terms of the Krylov-Duncan

functions as:

Then applying the cyclic properties of the Krylov-Duncan functions described

earlier, the four interface conditions translate to the following system of four
equations involving the eight coefficients .

,

which we write as a matrix equation
$(Matrix(4, 8, {(1, 1) = K__1(mu*L__i), (1, 2) = K__2(mu*L__i), (1, 3) = K__3(mu*L__i), (1, 4) = K__4(mu*L__i), (1, 5) = -1, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (2, 1) = K__4(mu*L__i), (2, 2) = K__1(mu*L__i), (2, 3) = K__2(mu*L__i), (2, 4) = K__3(mu*L__i), (2, 5) = 0, (2, 6) = -1, (2, 7) = 0, (2, 8) = 0, (3, 1) = K__3(mu*L__i), (3, 2) = K__4(mu*L__i), (3, 3) = K__1(mu*L__i), (3, 4) = K__2(mu*L__i), (3, 5) = 0, (3, 6) = -I*k__t/(mu*E), (3, 7) = -1, (3, 8) = 0, (4, 1) = K__2(mu*L__i), (4, 2) = K__3(mu*L__i), (4, 3) = K__4(mu*L__i), (4, 4) = K__1(mu*L__i), (4, 5) = -I*k__d/(mu^3*E), (4, 6) = 0, (4, 7) = 0, (4, 8) = -1}))*(Vector(8, {(1) = `a__i,1`, (2) = `a__i,2`, (3) = `a__i,3`, (4) = `a__i,4`, (5) = `a__i+1,1`, (6) = `a__i+1,2`, (7) = `a__i+1,3`, (8) = `a__i+1,4`})) = (Vector(8, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0}))$ .

That coefficient matrix plays a central role in solving

for modal shapes of multi-span beams. Let's call it .

Note that the value of enters that matrix only in combinations with
and . Therefore we introduce the new symbols

, .

The following proc generates the matrix . The parameters and

are optional and are assigned the default values of infinity and zero, which

corresponds to a pinned support.

The % sign in front of each Krylov function makes the function inert, that is, it
prevents it from expanding into trig functions. This is so that we can

see, visually, what our expressions look like in terms of the functions. To

force the evaluation of those inert function, we will apply Maple's value function,

as seen in the subsequent demos.

>

M_interface := proc(mu, L, {Kd:=infinity, Kt:=0})
        local row1, row2, row3, row4;
        row1 := %K[1](mu*L), %K[2](mu*L), %K[3](mu*L), %K[4](mu*L), -1, 0, 0, 0;
        row2 := %K[4](mu*L), %K[1](mu*L), %K[2](mu*L), %K[3](mu*L), 0, -1, 0, 0;
        row3 := %K[3](mu*L), %K[4](mu*L), %K[1](mu*L), %K[2](mu*L), 0, Kt/mu, -1, 0;
        if Kd = infinity then
                row4 := 0, 0, 0, 0, 1, 0, 0, 0 ;
        else
                row4 := %K[2](mu*L), %K[3](mu*L), %K[4](mu*L), %K[1](mu*L), Kd/mu^3, 0, 0, -1;
        end if:
                return < <row1> | <row2> | <row3> | <row4> >^+;
end proc:

Here is the interface matrix for a pinned support:

>	M_interface(mu, L);

And here is the interface matrix for a general springy support:

>	M_interface(mu, L, 'Kd'=a, 'Kt'=b);

$Matrix(4, 8, {(1, 1) = %K[1](L*mu), (1, 2) = %K[2](L*mu), (1, 3) = %K[3](L*mu), (1, 4) = %K[4](L*mu), (1, 5) = -1, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (2, 1) = %K[4](L*mu), (2, 2) = %K[1](L*mu), (2, 3) = %K[2](L*mu), (2, 4) = %K[3](L*mu), (2, 5) = 0, (2, 6) = -1, (2, 7) = 0, (2, 8) = 0, (3, 1) = %K[3](L*mu), (3, 2) = %K[4](L*mu), (3, 3) = %K[1](L*mu), (3, 4) = %K[2](L*mu), (3, 5) = 0, (3, 6) = b/mu, (3, 7) = -1, (3, 8) = 0, (4, 1) = %K[2](L*mu), (4, 2) = %K[3](L*mu), (4, 3) = %K[4](L*mu), (4, 4) = %K[1](L*mu), (4, 5) = a/mu^3, (4, 6) = 0, (4, 7) = 0, (4, 8) = -1})$

Note: In Maple's Java interface, inert quantities are shown in gray.

Note: The in this matrix is the length of the span to the left of the interface.
Recall that it is , not , in the derivation that leads to that matrix.

In a beam consisting of spans, we write the th span's deflection as

Solving the beam amounts to determining the unknowns

which we order as

At each of the interface supports we have a set of four equations as derived
above, for a total of equations. Additionally, we have four user-supplied

boundary conditions -- two at the extreme left and two at the extreme right of the

overall beam. Thus, altogether we have equations which then we solve for the

unknown coefficients .

The user-supplied boundary conditions at the left end are two equations, each in the
form of a linear combination of the coefficients . We write for the
coefficient matrix of that set of equations. Similarly, the user-supplied boundary
conditions at the right end are two equations, each in the form of a linear combination
of the coefficients . We write for the coefficient matrix of
that set of equations. Putting these equations together with those obtained at the interfaces,

we get a linear set of equations represented by a matrix which can be assembled

easily from the matrices , and . In the case of a 4-span beam the

assembled matrix looks like this:

The pattern generalizes to any number of spans in the obvious way.

For future use, here we record a few frequently occurring and matrices.

>	M_left_pinned := < 1, 0, 0, 0; 0, 0, 1, 0 >;

Matrix(2, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 1, (2, 4) = 0})

>	M_right_pinned := (mu,L) -> < %K[1](muL), %K[2](muL), %K[3](muL), %K[4](muL); %K[3](muL), %K[4](muL), %K[1](muL), %K[2](muL) >;

>	M_left_clamped := < 1, 0, 0, 0; 0, 1, 0, 0 >;

Matrix(2, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0})

>	M_right_clamped := (mu,L) -> < %K[1](muL), %K[2](muL), %K[3](muL), %K[4](muL); %K[4](muL), %K[1](muL), %K[2](muL), %K[3](muL) >;

>	M_left_free := < 0, 0, 1, 0; 0, 0, 0, 1 >;

Matrix(2, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 1, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 1})

>	M_right_free := (mu,L) -> < %K[3](muL), %K[4](muL), %K[1](muL), %K[2](muL); %K[2](muL), %K[3](muL), %K[4](muL), %K[1](muL) >;

The following proc builds the overall matrix in the general case. It takes
two or three arguments. The first two arguments are the matrices
which are called and in the discussion above. If the beam
consists of a single span, that's all the information that need be supplied.
There is no need for the third argument.

In the case of a multi-span beam, in the third argument we supply the
list of the interface matrices , as in , listed in order
of the supports, from left to right. An empty list is also
acceptable and is interpreted as having no internal supports,
i.e., a single-span beam.

>

build_matrix := proc(left_bc::Matrix(2,4), right_bc::Matrix(2,4), interface_matrices::list)
        local N, n, i, M;

        # n is the number of internal supports
        n := 0;

        # adjust n if a third argument is supplied
        if _npassed = 3 then
                n := nops(interface_matrices);
                if n > 0 then
                        for i from 1 to n do
                                if not type(interface_matrices[i], 'Matrix(4,8)') then
                                        error "expected a 4x8 matrix for element %1 in the list of interface matrices", i;
                                end if;
                        end do;
                end if;
        end if;

        N := n + 1;                     # number of spans

        M := Matrix(4*N);
        M[1..2, 1..4] := left_bc;
        for i from 1 to n do
                M[4*i-1..4*i+2, 4*i-3..4*i+4] := interface_matrices[i];
        end do;
        M[4*N-1..4*N, 4*N-3..4*N] := right_bc;

        return M;
end proc:

For instance, for a single-span cantilever beam of length we get the following matrix:

>	build_matrix(M_left_clamped, M_right_free(mu,L));

Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, (3, 1) = %K[3](L*mu), (3, 2) = %K[4](L*mu), (3, 3) = %K[1](L*mu), (3, 4) = %K[2](L*mu), (4, 1) = %K[2](L*mu), (4, 2) = %K[3](L*mu), (4, 3) = %K[4](L*mu), (4, 4) = %K[1](L*mu)})

For a two-span beam with with span lengths of and , and all three
supports pinned, we get the following matrix:

>	build_matrix(M_left_pinned, M_right_pinned(mu,L[2]), [M_interface(mu, L[1])]);

Matrix(%id = 18446884696906262398)

The matrix represents a homogeneous linear system (i.e., the right-hand side vector

is zero.) To obtain a nonzero solution, we set the determinant of equal to zero.

That gives us a generally transcendental equation in the single unknown . Normally

the equation has infinitely many solutions. We call these

Remark: In the special case of pinned supports at the interfaces, that is, when
, the matrix depends only on the span lengths .
It is independent of the parameters that enter the Euler-Bernoulli
equations. The frequencies `ω__n` = sqrt(I*E/(rho*A))*`μ__n`^2 , however, depend on those parameters.

This proc plots the calculated modal shape corresponding to the eigenvalue .
The params argument is a set of equations which define the numerical values

of all the parameters that enter the problem's description, such as the span

lengths.

It is assumed that in a multi-span beam, the span lengths are named etc.,
and in a single-span beam, the length is named .

>

plot_beam := proc(M::Matrix,mu::realcons, params::set)
        local null_space, N, a_vals, i, j, A, B, P;
        eval(M, params);
        eval(%, :-mu=mu);
        value(%); #print(%);
        null_space := NullSpace(%); #print(%);
        if nops(null_space) <> 1 then
                error "Calculation failed. Increasing Digits and try again";
        end if;

        N := Dimension(M)[1]/4; # number of spans
        a_vals := convert([seq(seq(a[i,j], j=1..4), i=1..N)] =~ null_space[1], list);

        if N = 1 then
                eval(add(a[1,j]*K[j](mu*x), j=1..4), a_vals);
                P[1] := plot(%, x=0..eval(L,params));
        else
                A := 0;
                B := 0;
                for i from 1 to N do
                        B := A + eval(L[i], params);
                        eval(add(a[i,j]*K[j](mu*x), j=1..4), a_vals);
                        eval(%, x=x-A):
                        P[i] := plot(%, x=A..B);
                        A := B;
                end do;
                unassign('i');
        end if;
        plots:-display([P[i] $i=1..N]);

end proc:

>

A single-span pinned-pinned beam

Here we calculate the natural modes of vibration of a single span

beam, pinned at both ends. The modes are of the form

The matrix is:

>	M := build_matrix(M_left_pinned, M_right_pinned(mu,L));

Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 1, (2, 4) = 0, (3, 1) = %K[1](L*mu), (3, 2) = %K[2](L*mu), (3, 3) = %K[3](L*mu), (3, 4) = %K[4](L*mu), (4, 1) = %K[3](L*mu), (4, 2) = %K[4](L*mu), (4, 3) = %K[1](L*mu), (4, 4) = %K[2](L*mu)})

The characteristic equation:

>	Determinant(M); eq := simplify(value(%)) = 0;

>	solve(eq, mu, allsolutions);

We conclude that the eigenvalues are .

A non-trivial solution of the system is in the null-space of :

>	eval(value(M), mu=n*Pi/L) assuming n::integer; N := NullSpace(%);

$Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 1, (2, 4) = 0, (3, 1) = (1/2)*cosh(n*Pi)+(1/2)*(-1)^n, (3, 2) = (1/2)*sinh(n*Pi), (3, 3) = (1/2)*cosh(n*Pi)-(1/2)*(-1)^n, (3, 4) = (1/2)*sinh(n*Pi), (4, 1) = (1/2)*cosh(n*Pi)-(1/2)*(-1)^n, (4, 2) = (1/2)*sinh(n*Pi), (4, 3) = (1/2)*cosh(n*Pi)+(1/2)*(-1)^n, (4, 4) = (1/2)*sinh(n*Pi)})$

{Vector[column](%id = 18446884696899531350)}

Here are the weights that go with the Krylov functions:

>	a_vals := convert([a[1,j] $j=1..4] =~ N[1], set);

and here is the deflection:

>	add(a[1,j]K[j](mux), j=1..4); eval(%, a_vals); # plug in the a_vals calculated above eval(%, mu=n*Pi/L); # assert that n is an integer

a[1, 1]*((1/2)*cosh(mu*x)+(1/2)*cos(mu*x))+a[1, 2]*((1/2)*sinh(mu*x)+(1/2)*sin(mu*x))+a[1, 3]*((1/2)*cosh(mu*x)-(1/2)*cos(mu*x))+a[1, 4]*((1/2)*sinh(mu*x)-(1/2)*sin(mu*x))

We see that the shape functions are simple sinusoids.

>

A single-span free-free beam

Here we calculate the natural modes of vibration of a single span

beam, free at both ends. The modes are of the form
.

The reasoning behind the calculations is very similar to that in the

previous section, therefore we don't comment on many details.

>	M := build_matrix(M_left_free, M_right_free(mu,L));

Matrix(4, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 1, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 1, (3, 1) = %K[3](L*mu), (3, 2) = %K[4](L*mu), (3, 3) = %K[1](L*mu), (3, 4) = %K[2](L*mu), (4, 1) = %K[2](L*mu), (4, 2) = %K[3](L*mu), (4, 3) = %K[4](L*mu), (4, 4) = %K[1](L*mu)})

The characteristic equation:

>	Determinant(M); simplify(value(%)) = 0; eq_tmp := isolate(%, cos(L*mu));

Let . Then the characteristic equation takes the form

>	eq := algsubs(L*mu=lambda, eq_tmp);

Here are the graphs of the two sides of the characteristic equation:

>	plot([lhs,rhs](eq), lambda=0..4*Pi, color=["red","Green"]);

The first three roots are:

>	lambda__1, lambda__2, lambda__3 := fsolve(eq, lambda=Pi/2..4*Pi, maxsols=3);

>	params := { L=1 };

>	mu__1, mu__2, mu__3 := (lambda__1, lambda__2, lambda__3) /~ eval(L, params);

>	plots:-display([ plot_beam(M, mu__1, params), plot_beam(M, mu__2, params), plot_beam(M, mu__3, params)], color=["red","Green","blue"], legend=[mode1, mode2, mode3]);

>

A single-span clamped-free cantilever

We have a cantilever beam of length . It is clamped at the

left end, and free at the right end.

>	M := build_matrix(M_left_clamped, M_right_free(mu,L));

Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, (3, 1) = %K[3](L*mu), (3, 2) = %K[4](L*mu), (3, 3) = %K[1](L*mu), (3, 4) = %K[2](L*mu), (4, 1) = %K[2](L*mu), (4, 2) = %K[3](L*mu), (4, 3) = %K[4](L*mu), (4, 4) = %K[1](L*mu)})

>	Determinant(M); simplify(value(%)) = 0; eq_tmp := isolate(%, cos(L*mu));

Let . Then the characteristic equation takes the form

>	eq := algsubs(L*mu=lambda, eq_tmp);

Here are the graphs of the two sides of the characteristic equation:

>	plot([lhs,rhs](eq), lambda=0..3*Pi, color=["red","Green"]);

>	lambda__1, lambda__2, lambda__3 := fsolve(eq, lambda=Pi/2..3*Pi, maxsols=3);

>	params := { L=1 };

>	mu__1, mu__2, mu__3 := (lambda__1, lambda__2, lambda__3) /~ eval(L, params);

>	plots:-display([ plot_beam(M, mu__1, params), plot_beam(M, mu__2, params), plot_beam(M, mu__3, params)], color=["red","Green","blue"], legend=[mode1, mode2, mode3]);

>

A dual-span pinned-pinned-free cantilever beam

We have a two-span beam of span lengths and , with the left end of the
first span pinned, the right end of the second span free, and the interface

between the spans on a pinned support. .

>	M := build_matrix( M_left_pinned, M_right_free(mu,L[2]), [ M_interface(mu,L[1])] );

The characteristic equation:

>	Determinant(M); eq_tmp1 := simplify(value(%)) = 0;

(1/4)*(-cos(L[1]*mu)*sinh(L[2]*mu)*cos(L[2]*mu)+cos(L[1]*mu)*sin(L[2]*mu)*cosh(L[2]*mu)+2*sin(L[1]*mu)*cosh(L[2]*mu)*cos(L[2]*mu)+2*sin(L[1]*mu))*sinh(L[1]*mu)+(1/4)*sin(L[1]*mu)*cosh(L[1]*mu)*(sinh(L[2]*mu)*cos(L[2]*mu)-sin(L[2]*mu)*cosh(L[2]*mu)) = 0

That equation does not seem to be amenable to simplification. The special case of , however,

is much nicer:

>	eval(eq_tmp1, {L[1]=L, L[2]=L}): eq_tmp2 := simplify(%*4);

Let :

>	eq_tmp3 := algsubs(L*mu=lambda, eq_tmp2);

That expression grows like , so we divide through by that to obtain

a better-behaved equation

>	eq := eq_tmp3/cosh(lambda)^2;

((4*cosh(lambda)*cos(lambda)+2)*sinh(lambda)*sin(lambda)+cos(lambda)^2-cosh(lambda)^2)/cosh(lambda)^2 = 0

>	plot(lhs(eq), lambda=0..2*Pi);

Here are the first three roots:

>	lambda__1, lambda__2, lambda__3 := fsolve(eq, lambda=1e-3..2*Pi, maxsols=3);

>	params := { L[1]=1, L[2]=1 };

>	mu__1, mu__2, mu__3 := (lambda__1, lambda__2, lambda__3) /~ eval(L[1], params);

>	plots:-display([ plot_beam(M, mu__1, params), plot_beam(M, mu__2, params), plot_beam(M, mu__3, params)], color=["red","Green","blue"], legend=[mode1, mode2, mode3]);

>

A dual-span clamped-pinned-free cantilever beam

We have a two-span beam of span lengths and , with the left end of the
first span clamped, the right end of the second span free, and the interface

between the spans on a pinned support. This is different from the previous

case only in the nature of the left boundary condition.

>	M := build_matrix( M_left_clamped, M_right_free(mu,L[2]), [ M_interface(mu,L[1])] );

The characteristic equation:

>	Determinant(M); eq_tmp1 := simplify(value(%)) = 0;

(1/4)*((-cos(L[1]*mu)*sin(L[2]*mu)-sin(L[1]*mu)*cos(L[2]*mu))*cosh(L[2]*mu)+cos(L[1]*mu)*sinh(L[2]*mu)*cos(L[2]*mu)-sin(L[1]*mu))*cosh(L[1]*mu)+(1/4)*(sinh(L[1]*mu)*cos(L[1]*mu)*cos(L[2]*mu)+sin(L[2]*mu))*cosh(L[2]*mu)+(1/4)*sinh(L[1]*mu)*cos(L[1]*mu)-(1/4)*sinh(L[2]*mu)*cos(L[2]*mu) = 0

That equation does not seem to be amenable to simplification. The special case of , however,

is much nicer:

>	eval(eq_tmp1, {L[1]=L, L[2]=L}): eq_tmp2 := simplify(%*4);

Let :

>	eq_tmp3 := algsubs(L*mu=lambda, eq_tmp2);

That expression grows like , so we divide through by that to obtain

a better-behaved equation

>	eq := eq_tmp3/cosh(lambda)^2;

>	plot(lhs(eq), lambda=0..2*Pi);

Here are the first three roots:

>	lambda__1, lambda__2, lambda__3 := fsolve(eq, lambda=1e-3..2*Pi, maxsols=3);

>	params := { L[1]=1, L[2]=1 };

>	mu__1, mu__2, mu__3 := (lambda__1, lambda__2, lambda__3) /~ eval(L[1], params);

>	plots:-display([ plot_beam(M, mu__1, params), plot_beam(M, mu__2, params), plot_beam(M, mu__3, params)], color=["red","Green","blue"], legend=[mode1, mode2, mode3]);

>

A triple-span free-pinned-pinned-free beam

We have a triple-span beam with span lengths of . The beam is supported

on two internal pinned supports. The extreme ends of the beam are free.

The graphs of the first three modes agree with those

in Figure 3.22 on page 70 of the 2007 article of
Henrik Åkesson, Tatiana Smirnova, Thomas Lagö, and Lars Håkansson.

In the caption of Figure 2.12 on page 28 the span lengths are given
as 3.5, 5.0, 21.5.

>	interface(rtablesize=12):

>	M := build_matrix( M_left_free, M_right_free(mu,L[3]), [ M_interface(mu,L[1]), M_interface(mu,L[2])] );

>	params := { L[1]=3.5, L[2]=5.0, L[3]=21.5 };

The characteristic equation

>	simplify(Determinant(M)): value(%): eq := simplify(eval(%, params));

>	plot(eq, mu=0..0.4);

That graphs grows much too fast to be useful. We moderate it by dividing through
the fastest growing cosh term:

>	plot(eq/cosh(21.5*mu), mu=0..0.4);

Here are the first three roots:

>	mu__1, mu__2, mu__3 := fsolve(eq, mu=1e-3..0.4, maxsols=3);

>	plots:-display([ plot_beam(M, mu__1, params), plot_beam(M, mu__2, params), plot_beam(M, mu__3, params)], color=["red","Green","blue"], legend=[mode1, mode2, mode3]);

>

A triple-span free-spring-spring-free beam

We have a triple-span beam with span lengths of . The beam is supported

on two internal springy supports. The extreme ends of the beam are free.
The numerical data is from the worksheet posted on July 29, 2020 at
https://www.mapleprimes.com/questions/230085-Elasticfoundation-Multispan-EulerBernoulli-Beamthreespan#comment271586

The problem is pretty much the same as the one in the previous section, but the

pinned supports have been replaced by spring supports.

This section's calculations require a little more precision than

Maple's default of 10 digits:

>	Digits := 15;

>	interface(rtablesize=12):

>	M := build_matrix(M_left_free, M_right_free(mu,L[3]), [ M_interface(mu, L[1], 'Kd'=kd/(EI), 'Kt'=kt/(EI)), M_interface(mu, L[2], 'Kd'=kd/(EI), 'Kt'=kt/(EI)) ]);

$Matrix(12, 12, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 1, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (1, 10) = 0, (1, 11) = 0, (1, 12) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 1, (2, 5) = 0, (2, 6) = 0, (2, 7) = 0, (2, 8) = 0, (2, 9) = 0, (2, 10) = 0, (2, 11) = 0, (2, 12) = 0, (3, 1) = %K[1](L[1]*mu), (3, 2) = %K[2](L[1]*mu), (3, 3) = %K[3](L[1]*mu), (3, 4) = %K[4](L[1]*mu), (3, 5) = -1, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (3, 9) = 0, (3, 10) = 0, (3, 11) = 0, (3, 12) = 0, (4, 1) = %K[4](L[1]*mu), (4, 2) = %K[1](L[1]*mu), (4, 3) = %K[2](L[1]*mu), (4, 4) = %K[3](L[1]*mu), (4, 5) = 0, (4, 6) = -1, (4, 7) = 0, (4, 8) = 0, (4, 9) = 0, (4, 10) = 0, (4, 11) = 0, (4, 12) = 0, (5, 1) = %K[3](L[1]*mu), (5, 2) = %K[4](L[1]*mu), (5, 3) = %K[1](L[1]*mu), (5, 4) = %K[2](L[1]*mu), (5, 5) = 0, (5, 6) = -I*kt/(E*mu), (5, 7) = -1, (5, 8) = 0, (5, 9) = 0, (5, 10) = 0, (5, 11) = 0, (5, 12) = 0, (6, 1) = %K[2](L[1]*mu), (6, 2) = %K[3](L[1]*mu), (6, 3) = %K[4](L[1]*mu), (6, 4) = %K[1](L[1]*mu), (6, 5) = -I*kd/(E*mu^3), (6, 6) = 0, (6, 7) = 0, (6, 8) = -1, (6, 9) = 0, (6, 10) = 0, (6, 11) = 0, (6, 12) = 0, (7, 1) = 0, (7, 2) = 0, (7, 3) = 0, (7, 4) = 0, (7, 5) = %K[1](L[2]*mu), (7, 6) = %K[2](L[2]*mu), (7, 7) = %K[3](L[2]*mu), (7, 8) = %K[4](L[2]*mu), (7, 9) = -1, (7, 10) = 0, (7, 11) = 0, (7, 12) = 0, (8, 1) = 0, (8, 2) = 0, (8, 3) = 0, (8, 4) = 0, (8, 5) = %K[4](L[2]*mu), (8, 6) = %K[1](L[2]*mu), (8, 7) = %K[2](L[2]*mu), (8, 8) = %K[3](L[2]*mu), (8, 9) = 0, (8, 10) = -1, (8, 11) = 0, (8, 12) = 0, (9, 1) = 0, (9, 2) = 0, (9, 3) = 0, (9, 4) = 0, (9, 5) = %K[3](L[2]*mu), (9, 6) = %K[4](L[2]*mu), (9, 7) = %K[1](L[2]*mu), (9, 8) = %K[2](L[2]*mu), (9, 9) = 0, (9, 10) = -I*kt/(E*mu), (9, 11) = -1, (9, 12) = 0, (10, 1) = 0, (10, 2) = 0, (10, 3) = 0, (10, 4) = 0, (10, 5) = %K[2](L[2]*mu), (10, 6) = %K[3](L[2]*mu), (10, 7) = %K[4](L[2]*mu), (10, 8) = %K[1](L[2]*mu), (10, 9) = -I*kd/(E*mu^3), (10, 10) = 0, (10, 11) = 0, (10, 12) = -1, (11, 1) = 0, (11, 2) = 0, (11, 3) = 0, (11, 4) = 0, (11, 5) = 0, (11, 6) = 0, (11, 7) = 0, (11, 8) = 0, (11, 9) = %K[3](L[3]*mu), (11, 10) = %K[4](L[3]*mu), (11, 11) = %K[1](L[3]*mu), (11, 12) = %K[2](L[3]*mu), (12, 1) = 0, (12, 2) = 0, (12, 3) = 0, (12, 4) = 0, (12, 5) = 0, (12, 6) = 0, (12, 7) = 0, (12, 8) = 0, (12, 9) = %K[2](L[3]*mu), (12, 10) = %K[3](L[3]*mu), (12, 11) = %K[4](L[3]*mu), (12, 12) = %K[1](L[3]*mu)})$

Calculate the determinant of . The result is quite large, so we terminate the command
with a colon so that not to have to look at the result. If we bothered to peek, however, we
will see that the determinant has a factor of . But that quite obvious by looking at the
entries of the matrix shown above. Two of its rows have in them and another two have
. When multiplied, they produce the overall factor of .

>	DET := Determinant(M):

Here are the parameters that the determinant depends on:

>	indets(DET, name); # the parameters that make up M

So we provide values for those parameters:

>	params := { L[1]=3.5, L[2]=5.0, L[3]=21.5, kd=4.881e9, kt=1.422e4, E = 2.05e11, I = 1.1385e-7 };

Here is the characteristic equation. We multiply it by to remove the singularity at

>	mu^8 * value(DET): eq := eval(%, params):

>	plot(eq, mu=0..0.6);

We can't see anything useful in that graph. Let's limit the vertical range:

>	plot(eq, mu=0..0.6, view=-1e8..1e8);

>	mu__1, mu__2, mu__3 := fsolve(eq, mu=1e-3..0.6, maxsols=3);

>	plots:-display([ plot_beam(M, mu__1, params), plot_beam(M, mu__2, params), plot_beam(M, mu__3, params)], color=["red","Green","blue"], legend=[mode1, mode2, mode3]);

>	Digits := 10; # restore the default

>

Download krylov-duncan.mw

Introducing Maple Learn

by: Karishma 885 Product: Maple Learn

August 18 2020

9 0

I am very pleased to announce that we have just begun a free public beta for a new online product, Maple Learn! Maple Learn is a dynamic online environment designed specifically for teaching and learning math and solving math problems, from mid-high school to second year university.

Maple Learn is much more than just a sophisticated online graphing calculator. We tried to create an environment that focuses on the things instructors and students in those courses have told us that they want/need in a math tool. Here are some of my personal favorites:

You can get the answer directly if you want it, but you can also work out problems line-by-line as you would on paper, or use a combination of manual steps and computations performed by Maple Learn.
The plot of your expression shows up as soon as you start typing, so plotting is super easy
You can parameterize expressions with a single mouse click, and then watch the plots and results change as you modify the values using sliders
It’s really easy to share your work, so when a student asks for help, the helper can always see exactly what they’ve done so far (and it will be legible, unlike a lot of the tutoring I’ve done!)

Free public beta: Maple Learn is freely available to instructors and students as part of an on-going public beta program. Please try it out, and feel free to use with your classes this fall.

Visit Maple Learn for more information, and to try it out. We hope you find it useful, and we’d love to know what you think.

Maple's user-interface deters me from doing math

by: DoingMath2018 35

August 10 2020

1 17

I don't stand by all I said here, so I'm deleting. I think there was valid criticism to be made, but I didn't make it any valid way.

Sampling linearly correlated random variables

by: mmcdara 7896 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 14092 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 1370 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 7896 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 1370 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

Creating Help - using HelpTools

by: dharr 8917 Product: Maple

July 19 2020

5 3

@ianmccr posted here about making help for a package using makehelp. Here I show how to do this with the HelpTools package.

The attached worksheet shows how to create the help database for the Orbitals package available at the Applications Center or in the Maple Cloud. The help pages were created as worksheets - start using an existing help page as a model - use View/Open Page As Worksheet and then save from there. The topics and other information are entered by adding Attributes under File/Document Properties - for example for the realY help page these are:

Active=false means a regular help page; Active=true means an example worksheet.

There may be several aliases, for example the cartesion help page also describes the fullcartesian command and so Alias is: Orbitals[fullcartesian],cartesian,fullcartesian and Keyword is: Orbitals,cartesian,fullcartesian

Once a worksheet is created for each help page they are assembled into the help database with the attached file. More information is in the attached file

Orbitals_Make_help_database.mw

Creating Help - An Exploration of makehelp

by: ianmccr 170 Product: Maple 2019

July 18 2020

3 0

MacDude posted a worksheet for creating help files using the makehelp command that I found very useful. However, his worksheet created a table of contents that organized the command help pages under a single folder. I wanted to create a help file table of contents with the commands organized into sub-folders under the main application folder. ie. Package Name,Folder Name, command. The help page for makehelp is not very informative; in particular the example that shows (purportedly) how to override the existing help pages with your own help file is very misleading. Also, the description of the parameters to the browser option of the makehelp command is too vague. I needed an error message to tell me that the browser option expects parameters in the form List(Name,String). In the end, I was able to get a folder/sub-folder structure using a structure [`Folder Name`,`Subfolder Name`,"Command"]. Sub-folders are sorted alphabetically. If anyone has a need, the attached worksheet shows how I created the table of contents structure.

helpcode.mw

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

by: Samir Khan 2181 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 7896 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 7896

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

Maple Animation of 150 Days of Corona in the US ...

by: MathField 40 Product: Maple 2020

June 23 2020

7 2

This is an animation of the spread of the COVID-19 over the U.S. in the first 150 days. It was created in Maple 2020, making extensive use of DataFrames.

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

The animation of 150 Day history includes COVID-19 data published by the NY Times and geographic data assembled from other sources. Each cylinder represents a county or in two special cases New York City and Kansas City. The cross-sectional area of each cylinder is the area in square miles of the corresponding county. The height of each cylinder is on a logarithmic scale (in particular it is 100*log base 2 of the case count for the county. The argument of the logarithm function is the number of cases per county divided by the are in square miles-so an areal density. Using a logarithmic scale facilitates showing super high density areas (e.g., NYC) along with lower density areas. The heights are scaled by a prefactor of 100 for visualization.

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!

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

MaplePrimes Posts

MaplePrimes Posts are for sharing your experiences, techniques and opinions about Maple, MapleSim and related products, as well as general interests in math and computing.

Latest Post
Latest Posts Feed

Eigenfunctions of a multi-span Euler-Bernoulli...

Introducing Maple Learn

Maple's user-interface deters me from doing math

Sampling linearly correlated random variables

International Mathematics Competition for University...

Circumscribed ellipse

Doing plots "a la XKCD"

An ellipse inscribed in a given triangle

Creating Help - using HelpTools

Creating Help - An Exploration of makehelp

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

Realistic 3D rendering of a complex system

Numerical evaluation of integrals over polygons...

Maple Animation of 150 Days of Corona in the US ...

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

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again

Ask a Question

Create a Post

MaplePrimes Posts

MaplePrimes Posts are for sharing your experiences, techniques and opinions about Maple, MapleSim and related products, as well as general interests in math and computing.

Latest Post Latest Posts Feed

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Latest Post
Latest Posts Feed