tomleslie

Same error...

Answered: tomleslie 13876

November 10 2022

0 0

in

Maple 2021
Maple 2020

No error in Maple 2019, although not the most *useful* solution I've seen - see attached

>	interface(version);

(1)

>

restart;

>

ode := diff(f(Z), Z$2) = 3602879701896397*(3860741894751515/1125899906842624 + (5048392920194975*f(Z)^2)/1073741824 - (2484212754397225*f(Z)^4)/1073741824 - (321579057930643*f(Z)^6)/2147483648 - (4936153155163025*f(Z)^8)/2251799813685248)/(4611686018427387904*f(Z)^3*(4178268760908915/1073741824 + 315761000*f(Z)^2));
dsolve([ode,f(0)=1,D(f)(0)=0])

diff(diff(f(Z), Z), Z) = (3602879701896397/4611686018427387904)*(3860741894751515/1125899906842624+(5048392920194975/1073741824)*f(Z)^2-(2484212754397225/1073741824)*f(Z)^4-(321579057930643/2147483648)*f(Z)^6-(4936153155163025/2251799813685248)*f(Z)^8)/(f(Z)^3*(4178268760908915/1073741824+315761000*f(Z)^2))

(2)

>

Download odeErr.mw

You seem...

Answered: tomleslie 13876

November 09 2022

4 0

to have picked a "non-robust" method of computing tangent planes.

A couple of "better" methods are shown in the attached

>	# # From Google translate # # Given the hemisphere with the equation z=sqrt(4-x^2-y^2). # Verify that the given point lies on the hemisphere and # find the equation of the tangent plane to the hemisphere # at the given point #

>

  restart;
  with(Student[MultivariateCalculus]):
  surf:= z^2+x^2+y^2=4:
  tp:= pt->expand(Gradient((rhs-lhs)(surf), [x,y,z]=convert(pt, list))[]^%T.(<x,y,z>-pt))=0:
  onSur:= pt-> evalb(eval(surf, [x=pt[1], y=pt[2], z=pt[3]])):
#
# Check point is on sphere
#
  onSur( <1,1,sqrt(2)> );
#
# Compute tangent plane
#
  tp( <1,1,sqrt(2)> );
#
# Check point is on sphere
#
  onSur( <sqrt(2), sqrt(2), 0> );
#
# Compute tangent plane
#
  tp( <sqrt(2), sqrt(2), 0> );
#
# Check point is on sphere
#
  onSur( <2,0,0> );
#
# Compute tangent plane
#
  tp( <2,0,0> );

(1)

>

  restart;
  with(geom3d):
  _EnvXName := 'x': _EnvYName := 'y': _EnvZName := 'z':
  sphere(sph, z^2+x^2+y^2=4):
  point(P1, [1, 1, sqrt(2)]):
#
# Check point is on sphere
#
  evalb(eval( Equation(sph), [x=xcoord(P1), y=ycoord(P1), z=zcoord(P1)]));
#
# Check point is on sphere
#
  point(P2, [sqrt(2), sqrt(2), 0]):
  evalb(eval( Equation(sph), [x=xcoord(P2), y=ycoord(P2), z=zcoord(P2)]));
#
# Check point is on sphere
#
  point(P3, [2, 0, 0]):
  evalb(eval( Equation(sph), [x=xcoord(P3), y=ycoord(P3), z=zcoord(P3)]));
#
# Compute tangent planes
#
  Equation( TangentPlane(pl1, P1, sph));
  Equation( TangentPlane(pl2, P2, sph));
  Equation( TangentPlane(pl3, P3, sph));
#
# Draw everything - just for fun
#
  draw( [ sph(color=red, style=surface),
          P1(color=blue, symbol=solidsphere, symbolsize=20),
          P2(color=green, symbol=solidsphere, symbolsize=20),
          P3(color=yellow, symbol=solidsphere, symbolsize=20),
          pl1(color=blue, style=surface, transparency=0.4),
          pl2(color=green, style=surface, transparency=0.4),
          pl3(color=yellow, style=surface, transparency=0.4)
        ]
      );

>

Download tanPlanes.mw

Something similar...

Answered: tomleslie 13876

November 08 2022

1 0

to the ecxample shown in the attached, maybe?

