Advanced Search

tomleslie

13836 Reputation

20 Badges

14 years, 342 days
Contact tomleslie

MaplePrimes Activity


These are replies submitted by tomleslie

I suspect...

tomleslie 13836
September 22 2022
0 0

@JAMET 

that this code does not produce what you think it does - where the lines L1 and L2 are tangent to the circles C1, C2 and the ellipse p1. As you have constructed it, neither L1 nor \L2 is tangent to the ellipse p1.

The "zoomed-in" plots in the first file attached show that L1 does not touch the ellipse at all, and L2 cuts the ellipse in two places - very close, but not close enough for a prize!!!

The second file attached shows my caculation of the same general construction, but this time the lines are tangent to both circles and the ellipse

> 

  restart:
  with(plots):
  with(geometry):
  _EnvHorizontalName := x:
  _EnvVerticalName := y:
  R := 11:
  r := 7:
  a := sqrt(R*r):

  b := 2:
  circle(C1, [point(P1, [0, 0]), R]):
  circle(C2, [point(P2, [R + 2*b + r, 0]), r]):
  ellipse(p, (x - R - b)^2/b^2 + y^2/a^2 = 1):
  draw([C1(color = yellow, filled = true),
  C2(color = red, filled = true), p(color = blue, filled = true),
  C1(color = black), C2(color = black), p(color = black)],
  axes = none, view = [-15 .. 35, -15 .. 15], scaling = constrained):
  alpha := arctan((R - r)/(R + 2*b + r)):
  long := cos(alpha)*(R + 2*b + r):
  evalf(%):
  circle(C2, [point(P2, [long, r - R]), r]):
  rotation(p1, p, alpha, 'clockwise'):
  detail(p1):
  point(A, 0, -R):
  point(B, long, -R):
  line(L1, [A, B]):
  point(cen, [(143*sqrt(5))/25, -(26*sqrt(5))/25]):
  line(L, [P1, P2]):
  reflection(L2, L1, L):
  detail(L2):
  draw( [ C1(color = yellow, filled = true),
          C2(color = red, filled = true),
          p1(color = blue, filled = true),
          C1(color = black),
          C2(color = black),
          p1(color = black),
          L1(color = black),
          L2(color = black)
        ],
        axes = none,
        view = [-15 .. 35, -15 .. 15],
        scaling = constrained
      );
  draw([p1(color = blue, filled = true, transparency=0.5, numpoints=1000),
        p1(color = black, numpoints=1000),
        L2(color = black)
       ],
       axes = none,
       view = [14..15, 6..6.5],
       scaling = constrained,
       title="Zoomed plot of L2 and p1",
       titlefont=[times, bold, 20]
     );
  draw([p1(color = blue, filled = true, transparency=0.5, numpoints=1000),
        p1(color = black, numpoints=1000),
        L1(color = black)
       ],
       axes = none,
       view = [10..12, -11.5..-10.5],
       scaling = constrained,
       title="Zoomed plot of L1 and p1",
       titlefont=[times, bold, 20]
     );

 

 

 
> 

 

Download badGeo.mw

Corrected???

> 

  restart:
  with(plots):
  with(geometry):
  _EnvHorizontalName := x:
  _EnvVerticalName := y:
  R := 11:
  r := 7:
  a:=sqrt(R*r):
  b := 2:
  circle(C1, [point(P1, [0, 0]), R]):
  circle(C2, [point(P2, [R + 2*b + r, 0]), r]):
#
# Construct external tangent lines
#
  circle( C3, [P1, R-r]):
  circle( C4, [ point
                ( P4,
                  1/2*~coordinates(center(C2))
                ),
                1/2*coordinates(center(C2))[1]
              ]
        ):
  intersection( Px, C3, C4, [Px1, Px2]):
  line(L1, [center(C1), Px1]):
  line(L2, [center(C1), Px2]):
  intersection( Py, L1, C1, [Py1, Py2]):
  intersection( Pz, L2, C1, [Pz1, Pz2]):
  PerpendicularLine( t1,
                     `if`( coordinates(Py1)[1]>0,Py1,Py2),
                     L1
                   ):
  PerpendicularLine( t2,
                     `if`( coordinates(Pz1)[1]>0,Pz1,Pz2),
                     L2
                   ):
#
# Compute the ellipse
#
  ellipse( e1,
           (x - R - b)^2/b^2 + y^2/a^2 = 1
          ):
#
# Rotation
#
  angle:=arctan(solve(eval(Equation(t2), x=1),y)-solve( eval(Equation(t2), x=0),y)):
  rotation( t2r,
            t2,
            angle,
            clockwise,
            P1
          ):
  rotation( t1r,
            t1,
            angle,
            clockwise,
            P1
          ):
  rotation( C2r,
            C2,
            angle,
            clockwise,
            P1
          ):
  rotation( e1r,
            e1,
            angle,
            clockwise,
            P1
          ):
