Featured Post

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


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:

  1. Each sick cat infects two other cats.
  2. Each sick cat infects three other cats.
  3. 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!

Featured Post

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



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

"rho ((∂)^2u)/((∂)^( )t^2)+ ((∂)^2)/((∂)^( )x^2)(EI ((∂)^(2)u)/((∂)^( )x^(2)))=f"

subject to initial and boundary conditions.  The solution u = u(x, t) is the

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

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

the beam's material, and I 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 rho, E, I are constants for simplicity. Since only
the product of the coefficients E and I enters the calculations, we lump the two

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

Commonly, EI 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 u, u__x, u__xx, u__xxx 
(that's 36 possible combinations!) where I have used subscripts to indicate

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

set of boundary conditions:
u(0, t) = 0, u__xx(0, t) = 0, u__xx(L, t) = 0, u__xxx(t) = sin*t.   (Oops, coorection, that last
condition was meant to be uxxx(L,t) = sin t.)

Additionally, the PDE needs to be supplied with initial conditions which express

the initial displacement and the initial velocity:
"u(x,0)=phi(x),   `u__t`(x,0)=psi(x)."


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(L, n, options)



        L:  List of span lengths, in order from left to right, as in [L__1, L__2 .. (), `L__ν`].

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




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


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 L may be entered as a number
rather than as a list. That is, L__1 is equivalent to [L__1].



        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 0 < t and t < T (default: T = 1)

        F: an expression in x and t that describes the applied force f(x, t)  (default: F = 0)
        IC: the list [u(x, 0), u__t(x, 0)]of the initial displacement and velocity,  as
                expressions in x (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
                u = alpha(t), u_x = beta(t), u_xx = gamma(t), u_xxx = delta(t). The right-hand sides of these equations

                can be any expression in t.  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 = [u, u_xx] which

                corresponds to a simple support).

        BC_right: like BC_left, but for the right end of the overall beam (default: BC_right = "[u,u_xx])"


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 u(x, t) at time t.  If the discretization in the x 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 u(x, t).  The optional m = val specification requests

                a grid consisting of val subintervals in the time direction (default: "m=25)"

                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.


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

        is considered to be horizontal, then u is the vertical displacement, `u__x ` is the slope,  EI*u__xx

        is the bending moment, and EI*u__xxx is the transverse shear force.


A single-span simply-supported beam with initial velocity


The function u(x, t) = sin(Pi*x)*sin(Pi^2*t) is an exact solution of a simply supported beam with

"u(x,0)=0,   `u__t`(x,0)=Pi^(2)sin(Pi x)."  The solution is periodic in time with period 2/Pi.

sol := beamsolve(1, 25, 'T'=2/Pi, 'IC'=[0, Pi^2*sin(Pi*x)]):

The initial condition u(x, 0) = 0, u__t(x, 0) = 1  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 "u(x,0)=g(x),  `u__t`(x,0)=0," where g(x) is the
first eigenmode of its free vibration (calculated in another spreadsheet).  The motion is
periodic in time, with period "1.787018777."

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

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.1*sin(10*t),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):
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):
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]);



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]):
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


Get eigenvectors in Maple

Maple asked by WYing 5 June 29