>

  restart;
  with(plots):
  doCol:= (c1, c2, c3)-> display
                         ( [ arrow
                             ( <0,0,255>,
                                 color=ColorTools:-Color( "RGB", [0,0,255])
                             ),
                             arrow
                             ( <0,255,0>,
                                color=ColorTools:-Color( "RGB", [0,255,0])
                             ),
                             arrow
                             ( <255,0,0>,
                                color=ColorTools:-Color( "RGB", [255,0,0])
                             ),
                             arrow
                             ( <c1,c2,c3>,
                                color=ColorTools:-Color( "RGB", [c1,c2,c3])
                             ),
                             textplot3d
                             ( [ c1, c2, c3,
                                 StringTools:-Join(convert~([c1, c2, c3], string), ", ")
                               ],
                               align=right,
                               font=[times, cold, 20]
                             )
                           ]
                         ):
#
# An example
#
  doCol(217, 43, 143);

>

Download rgbplt.mw

Well...

Answered: tomleslie 13876

November 08 2022

1 0

I haven't got Maple 17 (released 2013) installed - too old, even for me. Note that in your worksheet, with output=plot, only five iterations are used and the figure caption includes the statement

"The stopping criterion is not met"

The help for this command does not make clear how many iterations will be performed for each of the output options, so I assumethat with output=plot, Maple tries five iterations, then gives up. You can get around this by using the maxiterations option along with output=plot option, ie the code

Newton(P(x), x = 0.5, output = plot, tolerance = 10^(-6), stoppingcriterion = function_value, maxiterations = 10)

This shows that the stopping criterion is met after 6 iterations - whihc is the same as the command with the output=animation option

Why does output=plot only try five iterations, when output=animation apparently tries more - no idea!

Not clear to me...

Answered: tomleslie 13876

November 08 2022

1 4

exactly what you want - maybe the final plot in the attached?

Download contPlot2.mw

Hmmm...

Answered: tomleslie 13876

November 07 2022

0 0

A few points

'print' statements do not return anything from a procedure - you should be using a return statement
if you wish to return a vector, why do you (implicitly) define it as a table? Wouldn't it be a better idea to (explcitly) define it as a vector?
In the attached I fixed the above, and also deleted a few things that you don't need. The output is a four column matrix, mainly because the procedure exited on the first iteration, and I wanted to check that the error criterion was satisfied.
I have made no attempt to check your algorithm, but I have added an execution group at the end which returns the real eigenvalues and eigenvectors, using Maple's built-in LinearAlgebra:-Eigenvectors() command. The values returned by your procedure do not seem to agree very well with this :-(

ev.mw

The basic problem...

Answered: tomleslie 13876

November 06 2022

1 0

with your original worksheet is that in the definition of rot_mat, you have used the command MatrixPower(), which should have been LinearAlgebra[MatrixPower](). This is fixed in the attached.

You can avoid this issue (and make you code a little cleaner) by using the commands with(LinearAlgebra) and with(ArrayTools), so that you don't have to use 'long-form' package commands.

Both approaches are are shown in the attached

$proc (x::Vector) options operator, arrow; rtable(1 .. 3, 1 .. 3, {(1, 2) = -x[3], (1, 3) = x[2], (2, 1) = x[3], (2, 3) = -x[1], (3, 1) = -x[2], (3, 2) = x[1]}, datatype = anything, subtype = Matrix, storage = rectangular, order = Fortran_order) end proc$

(1)

Matrix(%id = 36893488147927011804)

(2)

$proc (x::Vector, theta) options operator, arrow; rtable(1 .. 3, 1 .. 3, {(1, 1) = 1+(1-cos(theta))*(-x[2]^2-x[3]^2), (1, 2) = -sin(theta)*x[3]+(1-cos(theta))*x[2]*x[1], (1, 3) = sin(theta)*x[2]+(1-cos(theta))*x[3]*x[1], (2, 1) = sin(theta)*x[3]+(1-cos(theta))*x[2]*x[1], (2, 2) = 1+(1-cos(theta))*(-x[1]^2-x[3]^2), (2, 3) = -sin(theta)*x[1]+(1-cos(theta))*x[3]*x[2], (3, 1) = -sin(theta)*x[2]+(1-cos(theta))*x[3]*x[1], (3, 2) = sin(theta)*x[1]+(1-cos(theta))*x[3]*x[2], (3, 3) = 1+(1-cos(theta))*(-x[1]^2-x[2]^2)}, datatype = anything, subtype = Matrix, storage = rectangular, order = Fortran_order) end proc$

(3)

Matrix(%id = 36893488147962579172)

(4)

