tomleslie

12497 Reputation

19 Badges

12 years, 231 days

MaplePrimes Activity


These are answers submitted by tomleslie

you can get a solution as shown in the attached.

alpha[1](t), alpha[2](t), m[1](t) and m[2](t) all appear to be complex whereas z[1](t) and z[2](t) appear to be real.

(NB the plots are much nicer  in the worksheet than they render on this site - honest)
 

  restart;
  eqns := { diff(alpha[1](t), t) = -w[r]*alpha[1](t)*I - g*m[1](t)*I + J*alpha[2](t)*I,
            diff(alpha[2](t), t) = -w[r]*alpha[2](t)*I - g*m[2](t)*I + J*alpha[1](t)*I,
            diff(m[1](t), t) = -2*I*w[a]*m[1](t) + g*alpha[1](t)*z[1](t)*I,
            diff(m[2](t), t) = -2*I*w[a]*m[2](t) + g*alpha[2](t)*z[2](t)*I,
            diff(z[1](t), t) = -2*I*g*(alpha[1](t)*conjugate(m[1](t)) - conjugate(alpha[1](t))*m[1](t)),
            diff(z[2](t), t) = -2*I*g*(alpha[2](t)*conjugate(m[2](t)) - conjugate(alpha[2](t))*m[2](t))
          };
  cons := {J = 0.1, g = 1, w[a] = 100, w[r] = 100};
  ics := {alpha[1](0) = sqrt(20), alpha[2](0) = 0, m[1](0) = 0, m[2](0) = 0, z[1](0) = -1, z[2](0) = -1};

{diff(alpha[1](t), t) = -I*w[r]*alpha[1](t)-I*g*m[1](t)+I*J*alpha[2](t), diff(alpha[2](t), t) = -I*w[r]*alpha[2](t)-I*g*m[2](t)+I*J*alpha[1](t), diff(m[1](t), t) = -(2*I)*w[a]*m[1](t)+I*g*alpha[1](t)*z[1](t), diff(m[2](t), t) = -(2*I)*w[a]*m[2](t)+I*g*alpha[2](t)*z[2](t), diff(z[1](t), t) = -(2*I)*g*(alpha[1](t)*conjugate(m[1](t))-conjugate(alpha[1](t))*m[1](t)), diff(z[2](t), t) = -(2*I)*g*(alpha[2](t)*conjugate(m[2](t))-conjugate(alpha[2](t))*m[2](t))}

 

{J = .1, g = 1, w[a] = 100, w[r] = 100}

 

{alpha[1](0) = 2*5^(1/2), alpha[2](0) = 0, m[1](0) = 0, m[2](0) = 0, z[1](0) = -1, z[2](0) = -1}

(1)

  dsol:=dsolve( eval( `union`(eqns, ics), cons), numeric);
  dsol(0.1);
#
# Hmmmm - lookks as if alpha[1](t), alpha[2](t), m[1](t), m[2](t)
# are all complex whereas z[1](t) and z[2](t) are real.
#
  plots:-odeplot( dsol,
                  [ [t, Re(alpha[1](t))],
                    [t, Im(alpha[1](t))]
                  ],
                  t=0..0.2,
                  color=[red, blue],
                  title=typeset(Re(alpha[1](t)), " and ", Im(alpha[1](t))),
                  titlefont=[times, italic, 16]
                );
  plots:-odeplot( dsol,
                  [ [t, Re(alpha[2](t))],
                    [t, Im(alpha[2](t))]
                  ],
                  t=0..0.2,
                  color=[red, blue],
                  title=typeset(Re(alpha[2](t)), " and ", Im(alpha[2](t))),
                  titlefont=[times, italic, 16]
                );
  plots:-odeplot( dsol,
                  [ [t, Re(m[1](t))],
                    [t, Im(m[1](t))]
                  ],
                  t=0..0.2,
                  color=[red, blue],
                  title=typeset(Re(m[1](t)), " and ", Im(m[1](t))),
                  titlefont=[times, italic, 16]
                );
  plots:-odeplot( dsol,
                  [ [t, Re(m[2](t))],
                    [t, Im(m[2](t))]
                  ],
                  t=0..0.2,
                  color=[red, blue],
                  title=typeset(Re(m[2](t)), " and ", Im(m[2](t))),
                  titlefont=[times, italic, 16]
                );
  plots:-odeplot( dsol,
                  [ [t, z[1](t)],
                    [t, z[2](t)]
                  ],
                  t=0..0.2,
                  color=[black, green],
                  title=typeset( z[1](t), " and ", z[2](t)),
                  titlefont=[times, italic, 16]
                );