#
# Draw everything
#
  draw( [ C1(color=yellow, filled=true),
          C2r(color=red, filled=true),
          C1(color = black),
          C2r(color = black),
          t1r(color=black),
          t2r(color=black),
          e1r(color=blue, filled=true, numpoints=1000)
        ],
        axes=none,
        scaling=constrained
      );
#
# Zoom in on tangents to the ellipse
#
  draw( [ t1r(color=black),
          e1r(color=blue,filled=true, numpoints=1000)
        ],
        axes=none,
        scaling=constrained,
        view=[13..15, 5..7],
        title="Zoomed plot of t1r and e1r",
        titlefont=[times, bold, 20]
      );
  draw( [ t2r(color=black),
          e1r(color=blue,filled=true, numpoints=1000)
        ],
        axes=none,
        scaling=constrained,
        view=[10.5..12.5, -11..-9],
        title="Zoomed plot of t2r and e1r",
        titlefont=[times, bold, 20]
      );

 

 

 
> 

 

Download geoRFl4.mw

Everything...

tomleslie 13836
September 21 2022
0 0

@C_R 

Everything in my original response is documented - I'm really not doing anything clever!

From the help page for indets() command

Calling Sequence
indets(expr)
indets(expr, typename)
indets[flat](expr)
indets[flat](expr, typename)

I simply used the second of these calling sequences, highlighted in red above.

And 'suffixed is a type. From the help page for type/suffixed

Calling Sequence
type(x, suffixed(y))
type(x, suffixed(y, stype))

again I used the highlighted sequence; Later the same help page says (emphasis added)

With two parameters, x is checked to see if it has the same type of and begins with the characters of y and that the balance, when parsed, is of type stype.

so suffixed(t, nonnegint) matches any name which starts with a 't' and is followed by a non-negative integer

Is there a reason...

tomleslie 13836
September 17 2022
0 0

why you don't want the exact answer, as shown in the attached

> 

  restart;
  odeSys:= diff(i(t),t)=20000-50*i(t)-150*w(t),
           diff(w(t),t)=1.8*i(t)-29.16;
  inits:=w(0)=0, i(0)=0;

  ans:=dsolve([odeSys, inits]);
  plot( [ eval(i(t), ans),
          eval(w(t), ans)
        ],
        t=0..1,
        color=[red, blue],
        legend=[i(t), w(t)]
      );

diff(i(t), t) = 20000-50*i(t)-150*w(t), diff(w(t), t) = 1.8*i(t)-29.16

  

w(0) = 0, i(0) = 0

  

{i(t) = (-(5/9)*355^(1/2)-125/9)*((121031/26625)*355^(1/2)-1919/30)*exp(-(25+355^(1/2))*t)+81/5+((5/9)*355^(1/2)-125/9)*(-(121031/26625)*355^(1/2)-1919/30)*exp((-25+355^(1/2))*t), w(t) = exp((-25+355^(1/2))*t)*(-(121031/26625)*355^(1/2)-1919/30)+exp(-(25+355^(1/2))*t)*((121031/26625)*355^(1/2)-1919/30)+1919/15}

  

 
> 

 

Download simpleODE.mw

Well...

tomleslie 13836
September 16 2022
0 0

the equation

y(x)=x*f(diff(y(x),x))+g(diff(y(x),x))

for arbitrary functions f() and g() -  so what about the particular case f()=0, g() arbitrary. Is that a D'Alembert equation or not?

*Edit*

Another interesting question - if one has the particular case f() arbitrary, g()=0. Is that a D'Alembert equation or not?

Wrong question...

tomleslie 13836
September 15 2022
0 0

@PhearunSeng 

You state

". eqs is generated by by another equaction, the question is how to convert......"

In other words,

You have produced crap and I would like someone to fix it

Have you even considered the correct approach - which is

Don't produce crap in the first place

Your procedure...

tomleslie 13836
September 14 2022
0 0

@ijuptilk 

doCalc(j) returns three curves (call these c1, c2, c3) in a single plot for each value of j. For the "stipped down example I supplied with only two values of j, the ordering in myData is

 <--------------------  j=1  ------------------------>|<----------------------- j=2  ---------------------->

c1(x)   c1(y)   c2(x)   c2(y)    c3(x)    c3(y)  |  c1(x)   c1(y)   c2(x)   c2(y)   c3(x)   c3(y)

Rather obviously...

tomleslie 13836
September 14 2022
0 0