trans_mat_ang := unapply(ArrayTools[Concatenate](1, ArrayTools[Concatenate](2, rot_mat(`<,>`(x[1], x[2], x[3]), theta), (LinearAlgebra[IdentityMatrix](3)-rot_mat(`<,>`(x[1], x[2], x[3]), theta)).`<,>`(p[1], p[2], p[3])), `<|>`(0, 0, 0, 1)), x::Vector, theta, p::Vector)

$proc (x::Vector, theta, p::Vector) options operator, arrow; rtable(1 .. 4, 1 .. 4, {(1, 1) = 1+(1-cos(theta))*(-x[2]^2-x[3]^2), (1, 2) = -sin(theta)*x[3]+(1-cos(theta))*x[2]*x[1], (1, 3) = sin(theta)*x[2]+(1-cos(theta))*x[3]*x[1], (1, 4) = -(1-cos(theta))*(-x[2]^2-x[3]^2)*p[1]+(sin(theta)*x[3]-(1-cos(theta))*x[2]*x[1])*p[2]+(-sin(theta)*x[2]-(1-cos(theta))*x[3]*x[1])*p[3], (2, 1) = sin(theta)*x[3]+(1-cos(theta))*x[2]*x[1], (2, 2) = 1+(1-cos(theta))*(-x[1]^2-x[3]^2), (2, 3) = -sin(theta)*x[1]+(1-cos(theta))*x[3]*x[2], (2, 4) = (-sin(theta)*x[3]-(1-cos(theta))*x[2]*x[1])*p[1]-(1-cos(theta))*(-x[1]^2-x[3]^2)*p[2]+(sin(theta)*x[1]-(1-cos(theta))*x[3]*x[2])*p[3], (3, 1) = -sin(theta)*x[2]+(1-cos(theta))*x[3]*x[1], (3, 2) = sin(theta)*x[1]+(1-cos(theta))*x[3]*x[2], (3, 3) = 1+(1-cos(theta))*(-x[1]^2-x[2]^2), (3, 4) = (sin(theta)*x[2]-(1-cos(theta))*x[3]*x[1])*p[1]+(-sin(theta)*x[1]-(1-cos(theta))*x[3]*x[2])*p[2]-(1-cos(theta))*(-x[1]^2-x[2]^2)*p[3], (4, 4) = 1}, datatype = anything, subtype = Matrix, storage = rectangular, order = Fortran_order) end proc$

(5)

$proc (x::Vector, p::Vector, theta) options operator, arrow; rtable(1 .. 3, 1 .. 4, {(1, 1) = 1+(1-cos(theta))*(-x[2]^2-x[3]^2), (1, 2) = -sin(theta)*x[3]+(1-cos(theta))*x[2]*x[1], (1, 3) = sin(theta)*x[2]+(1-cos(theta))*x[3]*x[1], (1, 4) = -(1-cos(theta))*(-x[2]^2-x[3]^2)*p[1]+(sin(theta)*x[3]-(1-cos(theta))*x[2]*x[1])*p[2]+(-sin(theta)*x[2]-(1-cos(theta))*x[3]*x[1])*p[3], (2, 1) = sin(theta)*x[3]+(1-cos(theta))*x[2]*x[1], (2, 2) = 1+(1-cos(theta))*(-x[1]^2-x[3]^2), (2, 3) = -sin(theta)*x[1]+(1-cos(theta))*x[3]*x[2], (2, 4) = (-sin(theta)*x[3]-(1-cos(theta))*x[2]*x[1])*p[1]-(1-cos(theta))*(-x[1]^2-x[3]^2)*p[2]+(sin(theta)*x[1]-(1-cos(theta))*x[3]*x[2])*p[3], (3, 1) = -sin(theta)*x[2]+(1-cos(theta))*x[3]*x[1], (3, 2) = sin(theta)*x[1]+(1-cos(theta))*x[3]*x[2], (3, 3) = 1+(1-cos(theta))*(-x[1]^2-x[2]^2), (3, 4) = (sin(theta)*x[2]-(1-cos(theta))*x[3]*x[1])*p[1]+(-sin(theta)*x[1]-(1-cos(theta))*x[3]*x[2])*p[2]-(1-cos(theta))*(-x[1]^2-x[2]^2)*p[3]}, datatype = anything, subtype = Matrix, storage = rectangular, order = Fortran_order) end proc$

(6)

Matrix(%id = 36893488148141183572)

(7)

restart; with(LinearAlgebra); with(ArrayTools); unit_axis_cross_mat := unapply(`<,>`(`<|>`(0, -x[3], x[2]), `<|>`(x[3], 0, -x[1]), `<|>`(-x[2], x[1], 0)), x::Vector)