`Non-fatal error while reading data from kernel.`

 

[t = .1, alpha[1](t) = -3.75412937362476+2.42843979680634*I, alpha[2](t) = -0.242878255083043e-1-0.375392393317978e-1*I, m[1](t) = 0.536167658802068e-1-0.653827773144382e-1*I, m[2](t) = 0.187160998874069e-3+0.440484190777856e-3*I, z[1](t) = -.985597053964343+0.*I, z[2](t) = -.999999541889001+0.*I]

 

 

 

 

 

 

 

Download odeComp.mw

I have no idea what you are trying to achieve. For example in the code you present

Cen := proc(M, N, R) local eq1, eq2, sol;
eq1 := EqBIS(M, N, R) = 0;
eq2 := EqBIS(N, M, R) = 0;
sol := simplify(solve({eq1, eq2}, {x, y}));
RETURN([subs(sol, x), subs(sol, y)]);
end proc
  Ce := Cen(P, Dd, Q);

what are P, Dd, Q?

Am I expected to guess???

If I do make some more-or-less random guesses - required because of your inability to ask a sensible question or post reasonably functional code using the big green up-arrow in the Mapleprimes toolbar, then I'd probably come up with what is shown in the attached

I assume that if this answer is *satisfactory* then you will delete the whole thread?

  restart:
  with(geometry):
  _EnvHorizontalName:= x:
  _EnvVerticalName:= y:
  triangle( ABC,
            [ point(A, [4, 5]),
              point(B, [11, 7/3]),
              point(C, [11, 5])
            ]
          ):
  bisector( bA, A, ABC ):
  incircle( inc, ABC, centername=incc):
  evalf(coordinates(incc));

  draw( [  A(symbol=solidcircle, symbolsize=20, color=red),
           B(symbol=solidcircle, symbolsize=20, color=red),
           C(symbol=solidcircle, symbolsize=20, color=red),
           bA(color=black),
           ABC(color=red),
           inc(color=blue),
           incc(symbol=solidcircle, symbolsize=20, color=blue)
        ],
        axes=boxed,
        view=[0..12, 0..6],
        scaling=constrained
      );

[9.912034177, 3.912034175]

 

 

 

Download useGeom2.mw

same calculation using the geometry() package. See the attached

  restart;
  EqBIS:= proc(P, U, V)
               uses LinearAlgebra:
               local a, M1, t;
               a := (P - U)/Norm(P - U, 2) + (P - V)/Norm(P - V, 2);
               M1 := P + a*t;
               return eliminate( { x = M1[1], y = M1[2]},
                                 t
                               )[2][]=0;
          end proc:
  evalf
  ( isolate
    ( EqBIS( Vector([4, 5]),
             Vector([11, 7/3]),
             Vector([11, 5])
           ),
      y
    )
  );

y = 5.736102529-.1840256318*x

(1)

  restart:
  with(geometry):
  _EnvHorizontalName := x:
  _EnvVerticalName := y:
  triangle( ABC, [ point(A, [4, 5]),
                   point(B, [11, 7/3]),
                   point(C, [11, 5])
                 ]
          ):
  evalf
  ( isolate
    ( Equation
      ( bisector( bA, A, ABC )
      ),
      y
    )
  );

y = -.1840256319*x+5.736102529

(2)

 

Download useGeom.mw

Unfortunately your system of equations is inconsistent, which is demonstrated in a couple of ways in the attached

restart

Eq[0, 0] := 1.33120000000000000000000000000*lambda[0, 1]+1.33120000000000000000000000000*lambda[0, 2]+1.33120000000000000000000000000*lambda[0, 3]+1.33120000000000000000000000000*lambda[0, 4] = .916487142969312002551492271668

Eq[0, 1] := 1.32901933598375616624661615670*lambda[0, 1]+1.32901933598375616624661615670*lambda[0, 2]+1.32901933598375616624661615670*lambda[0, 3]+1.32901933598375616624661615670*lambda[0, 4] = 1.09232395220587507357427365904

Eq[0, 2] := 1.37120000000000000000000000000*lambda[0, 1]+1.37120000000000000000000000000*lambda[0, 2]+1.37120000000000000000000000000*lambda[0, 3]+1.37120000000000000000000000000*lambda[0, 4] = 1.25415129307905065856083635281