@PhearunSeng 

you have a restart command, which crears Maple memory. You then attempt to use 'eqs', whihc has not been defined

Hmmmm...

tomleslie 13836
September 14 2022
0 0

The only "ExcelTools()" command in the worksheet you supply is

#ExcelTools:-Export(op([1, 1], plot(theta_f, t=0..100)));

and it is commented out. Out of interest, the argument for the above command, ie op([1, 1], plot(theta_f, t=0..100)), returns a 200X2 matrix, which ExcelTools:-Export() will happily export.

Fixed...

tomleslie 13836
September 12 2022
0 0

the syntax errors in your latest upload. A lot of these were te same as I fixed in the answer to your original upload.

See the attached

 


 

restart

with(VariationalCalculus)NULL

with(plots)
with(linalg)

We will now declare our first equations

 

 

x[1] := l[1]*sin(theta(t))

y[1] := -l[1]*cos(theta(t))``

x[2] := -l[2]*sin(theta(t)); y[2] := l[2]*cos(theta(t))

x[4] := x[1]-l[4]*sin(theta(t)+phi(t)); y[4] := y[1]+l[4]*cos(theta(t)+phi(t)); x[3] := x[2]+l[3]*sin(theta(t)-psi(t)); y[3] := y[2]-l[3]*cos(theta(t)-psi(t))

R[1] := vector(2, [x[1], y[1]]); R[2] := vector(2, [x[2], y[2]]); R[4] := vector(2, [x[4], y[4]]); R[3] := vector(2, [x[3], y[3]])

R[1] := map(diff, R[1], t); R[2] := map(diff, R[2], t); R[4] := map(diff, R[4], t); R[3] := map(diff, R[3], t)

NULL

NULL

T := combine(simplify(collect((1/2)*m[1]*innerprod(R[4], R[4])+(1/2)*m[2]*innerprod(R[3], R[3]), [l[1], l[2], l[3], l[4], m[1], m[2], diff(theta(t), t), diff(phi(t), t), diff(psi(t), t), cos(psi(t)), sin(psi(t))])), trig)

NULL

U := collect(simplify(expand(g*m[1]*y[4]+g*m[2]*y[3])), [l[1], l[2], l[3], l[4], m[1], m[2], diff(theta(t), t), diff(phi(t), t), diff(psi(t), t), cos(psi(t)), sin(psi(t)), g])

NULL

L := simplify(collect(T-U, [l[1], l[2], l[3], l[4], m[1], m[2], diff(theta(t), t), diff(phi(t), t), diff(psi(t), t), cos(psi(t)), sin(psi(t))]))

NULL

NULL

NULL

Eul := remove(has, EulerLagrange(L, t, [theta(t), phi(t), psi(t)]), K[1]); eq1 := op(select(has, Eul, (diff(psi(t), t))^2)); eq2 := op(select(has, remove(has, Eul, (diff(psi(t), t))^2), sin(psi(t)))); eq3 := op(remove(has, remove(has, Eul, (diff(psi(t), t))^2), sin(psi(t))))

NULL

NULL

INITS := theta(0) = 3*Pi*(1/4), (D(theta))(0) = 0, phi(0) = (1/4)*Pi, (D(phi))(0) = 0, psi(0) = (1/4)*Pi, (D(psi))(0) = 0

NULL

PARAM := [g = 9.8, l[1] = 10, l[2] = 100, l[3] = 100, l[4] = 21, m[1] = 1000, m[2] = 1]

NULL

Evaluating the general equation for motion with respect to our
parameters.

NULL

sys := (eval([eq1, eq2, eq3], PARAM))[]

NULL``

 

We use a numeric ODE solver to approximate the change in theta and
theta’ over time.

 

sol := dsolve([sys, INITS], numeric)

NULL

Plotting time vs angular velocity

 

odeplot(sol, [t, diff(psi(t), t)], 1 .. 3, numpoints = 1000, labels = ["times(sec)", "Psi(deg/s)"], labeldirections = [horizontal, vertical], font = [TIMES, ROMAN, 20])

  

odeplot(sol, [t, diff(theta(t), t)], t = 1 .. 3, numpoints = 1000, labels = ["time(sec)", "Theta(deg/s)"], labeldirections = [horizontal, vertical], font = [TIMES, ROMAN, 20])

  

NULL


Download test2_again.mw

 

Please...

tomleslie 13836
September 11 2022
0 0

upload your worksheet using the big green up-arrow in the Mapleprimes toolbar

Until you answer two simple questions...

tomleslie 13836
September 11 2022
0 0

you will make absolutely no progress. Questions are

  1. Which symbolic toolbox are you using within the Matlab environment?
  2. Which "symbolic formats" can the above toolbox import?