$proc (x::Vector) options operator, arrow; rtable(1 .. 3, 1 .. 3, {(1, 2) = -x[3], (1, 3) = x[2], (2, 1) = x[3], (2, 3) = -x[1], (3, 1) = -x[2], (3, 2) = x[1]}, datatype = anything, subtype = Matrix, storage = rectangular, order = Fortran_order) end proc$

(8)

Matrix(%id = 36893488147926986612)

(9)

$proc (x::Vector, theta) options operator, arrow; rtable(1 .. 3, 1 .. 3, {(1, 1) = 1+(1-cos(theta))*(-x[2]^2-x[3]^2), (1, 2) = -sin(theta)*x[3]+(1-cos(theta))*x[2]*x[1], (1, 3) = sin(theta)*x[2]+(1-cos(theta))*x[3]*x[1], (2, 1) = sin(theta)*x[3]+(1-cos(theta))*x[2]*x[1], (2, 2) = 1+(1-cos(theta))*(-x[1]^2-x[3]^2), (2, 3) = -sin(theta)*x[1]+(1-cos(theta))*x[3]*x[2], (3, 1) = -sin(theta)*x[2]+(1-cos(theta))*x[3]*x[1], (3, 2) = sin(theta)*x[1]+(1-cos(theta))*x[3]*x[2], (3, 3) = 1+(1-cos(theta))*(-x[1]^2-x[2]^2)}, datatype = anything, subtype = Matrix, storage = rectangular, order = Fortran_order) end proc$

(10)

Matrix(%id = 36893488148113533156)

(11)

trans_mat_ang := unapply(Concatenate(1, Concatenate(2, rot_mat(`<,>`(x[1], x[2], x[3]), theta), (IdentityMatrix(3)-rot_mat(`<,>`(x[1], x[2], x[3]), theta)).`<,>`(p[1], p[2], p[3])), `<|>`(0, 0, 0, 1)), x::Vector, theta, p::Vector)

$proc (x::Vector, theta, p::Vector) options operator, arrow; rtable(1 .. 4, 1 .. 4, {(1, 1) = 1+(1-cos(theta))*(-x[2]^2-x[3]^2), (1, 2) = -sin(theta)*x[3]+(1-cos(theta))*x[2]*x[1], (1, 3) = sin(theta)*x[2]+(1-cos(theta))*x[3]*x[1], (1, 4) = -(1-cos(theta))*(-x[2]^2-x[3]^2)*p[1]+(sin(theta)*x[3]-(1-cos(theta))*x[2]*x[1])*p[2]+(-sin(theta)*x[2]-(1-cos(theta))*x[3]*x[1])*p[3], (2, 1) = sin(theta)*x[3]+(1-cos(theta))*x[2]*x[1], (2, 2) = 1+(1-cos(theta))*(-x[1]^2-x[3]^2), (2, 3) = -sin(theta)*x[1]+(1-cos(theta))*x[3]*x[2], (2, 4) = (-sin(theta)*x[3]-(1-cos(theta))*x[2]*x[1])*p[1]-(1-cos(theta))*(-x[1]^2-x[3]^2)*p[2]+(sin(theta)*x[1]-(1-cos(theta))*x[3]*x[2])*p[3], (3, 1) = -sin(theta)*x[2]+(1-cos(theta))*x[3]*x[1], (3, 2) = sin(theta)*x[1]+(1-cos(theta))*x[3]*x[2], (3, 3) = 1+(1-cos(theta))*(-x[1]^2-x[2]^2), (3, 4) = (sin(theta)*x[2]-(1-cos(theta))*x[3]*x[1])*p[1]+(-sin(theta)*x[1]-(1-cos(theta))*x[3]*x[2])*p[2]-(1-cos(theta))*(-x[1]^2-x[2]^2)*p[3], (4, 4) = 1}, datatype = anything, subtype = Matrix, storage = rectangular, order = Fortran_order) end proc$

(12)