Eq[0, 3] := .966980664016243833753383843299*lambda[0, 1]+.966980664016243833753383843299*lambda[0, 2]+.966980664016243833753383843299*lambda[0, 3]+.966980664016243833753383843299*lambda[0, 4] = 1.37114174964252179339832329224

Eq[0, 0]/(1.33120000000000000000000000000); Eq[0, 1]/(1.32901933598375616624661615670); Eq[0, 2]/(1.37120000000000000000000000000); Eq[0, 3]/(.966980664016243833753383843299)

1.000000000*lambda[0, 1]+1.000000000*lambda[0, 2]+1.000000000*lambda[0, 3]+1.000000000*lambda[0, 4] = .6884669043

 

1.000000000*lambda[0, 1]+1.000000000*lambda[0, 2]+1.000000000*lambda[0, 3]+1.000000000*lambda[0, 4] = .8219022270

 

1.000000000*lambda[0, 1]+1.000000000*lambda[0, 2]+1.000000000*lambda[0, 3]+1.000000000*lambda[0, 4] = .9146377574

 

1.000000000*lambda[0, 1]+1.000000000*lambda[0, 2]+1.000000000*lambda[0, 3]+1.000000000*lambda[0, 4] = 1.417961911

(1)

with(LinearAlgebra); A, b := GenerateMatrix([Eq[0, 0], Eq[0, 1], Eq[0, 2], Eq[0, 3]], indets(Eq[0, 0])); LinearSolve(A, b)

Matrix(%id = 36893488148293117220), Vector[column](%id = 36893488148293117100)

 

Error, (in LinearAlgebra:-BackwardSubstitute) inconsistent system

 

``

Download inconSys.mw

precisely what you want (or why!), Rather obviously, if

xi=xi-1+b

then the answer will be

{x__0, x__0 + b, x__0 + 2*b, x__0 + 3*b, x__0 + 4*b, x__0 + 5*b, x__0 + 6*b, x__0 + 7*b, x__0 + 8*b, x__0 + 9*b, x__0 + 10*b}

While you can generate this set recursively (as shown in the attached) it really seems like doing it the *hard way*!

  restart:
  n:=10:
  rProc:= proc(N, xv:=x__0)
               if   N=0
               then return xv;
               else return xv, rProc(N-1, xv+b)
               fi;
          end proc:

{rProc(n)};

{x__0, x__0+b, x__0+2*b, x__0+3*b, x__0+4*b, x__0+5*b, x__0+6*b, x__0+7*b, x__0+8*b, x__0+9*b, x__0+10*b}

(1)

 

Download recur.mw

in what context you want to do this - it just *seems* a bit weird - but the attached will do it

  restart:

  nArgs:= expr-> numelems([op([1..-1], expr)]);

  nArgs('f'(x1, x2, x3));
  nArgs('g'(x1, x2, x3, x4, x5, x6, x7));

proc (expr) options operator, arrow; numelems([op([1 .. -1], expr)]) end proc

 

3

 

7

(1)

 

Download args.mw

what you are trying to achieve here.- so the attached just contains some hints and tips which you might find useful - or maybe not!

There are many, ways in Maple to  perform multiple assignments - the attached shows a couple. In general I don't think it is a good idea to ask how to replicate the behaviour of a command in some other language. Much better to specify exactly what you want Maple to do.

The *simplest* way of getting "user input" is probably via the readline() command. Provided that you realise that anything entered in response to a  readline() command will be a *string*, and therefore you may have to "convert" appropriately, this approach can be surprisingly flexible. The attached shows how to do a simple variable assignment based on user input and also a simple guessing game, which I have shamelessly replicated from the  "help" for the readline() command.

So if these do not satisfy your requirement(s), then ithink yu need to specify your requirements more carefully without reference to other software packages.

  restart:
#
# Assigning multiple matrixes of the same size
#
  with(LinearAlgebra):
  randomize():
  n:=3:
  m:=5:
  A, B, C:= seq( RandomMatrix(n, m), j=1..3):
  A;
  B;
  C;

Matrix(3, 5, {(1, 1) = 2, (1, 2) = -22, (1, 3) = 47, (1, 4) = -44, (1, 5) = 31, (2, 1) = -37, (2, 2) = -1, (2, 3) = -96, (2, 4) = 21, (2, 5) = -90, (3, 1) = 6, (3, 2) = 90, (3, 3) = 27, (3, 4) = 45, (3, 5) = -46})

 