Which...

tomleslie 13836
September 10 2022
0 0

@mehdibgh 

"symbolic format" can Matlab import?

Until you specify this, no progress can be made

Difficult (for me)...

tomleslie 13836
September 07 2022
0 0

to work out what you are trying to achieve, but I think that the fundamental problem with you code is in the lines

circle(C5,[reflection(P5, center(C3), L2), radius(C2)]):

detail(C5):
circle(C6, [point(P5, [0, R/3 - (2*R)/9]), radius(C2)]):

where P5 is defined in the creation of C5, and then redefined in the definition of C6 - so both C5 and C6 will now have the same centre, which will be the last one defined above.

Fixing this, guessing a little about what you are trying to achieve, and tidying up a bit results in the attached

> 

  restart;
  with(plots):
  with(geometry):
  _EnvHorizontalName := x:
  _EnvVerticalName := y:
  R := 7:
  point(A, [0, R]):
  line(L1, y = sqrt(3)*x + R):
  line(L2, y = -sqrt(3)*x + R):
  line(L3, y = R/3):
  intersection(B, L1, L3):
  intersection(C, L2, L3):

  triangle(ABC, [A, B, C]):
  circle(C1, [point(P1, [0, 0]), R]):

  incircle(C2, ABC, 'centername'=P2);
  dsegment( ds1, [ point(sp1, [0, radius(C2)]), point(sp2, [0, -radius(C2)]) ]):

  circle(C3, [reflection(P3, center(C2), C),  radius(C2)]):
  circle(C4, [reflection(P4, center(C2), B),  radius(C2)]):
  circle(C5, [reflection(P5, center(C3), L2), radius(C2)]):
  circle(C6, [reflection(P6, center(C4), L1), radius(C2)]):

  translation(C7, C2, ds1):
  translation(C8, C7, ds1):
  translation(C9, C8, ds1):

  circle(Cir, [center(C8), 3*radius(C8)]):

  draw([ L1(color = blue),
         ABC(color = red, transparency = 0.5, filled = true),
         L2(color = blue), L3(color = blue),
         C1(color = blue, thickness = 3),
         C1(color = yellow, transparency = 0.8, filled = true),
         C2(color = blue, filled = true),
         C3(color = blue, filled = true),
         C4(color = blue, filled = true),
         C5(color = blue, filled = true),
         C6(color = blue, filled = true),
         C7(color = blue, filled = true),
         C8(color = blue, filled = true),
         C9(color = blue, filled = true),
         Cir(color = red,transparency = 0.5,filled=true)
       ],
       axes = normal,
       view = [-R .. R, -R .. R],
       scaling = constrained
     );

C2

  

 
> 

 

Download refLCir3.mw

 

FYI...

tomleslie 13836
September 07 2022
0 0
  1. I decide that this is a bug amd have reported it to Maplesoft, so (hopefully?!) it will be fixed in a future Maple release. However since there is such a simple workaround which achieves the same transformation (see my original worksheet), I did categorize this "bug" as minor
  2. There is a better(?) way to compute the "reference" circle C2, in your/my original worksheet - rather than having to supply explicit expressions for the centre and the radius, use the incircle() command appied to the triangle ABC - see the attached
> 

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

  R := 7:
  point(A, [0, R]):
  line(L1, y = sqrt(3)*x + R):
  line(L2, y = -sqrt(3)*x + R):
  line(L3, y = R/3):
  intersection(B, L1, L3):
  intersection(C, L2, L3):
  triangle(ABC, [A, B, C]):

  circle(C1, [point(P1, [0, 0]), R]):
  incircle(C2, ABC, 'centername'=P2):

  circle(C3, [ reflection(P3, center(C2), C), radius(C2)]):
  circle(C4, [ reflection(P4, center(C2), B), radius(C2)]):

  draw( [ L1(color = blue),
          ABC(color = red, transparency = 0.5, filled = true),
          L2(color = blue),
          L3(color = blue),
          C1(color = blue, thickness = 3),
          C1(color = yellow, transparency = 0.8, filled = true),
          C2(color = green, filled = true),
          C3(color = green, filled = true),
          C4(color = green, filled = true)
        ],
        axes = normal,
        view = [-R .. R, -R .. R],
        scaling = constrained
      );

 
> 

 
> 

 

Download reflCir2.mw

@David Mazziotti  The brackets [[ ...

tomleslie 13836
September 06 2022
0 0

@David Mazziotti 

The brackets [[ ]]  do not appear in the worksheet - they are an artefact of how this website renders "unit" quantities.

First 11 12 13 14 15 16 17 Last Page 13 of 207