$proc (x::Vector, p::Vector, theta) options operator, arrow; rtable(1 .. 3, 1 .. 4, {(1, 1) = 1+(1-cos(theta))*(-x[2]^2-x[3]^2), (1, 2) = -sin(theta)*x[3]+(1-cos(theta))*x[2]*x[1], (1, 3) = sin(theta)*x[2]+(1-cos(theta))*x[3]*x[1], (1, 4) = -(1-cos(theta))*(-x[2]^2-x[3]^2)*p[1]+(sin(theta)*x[3]-(1-cos(theta))*x[2]*x[1])*p[2]+(-sin(theta)*x[2]-(1-cos(theta))*x[3]*x[1])*p[3], (2, 1) = sin(theta)*x[3]+(1-cos(theta))*x[2]*x[1], (2, 2) = 1+(1-cos(theta))*(-x[1]^2-x[3]^2), (2, 3) = -sin(theta)*x[1]+(1-cos(theta))*x[3]*x[2], (2, 4) = (-sin(theta)*x[3]-(1-cos(theta))*x[2]*x[1])*p[1]-(1-cos(theta))*(-x[1]^2-x[3]^2)*p[2]+(sin(theta)*x[1]-(1-cos(theta))*x[3]*x[2])*p[3], (3, 1) = -sin(theta)*x[2]+(1-cos(theta))*x[3]*x[1], (3, 2) = sin(theta)*x[1]+(1-cos(theta))*x[3]*x[2], (3, 3) = 1+(1-cos(theta))*(-x[1]^2-x[2]^2), (3, 4) = (sin(theta)*x[2]-(1-cos(theta))*x[3]*x[1])*p[1]+(-sin(theta)*x[1]-(1-cos(theta))*x[3]*x[2])*p[2]-(1-cos(theta))*(-x[1]^2-x[2]^2)*p[3]}, datatype = anything, subtype = Matrix, storage = rectangular, order = Fortran_order) end proc$

(13)

Matrix(%id = 36893488148113533756)

(14)

Download HomTrans.mw

Probaby simpler...

Answered: tomleslie 13876

November 06 2022

2 1

to use the option 'useint' in the dsolve() command, as shown in the attached

Download odeint.mw

Use a plot array...

Answered: tomleslie 13876

November 06 2022

2 10

The code in the attached, uses

Q=[0.5617, 0.4392, 0.2564, 0.1645, 0.0659]
lambda=[0.0096, 0.0010, 0.0093, 0.0075, 0.0031]

to produce an array of 25 plots, organised as 5x5. Each individual plot contains nine curves, given by

k = 0.1*j, We = 0.1, Br = 0.3, x = 0 with j=1..9

So far as I can tell, you want a further five such plot arrays based on diifferent parameter variations - these are left as an exercise

Plot arrays do not display on this website, so you will have to download the following

pltArr.mw

The solve() command...

Answered: tomleslie 13876

November 04 2022

1 0

in your procedure is returning non-numeric values. The solve() command returns non-numeric values, whenever

[[i = 20., x = x, y = 0.], [i = 40., x = x, y = 0.], [i = 60., x = x, y = 0.], [i = 80., x = x, y = 0.]]

where 'i' is the index in the command

display([seq(Fig((2*Pi*i)/40), i = 1 .. 80)], insequence = true);

In other words the solve() command fails whenerver the argument passed to the procedure is an integer multiple of Pi

Many issues...

Answered: tomleslie 13876

November 03 2022

0 3

You can't usee square brackets, ie '[]' to group terms in Maple. [] brackets are used for lists. Fixed in the attached
Yo attempt to assign to the variable 'define' whichi is a Maple keyword. The only way to do this is use the 'local define' command,- fixed in the attached.
You want to evalue the ODE de1 with values give by 'testExps'. However the latter is defined nowhere in your worksheet - commented out in the attached
I think it unlikely that you will be able to obtain a usable 'power series' in the variable 'we', because so many terms containing 'we' are raised to the power (n-1)/2 where n is unknown. Now if 'n' is known, and is an odd integer, it may be possible

Fixes (so far) are shown in the attached

>

Warning, A new binding for the name `define` has been created. The global instance of this name is still accessible using the :- prefix, :-`define`. See ?protect for details.

>

2*(1+we*(diff(u(y), y))^2)^((1/2)*n-1/2)*((1/2)*n-1/2)*we*(diff(u(y), y))*(diff(diff(u(y), y), y))*u(y)/(1+we*(diff(u(y), y))^2)+(1+we*(diff(u(y), y))^2)^((1/2)*n-1/2)*(diff(u(y), y)) = diff(p(x), x)

(1)

>

define := [u(y) = u__0(y)+we*u__1(y)+we^2*u__2(y), p(x) = p__0(x)+we*p__1(x)+we^2*p__2(x)]

(2)

>

eval(ode, define)

(3)

>