Matrix(3, 5, {(1, 1) = 91, (1, 2) = 74, (1, 3) = -55, (1, 4) = -91, (1, 5) = -52, (2, 1) = 65, (2, 2) = -46, (2, 3) = 28, (2, 4) = -2, (2, 5) = -3, (3, 1) = -92, (3, 2) = -17, (3, 3) = -70, (3, 4) = 96, (3, 5) = 68})

 

Matrix(%id = 36893488147886702828)

(1)

  restart:
#
# Initialize one matrix, then define several of the
# same size
#
  with(LinearAlgebra):
  M:=Matrix(3, 5);
  A, B, C:= seq( Matrix( op(1, M)), j=1..3):
  A;
  B;
  C;
  

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

 

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

 

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

 

Matrix(%id = 36893488147952599156)

(2)

  restart:
#
# Really simple user input
#
 printf( "enter a value for the variable x\n"):
 x:= parse(readline()):
 x;

enter a value for the variable x

 

3.13

(3)

  restart;
#
# Using readline() to get user input
#
  GetNumber:= proc()
                   local number;
                   randomize();
                   number := rand( 1 .. 5 )();
                   printf( "Guess a number between 1 and 5:\n>> " );
                   while( parse( readline() ) <> number ) do
                          printf( "Try again!\n>> " );
                   end do:
                   printf( "Correct! The number was %d.\n", number );
              end proc:
   GetNumber();

Guess a number between 1 and 5:
>>

 

Try again!
>>

 

Try again!
>>

 

Correct! The number was 4.

 

 

 

Download tips.mw

 

 

without using the geometry() package.

  restart;
  with(plots):
  a := 11:
  b := 5:
  ell := (x*cos(theta) + y*sin(theta))^2/a^2 + (x*sin(theta) - y*cos(theta))^2/b^2 = 1:

  pltAll:=  proc(ang)
                 local p1,d;
                 uses plots, plottools:
                 p1:= implicitplot( eval(ell, theta=ang),
                                    x = -12 .. 12,
                                    y = -12 .. 12,
                                    thickness = 3,
                                    color = blue
                                  ):
                 d:= getdata(p1)[2]:
                 display( [ p1,
                            plot( [ [op(1, d[1]), op(1, d[2])],
                                    [op(1, d[1]), op(2, d[2])],
                                    [op(2, d[1]), op(2, d[2])],
                                    [op(2, d[1]), op(1, d[2])],
                                    [op(1, d[1]), op(1, d[2])]
                                  ],
                                  color=red
                                )
                          ]
                        );
            end proc:
  animate( display,
           ['pltAll'(theta)],
           theta=0..2*Pi
         );

 

 

Download sameAnim.mw

since this poster's questions have a habit of disappearing after they have been answered, I can come up with two possible animations which might be desired. Both are shown in the attached,

I'd guess the OP probably wants the first one (but my guesses are often wrong!)

  restart:
  with(plots):
  with(geometry):
  _EnvHorizontalName := 'x':
  _EnvVerticalName := 'y':

  point(OO, [0,0]):
  a := 11:
  b := 5:
  ellipse( ell, x^2/a^2+y^2/b^2=1):

  aProc:= proc(ang)
               local ellr,plt,d;
               uses geometry:
               rotation(ellr, ell, ang, counterclockwise, OO):
               plt:=draw(ellr):
               d:=plottools:-getdata(plt)[2];
               point(pt1,[op(1, d[1]), op(1, d[2])]):
               point(pt2,[op(1, d[1]), op(2, d[2])]):
               point(pt3,[op(2, d[1]), op(2, d[2])]):
               point(pt4,[op(2, d[1]), op(1, d[2])]):
               segment(s1, [pt1, pt2]):
               segment(s2, [pt2, pt3]):
               segment(s3, [pt3, pt4]):
               segment(s4, [pt4, pt1]):
               draw( [ ellr(color=blue, thickness=3),
                       s1(color=red),
                       s2(color=red),
                       s3(color=red),
                       s4(color=red)
                     ],
                      axes=none
                    ):    
          end proc:

  bProc:= proc(ang)
               uses geometry:
               rotation(ellr, ell, ang, counterclockwise, OO):
               point(p1, [ -a, -b]):
               point(p2, [ -a,  b]):
               point(p3, [  a,  b]):
               point(p4, [  a, -b]):
               segment( s1, [p1, p2]):
               segment( s2, [p2, p3]):
               segment( s3, [p3, p4]):
               segment( s4, [p4, p1]):
               rotation( s1r, s1, ang, counterclockwise, OO ):
               rotation( s2r, s2, ang, counterclockwise, OO ):
               rotation( s3r, s3, ang, counterclockwise, OO ):
               rotation( s4r, s4, ang, counterclockwise, OO ):
               draw( [ ellr(color=blue, thickness=3),
                       s1r(color=red),
                       s2r(color=red),
                       s3r(color=red),
                       s4r(color=red)
                     ],
                     axes=none
                   ):
         end proc:

  animate( display,
           ['aProc'(theta)],
           theta=0..2*Pi
         );
  animate( display,
           [  ['bProc'(theta)]],
           theta=0..2*Pi
         );

               

 

 

 

Download anims.mw

 

your pde is second order in time and second order in space so a total of four boundary/initial conditions would be needed for a complete solution. One could (in principle) get an analytic solution containing an arbitrary constant. But since Maple cannot find an analytic solution, then the only option is to solve numerically. For a numeric solution four boundary/initial condition is mandatory - so I made one up (higlighted in blue in the attached).

Now a numeric solution is possible - see the attached, which for some reason is refusing to display inline on this site!

Download pdSolve.mw

that there is a "feasible region" whihc satisfies all four inequalities - see the attached

  ineqs:=[  alpha <= 0.0002500000000*(-18000.*m^2 + 47271.*m + 39514. + sqrt(3.24000000*10^8*m^4 - 1.701756000*10^9*m^3 - 4.266980559*10^9*m^2 - 3.036299412*10^9*m - 6.95987804*10^8))/(9.*m^2 + 12.*m + 4.),
  0.00005000000000*(-90000.*m^2 + 237237.*m + 198158. + sqrt(8.100000000*10^9*m^4 - 4.270266000*10^10*m^3 - 1.069976858*10^11*m^2 - 7.612670111*10^10*m - 1.744924704*10^10))/(9.*m^2 + 12.*m + 4.) <= alpha,
  -0.6666666667 < m,
   m < -0.6665522013
  ];

[alpha <= 0.2500000000e-3*(-18000.*m^2+47271.*m+39514.+(324000000.0*m^4-1701756000.*m^3-4266980559.*m^2-3036299412.*m-695987804.0)^(1/2))/(9.*m^2+12.*m+4.), 0.5000000000e-4*(-90000.*m^2+237237.*m+198158.+(8100000000.*m^4-0.4270266000e11*m^3-0.1069976858e12*m^2-0.7612670111e11*m-0.1744924704e11)^(1/2))/(9.*m^2+12.*m+4.) <= alpha, -.6666666667 < m, m < -.6665522013]

(1)

  p1:=plots:-inequal( ineqs[1], m=lhs(ineqs[3])..rhs(ineqs[4]), alpha=15000..200000, color=red):
  p2:=plots:-inequal( ineqs[2], m=lhs(ineqs[3])..rhs(ineqs[4]), alpha=15000..200000, color=blue):
  plots:-display([p1,p2]);

 

 

Download ineqs.mw

but maybe somethig in the attached will help?

  restart:
  exact:=u(x,t)=sin(x)*cos(t);
  approx:=u(x,t)=sin(x)*sin(t);
  plot3d( [rhs(exact), rhs(approx)], x=0..2*Pi, t=0..3*Pi, color=[red, blue], style=surface);

u(x, t) = sin(x)*cos(t)

 

u(x, t) = sin(x)*sin(t)

 

 

 

Download plt3.mw

 

if I have to do it the hard way (why?), then it would probably look something like the attached

  restart:
  with(plottools):
  with(plots):
  _local(D):

  a:= 7:
  b:= a/2:
  r:= b/2:
  c:= 2*sqrt(10)*r:

  Cir1:= disk([c/2 - r, -c/2 + r], r, color = blue):
  Cir2:= disk([-c/2 + r, -c/2 + r], r, color = blue):
  Cir3:= disk([-c/2 + r, c/2 - r], r, color = blue):
  Cir4:= disk([c/2 - r, c/2 - r], r, color = blue):

  rot:= Matrix(2, 2, [[cos(theta), -sin(theta)], [sin(theta), cos(theta)]]):
  T:= rot.(Vector(2, [x, y])):
  Ell:= x^2/a^2 + y^2/b^2 = 1:
  doRot:= ang -> expand(eval(Ell, eval([x = T[1], y = T[2]], theta = ang))):
  e1:= doRot(Pi/4):
  e2:= doRot(-Pi/4):

  p1:= shadebetween( rhs(solve(e2,[y])[2][]), -c/2, x=-c/2..0, scaling=constrained):
  p2:= display([p1, reflect(p1, [[ 0, -c/2], [0, c/2]])]):
  p3:= display([p2, reflect(p2, [[ -c/2, 0], [c/2, 0]])]):
  p4:= display([p3, rotate(p3, Pi/2, [0,0])], color=red):

  Ell1:= implicitplot(e1, x = -6 .. 6, y = -6 .. 6, thickness = 3, color = blue):
  Ell2:= implicitplot(e2, x = -6 .. 6, y = -6 .. 6, thickness = 3, color = blue):

  display([p4, Ell1, Ell2, Cir1, Cir2, Cir3, Cir4], scaling = constrained, axes = none);

 

 

 