de2 := collect(lhs(%)-rhs(%), delta)

(4)

>

coeff(de2, we, 0)

Error, unable to compute coeff

>

"Extract the ODE from de2 which is the #` coefficient of (we)^1, then substitute the value` #` for u__0(y) from sol1 above ` "

Error, unable to parse

>

Download serode.mw

Agree with Rouben...

Answered: tomleslie 13876

November 02 2022

2 1

but you also seem to believe that the variables 'phi' and 'varphi' are identical - they are not!

Fixing all of the above leads to the attached

Download tayl.mw

Wow...

Answered: tomleslie 13876

November 02 2022

1 1

Document Mode + 2D Input + Units

You do like pile up error-prone complexity, just to make things look 'pretty'

Anyhow, maybe the attached is what you were hoping for. The double square brackets ( '[[..]]') whihc appear around the units are a quirk of this site - they do not appear in the actual worksheet, honest

Download unitsAgain.mw

Maybe...

Answered: tomleslie 13876

November 02 2022

2 0

you were hoping to chieve something similar to the animation shown in the attached?

In future please use the big, green up-arrow in the Mapleprimes toolbar to upload a actual worksheet. When you cut-and-paste as in your original question, responders have too go through it line-by-line to work out what is Maple inpu, what is Maple output, what is comment, what is plaintext etc etc. This is very tedious!

>

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

EQ := proc(M, N)
           local eq;
           eq := (y - M[2])/(x - M[1]) = (N[2] - M[2])/(N[1] - M[1]);
      end proc:

#On considère un trapèze dans lequel une base est fixe l'autre base a une longueur constante et la somme des 2 autres côtés est constaante.
#Trouver :
#1-. le lieu des sommets mobiles.
  A := [xA, 0]:
  B := [xA + a, 0]:
  D := [xD, yD]:
  C := [xD + b, yD]:
  EQ(B, C):
  E := [xA + a - b, 0]:
  Omega1 := (A + E)/2:
#Application numérique :
#Lieux des sommets C et D

  xA := -5:
  a := 13:#a>=b
  b := 7:
  c := -3:
  xD := -6:
  xC := xD + c:

  A:
  B:
  C:
  D:
  Ll:=11:
  aa:= Ll/2:
  cc := (a - b)/2:
  bb := sqrt(aa^2 - cc^2):
  el1 := (x - Omega1[1])^2/aa^2 + y^2/bb^2 = 1:
  sol := solve(subs(x = xD, (x - Omega1[1])^2/aa^2 + y^2/bb^2 = 1), y):
  yD := sol[1]:
  el2 := (x - Omega1[1] - b)^2/aa^2 + y^2/bb^2 = 1:
  EL1 := implicitplot(el1, x = -9 .. 4, y = -6 .. 6, color = blue):
  EL2 := implicitplot(el2, x = -9 .. 12, y = -6 .. 12, color = blue):
  Trap := polygonplot([A, B, C, D], color = blue, filled = true, transparency = 0.9):
  Points := pointplot([A[], B[], C[], C[], D[], E[], Omega1[]], symbol = solidcircle, color = [red], symbolsize = 6):
  tex := textplot([[A[], "A"], [B[], "B"], [C[], "C"], [D[], "D"], [E[], "E"], [Omega1[], "Ω1"]], align = ["above", "right"]):
  display([Trap, EL1, EL2, tex, Points], axes = normal, scaling = constrained):
  Fig := proc(xD)
              local yD, D, C,Points,tex;
              global A, B, b, Omega1, xA, xB, EL1, EL2;
              subs(x = xD, (x - Omega1[1])^2/aa^2 + y^2/bb^2 = 1);
              solve(subs(x = xD, (x - Omega1[1])^2/aa^2 + y^2/bb^2 = 1), y);
              yD := %[1]; D:= [xD, yD]; C := [xD + b, yD];
              Points := pointplot([A[], B[], C[], C[], D[], E[], Omega1[]], symbol = solidcircle, color = [red], symbolsize = 6):
              tex := textplot([[A[], "A"], [B[], "B"], [C[], "C"], [D[], "D"], [E[], "E"], [Omega1[], "Ω1"]], align = ["above", "right"]):
              display([polygonplot([A, B, C, D], color = blue, filled = true, transparency = 0.9), Points,tex,EL1, EL2], axes = normal, scaling = constrained);
         end proc:

Fig(2):
Fig(-4):
display( [seq(Fig(-6 + 3*i/10), i = 1..20)], insequence=true);