Download hardWay.mw

 

when you specify the range o the independent variable as 0..3*Pi, then the calculation is (probably?) done in steps which are a fraction of 'Pi' - whihc means it is very likely to exactly "land on" the value Pi, at which point you have reached the limit specified for the range of the independent variable, so the calculation stops.

Two workarounds

  1. Don't supply a range for the dependent variable - for many ODEs you may not be able work this out a priori in any case
  2. Or make sure that the specified range of the dependet variable is "larger" than any value it can possibly achieve

Both methods are shown in the attached

  restart;
  interface(version);
  ode := diff(y(t), t$2) + y(t)=0;
#
# No specified range for the dependent variable
#
  DEtools:-DEplot( ode,
                   y(t),
                   t=0 .. 3*Pi,
                   [[y(0)=1, D(y)(0)=0]],
                   linecolor=blue
                 );
#
# "Large" specified range for the dependent variable
#
  DEtools:-DEplot( ode,
                   y(t),
                   t=0 .. 3*Pi,
                   y=-1.1..1.1,
                   [[y(0)=1, D(y)(0)=0]],
                   linecolor=blue
                 );

`Standard Worksheet Interface, Maple 2022.1, Windows 7, May 26 2022 Build ID 1619613`

 

diff(diff(y(t), t), t)+y(t) = 0

 

 

 

 

Download deplot.mw

to use the Maple's geometry() package. One advantage of this approach is that its draw() command has the concept of a stacking order for geometric objects. So, for example, in the attached

  1. the red square is at the bottom of the stack
  2. then two filled (with color=white) ellipses are placed on top of the red square
  3. then two ellipse "outlines" are placed on top of 1+2 above
  4. then four filled circles are placed on top of 1+2+3 above.

See the attached.

(As usual, the figure renders rather better in a Maple worksheet than it does on this site- honest!)

  restart;
  with(geometry):
  _EnvHorizontalName := 'x':
  _EnvVerticalName := 'y':
  point(OO,0,0):

  a := 7:
  b := a/2:
  r := b/2:
  c := 2*sqrt(10)*r:
#
# The "background" square
#
  square( Sq,
          [ point(sq1, -c/2, -c/2),
            point(sq2, -c/2,  c/2),
            point(sq3,  c/2,  c/2),
            point(sq4,  c/2, -c/2)
          ]
        ):
#
# The four circles
#
  circle( cir1,
          [ point(c1, c/2 - r, -c/2 + r),
            r
          ]
        ):
  circle( cir2,
          [ point(c2, c/2 - r,  c/2 - r),
            r
          ]
        ):
  circle( cir3,
          [ point(c3, -c/2 + r, -c/2 + r),
            r
          ]
        ):
  circle( cir4,
          [ point(c4, -c/2 + r, c/2 - r),
            r
          ]
        ):
#
# The two ellipses
#
  ellipse(e1, x^2/a^2 + y^2/b^2 = 1):
  rotation(e2, e1, Pi/4, clockwise, OO):
  rotation(e3, e1, Pi/4, counterclockwise, OO):
#
# Draw everything - note that the "stacking"
# order is important
#
  draw( [ cir1( filled=true, color=blue),
          cir2( filled=true, color=blue),
          cir3( filled=true, color=blue),
          cir4( filled=true, color=blue),
          e2( color=blue, thickness = 3),
          e3( color=blue, thickness = 3),
          e2( filled=true, color=white),
          e3( filled=true, color=white),
          Sq( filled=true, color=red)
        ],
        axes=none,
        size=[800, 800],
        gridlines=false
      );
 

 

 

Download geomDraw.mw

 

3 4 5 6 7 8 9 Last Page 5 of 191