tomleslie

6247 Reputation

17 Badges

10 years, 107 days

MaplePrimes Activity


These are answers submitted by tomleslie

restart:
coef7 := [-1, 2, alpha[1, 2], alpha[2, 6], (17*RootOf(64*_Z^3+80*_Z^2+1104*_Z+561)+17)/(alpha[1, 2]*alpha[2, 6]), -17/alpha[2, 6], -33/(32*RootOf(64*_Z^3+80*_Z^2+1104*_Z+561)), -(163/32+RootOf(64*_Z^3+80*_Z^2+1104*_Z+561)^2+5*RootOf(64*_Z^3+80*_Z^2+1104*_Z+561)*(1/4))/(alpha[1, 2]*alpha[2, 6]), -RootOf(64*_Z^3+80*_Z^2+1104*_Z+561)-5/4, RootOf(64*_Z^3+80*_Z^2+1104*_Z+561)]:
select(j-> not has(j, I), evalc~(allvalues~(coef7)));

[-1, 2, alpha[1, 2], alpha[2, 6], (-(17/24)*(11908+36*25678877^(1/2))^(1/3)+(13651/6)/(11908+36*25678877^(1/2))^(1/3)+119/12)/(alpha[1, 2]*alpha[2, 6]), -17/alpha[2, 6], -(33/32)/(-(1/24)*(11908+36*25678877^(1/2))^(1/3)+(803/6)/(11908+36*25678877^(1/2))^(1/3)-5/12), (-439/96-(-(1/24)*(11908+36*25678877^(1/2))^(1/3)+(803/6)/(11908+36*25678877^(1/2))^(1/3)-5/12)^2+(5/96)*(11908+36*25678877^(1/2))^(1/3)-(4015/24)/(11908+36*25678877^(1/2))^(1/3))/(alpha[1, 2]*alpha[2, 6]), (1/24)*(11908+36*25678877^(1/2))^(1/3)-(803/6)/(11908+36*25678877^(1/2))^(1/3)-5/6, -(1/24)*(11908+36*25678877^(1/2))^(1/3)+(803/6)/(11908+36*25678877^(1/2))^(1/3)-5/12]

(1)

#
# or in floats
#
  evalf~(%);

[-1., 2., alpha[1, 2], alpha[2, 6], 8.167082567/(alpha[1, 2]*alpha[2, 6]), -17./alpha[2, 6], 1.984763255, -4.714237667/(alpha[1, 2]*alpha[2, 6]), -.7304166203, -.5195833797]

(2)

 

Download realGet.mw

 which created this problem it is difficult to be specific.

In future I recommend that you use the big green up-arrow in the Mapleprimes toolbar to upload worksheets.

One thing I would try is,  rather than define a function as

E__max:=V->V/(r__1+ln(r__2/r__1);

use

E__max:=unapply( V/(r__1+ln(r__2/r__1), V);

The latter "evaluates" its first argument, before creating the function: the former does not. Read the help to fully appreciate the difference.

There may still be an issue about whether you are going to call this function with a "unitless" variable (ie a scalar) or always cally it with a quantity to which units have been assigned. I am pretty sure that the precise nature of the function definition will vary depending on which of these alternatives you prefer

 

to produce a family of curves satisfying the requirement is shown in the attached.

  restart:
#
# Going to look for a solution of the form
#
#         (x-1)*(x-2)*(x-3)
#         -----------------
#             a*x^2+b*x+c
#
# which certainly has the correct roots and
# at least has the potential to produce an
# asymptote of the required form
#
  p0:=a*x^2+b*x+c;
  p1:=(x-1)*(x-2)*(x-3);
  asymp:=2-x;
#
# Determine the value of the coefficient 'a' such
# that the gradient of the asymptote is satisfied
#
  m:= solve( [ limit(p1/(p0*x), x=infinity)=coeff(asymp, x, 1),
               limit(p1/(p0*x), x=-infinity)=coeff(asymp, x, 1)
             ],
             a
           );
#
# Determine the value of the coefficient 'b' such
# that the intercept of the asymptote is satisfied
#
  n:= solve( [ eval( limit( p1/p0-x/a, x=infinity)=coeff(asymp, x, 0), m),
               eval( limit( p1/p0-x/a, x=-infinity)=coeff(asymp, x, 0), m)
             ],
             b
           );
#
# Ensure that the denominator has no roots by selecting
# the value of the parameter 'c' so that the discriminant
# of the denominator is negative
#
  solve( eval(discrim(p0,x), [m[], n[]])<0, c);       

a*x^2+b*x+c

 

(x-1)*(x-2)*(x-3)

 

2-x

 

{a = -1}

 

{b = 4}

 

RealRange(-infinity, Open(-4))

(1)

#
# Generate some curves for values of the parameter
# 'c' in the range above
#
  curves:= [ seq( eval(p1/p0, [a=-1, b=4, c=j]),
                  j=-5..-25,-5
                )
           ];
#
# Verify the asymptotes of these curves
#
  Student[Calculus1]:-Asymptotes~(curves);
#
# Plot the curves, along with the desired asymptote
#
  plot( [ curves[],
          asymp
        ],
        x=-10..10
      );

[(x-1)*(x-2)*(x-3)/(-x^2+4*x-5), (x-1)*(x-2)*(x-3)/(-x^2+4*x-10), (x-1)*(x-2)*(x-3)/(-x^2+4*x-15), (x-1)*(x-2)*(x-3)/(-x^2+4*x-20), (x-1)*(x-2)*(x-3)/(-x^2+4*x-25)]

 

[[y = 2-x], [y = 2-x], [y = 2-x], [y = 2-x], [y = 2-x]]

 

 

 

Download asym.mw

 

to check whether or not the points are coplanar first? Something like the attached.

( I'd be a little wary if using floats to determine 'point' locations - same reason as why one shoeld never really check the equality of two floats).

  restart;
  with(geom3d):
  f:= (p,q,r,s, tname)->`if`( not AreCoplanar(p, q, r, s),
                              gtetrahedron(tname, [p,q,r,s] ),
                              `no_tetrahedron_possible`
                             ):

  f( point(A, 0, 0, 0), point(B, 1, 0, 0), point(C, 0, 0, 1), point(F, 0, 1, 0), t1):
  detail(t1);

Geom3dDetail(["name of the object", t1], ["form of the object", gtetrahedron3d], ["the 4 vertices", [[0, 0, 0], [1, 0, 0], [0, 0, 1], [0, 1, 0]]], ["the 4 faces", [[[0, 0, 0], [1, 0, 0], [0, 0, 1]], [[0, 0, 0], [1, 0, 0], [0, 1, 0]], [[0, 0, 0], [0, 0, 1], [0, 1, 0]], [[1, 0, 0], [0, 0, 1], [0, 1, 0]]]])

(1)

  f( point(A, 0, 0, 0), point(B, 0, 0, 0), point(C, 0, 0, 1), point(F, 0, 1, 0), t2);

no_tetrahedron_possible

(2)

 


Download istet.mw

the operation of the inner loop for the first two values of the index 'i'

when i=0, then b=0, inner loop executes once, because the 'until b=0' is only evaluated at the end of the loop - at which point the inner loop termantes

when i=1,then b=1, inner loop executes; c=1 (which is <>2), so the 'next' statement will cause this iteration of the inner loop to "exit", but still with b=1 - and this process repeats indefinitely!

Now obviously this can be fixed, but wouldn't it be simpler just to use something like the attached

  restart;
  f:=m-> seq
         ( `if`
           ( not member(2, convert(k, base, 3)),
             k,
             NULL
           ),
           k=1..m
         );
  f(20);
  f(50);

proc (m) local k; options operator, arrow; seq(`if`(not member(2, convert(k, base, 3)), k, NULL), k = 1 .. m) end proc

 

1, 3, 4, 9, 10, 12, 13

 

1, 3, 4, 9, 10, 12, 13, 27, 28, 30, 31, 36, 37, 39, 40

(1)

 

Download base3.mw

 

because no combination of strictly positive integers for x1, x2, x3, x4, x5, x6, x7, x8 will satisfy the constraint

10*x1 + 9*x2 + 15*x3 + 3*x4 + 11*x5 + 6*x6 + 3*x7 + 4*x8 <= 22

On the other hand, if you are prepared to accept non-negative integers, then you could use the attached

restart; with(Optimization)

LPSolve(5*x1+2*x2+7*x3+6*x4+x5+6*x6+8*x7+6*x8, {10*x1+9*x2+15*x3+3*x4+11*x5+6*x6+3*x7+4*x8 <= 22}, assume = {integer, nonnegative}, maximize = true)

[56, [x1 = 0, x2 = 0, x3 = 0, x4 = 0, x5 = 0, x6 = 0, x7 = 7, x8 = 0]]

(1)

NULL

Download intOpt.mw

 

 

See if you can work out why the 'toy' example in the attachedd worksheet returns the answer it does!

When you do understand it you will be bale to apply the same thought process to yur original worksheet

By the way it is inconsiderate to delete previous versions of your problem, particularly after they have been answered. If you continue to do this, then fewer and fewer people will biother to answer any new question you pos

limit(1.0/x, x = 0)

Float(undefined)

(1)

 

Download toyLim.mw

coeeff() only *really* works with polynomials.

Try the DEtools:-convertAlg() command as shown in the attached

expr:=(diff(P(t),t)+sin(Pi*x/L)*a[1](t))/(L^4*k*A*G)+sin(Pi*x/L)*diff(a[1](t),t$4)+sin(Pi*x/L)*diff(P(t),t);
DEtools:-convertAlg(expr,a[1](t))[1,1];

(diff(P(t), t)+sin(Pi*x/L)*a[1](t))/(L^4*k*A*G)+sin(Pi*x/L)*(diff(diff(diff(diff(a[1](t), t), t), t), t))+sin(Pi*x/L)*(diff(P(t), t))

 

sin(Pi*x/L)/(L^4*k*A*G)

(1)

coeff(expr, a[1](t));

Error, unable to compute coeff

 

 


 

Download getC.mw

 

In Maple by default,  the name 'gamma' is an "initially known name" - it is Euler's constant, or 0.5772156649 (approximately). You can overrride this definition if you want to by inserting the command

local gamma;

at the top of your worksheet. Alternatively, just use a different name!

###########################################################

You always have four equations with three unknowns.

The simplest way to illustrate why you sometimes get no solutions, is to 'solve' the system [f[1]=0, f[2]=0, f[3]=0] for [alpha, beta, g], as illustrated in the attached.

In case(1) the solution of the reduced system  is alpha=0, beta=0, g=g: in other words 'g' is arbitrary, so the additional condition alpha^2+beta^2+g^2=1 can be satisfied with g=+1 or -1

In cases(2..3) the solution of the reduced system  is alpha=0, beta=0, g=0: so it is impossible to satisfy the additional condition alpha^2+beta^2+g^2=1


 

restart

````

tau := Matrix(3, 3, {(1, 1) = 100, (1, 2) = 200, (1, 3) = 0, (2, 1) = 200, (2, 2) = -500, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 400})

 

q := tau-lambda*(Matrix(3, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1})) = Matrix(%id = 18446744074369822102)"(->)"-lambda^3+250000*lambda-36000000"(->)"-lambda^3+250000*lambda-36000000 = 0"(->)"[[lambda = 400], [lambda = -200-100*13^(1/2)], [lambda = -200+100*13^(1/2)]]"(->)"cst

 

NULL

lambda[I] := map[2](eval, lambda, cst)[1]

lambda[II] := map[2](eval, lambda, cst)[2]

lambda[III] := map[2](eval, lambda, cst)[3]

NULL

p := sort([evalf(lambda[I]), evalf(lambda[II]), evalf(lambda[III])])

 

sigma[I] := p[3]

sigma[II] := p[2]

sigma[III] := p[1]

NULL

NULL

NULLNULL

1:

eqP := alpha^2+beta^2+g^2 = 1

NULL

lambda := sigma[I]

f := q.(Vector(3, {(1) = alpha, (2) = beta, (3) = g}))

Vector[column](%id = 18446744074407773238)

(1)

NULL

sol1 := solve({f[1] = 0, f[2] = 0, f[3] = 0}, {alpha, beta, g}); solve(eval(eqP, sol1), g)

1, -1

(2)

NULL

``

``

2:

``

lambda := sigma[II]

f := q.(Vector(3, {(1) = alpha, (2) = beta, (3) = g}))

Vector[column](%id = 18446744074407752878)

(3)

solve({f[1] = 0, f[2] = 0, f[3] = 0}, {alpha, beta, g})

{alpha = 0., beta = 0., g = 0.}

(4)

``

``

3:

NULL

lambda := sigma[III]

f := q.(Vector(3, {(1) = alpha, (2) = beta, (3) = g}))

Vector[column](%id = 18446744074407746246)

(5)

NULL

solve({f[1] = 0, f[2] = 0, f[3] = 0}, {alpha, beta, g})

{alpha = 0., beta = 0., g = 0.}

(6)

``


 

Download noSol.mw

just because I prefer to use 1-D input for everything

  with(LinearAlgebra):
  Matrix
  ( 3,
    3,
    [ [  200,    0,   0],
      [  -50, -800,   0],
      [    0,    0, 200]
    ]
  ):
  %-lambda*Matrix
           ( 3,
             3,
             [ [1, 0, 0],
               [0, 1, 0],
               [0, 1, 1]
             ]
           ):
  Determinant(%):
  sigma:=[solve(%)];
  sigma[1];
  sigma[2];
  sigma[3];

[-800, 200, 200]

 

-800

 

200

 

200

(1)

  restart:
  with(LinearAlgebra):
  sigma:= [ solve
            ( Determinant
              ( Matrix
                ( 3,
                  3,
                  [ [  200,    0,   0],
                    [  -50, -800,   0],
                    [    0,    0, 200]
                  ]
                )-lambda*Matrix
                         ( 3,
                           3,
                           [ [ 1, 0, 0 ],
                             [ 0, 1, 0 ],
                             [ 0, 1, 1 ]
                           ]
                         )
              )
            )
          ];
  sigma[1];
  sigma[2];
  sigma[3];

[-800, 200, 200]

 

-800

 

200

 

200

(2)

 

NULL

Download use1D.mw

 

but a quick skim of the Wikipedia page suggests that it always finds roots of a polynomial in pairs.

So the process becomes

  1. Run  Bairstow'w method to return two roots (R1, R2) of a polynomial P1
  2. Divide the polynomial  P1 by the factors (x-R1)*(x-R2) to obtain a new polynomial P2.
  3. Run Bairstow's method again to return another two roots (R3 and R4) of the polynomial P2
  4. Divide the polynomial P2 by the factors (x-R3)*(x-R4) to obtain a new polynomial P3.
  5. Repeat as often as necessary

Based on your worksheet, I have (crudely) implemented this approach in the attached - and it appears to work!!!

restart; Bairstow := proc (n, a, u0, v0, itmax, TOL) local u, v, b, c, j, k, DetJ, du, dv, s1, s2; u := u0; v := v0; b := Array(0 .. n); c := Array(0 .. n); b[n] := a[n]; c[n] := 0; c[n-1] := a[n]; for j to itmax do b[n-1] := a[n-1]+u*b[n]; for k from n-2 by -1 to 0 do b[k] := a[k]+u*b[k+1]+v*b[k+2]; c[k] := b[k+1]+u*c[k+1]+v*c[k+2] end do; DetJ := c[0]*c[2]+(-1)*c[1]*c[1]; du := (-b[0]*c[2]+b[1]*c[1])/DetJ; dv := (b[0]*c[1]-b[1]*c[0])/DetJ; u := u+du; v := v+dv; printf("%3d %12.7f %12.7f %12.4g %12.4g\n", j, u, v, du, dv); if max(abs(du), abs(dv)) < TOL then break end if end do; printf("\nQ(x)=(%g)x^%d", b[n], n-2); for k from n-3 by -1 to 1 do printf(" + (%g)x^%d", b[k+2], k) end do; printf(" + (%g)\n", b[2]); printf("Remainder: %g(x-(%g))+(%g)\n", b[1], u, b[0]); printf("Quadratic Factor: x^2-(%g)x-(%g)\n", u, v); s1 := evalf((1/2)*u+(1/2)*sqrt(u*u+4*v)); if u^2+4*v < 0 then printf(" Zeros: %.13g +-(%.13g)i\n", Re(s1), abs(Im(s1))) else s2 := evalf((1/2)*u-(1/2)*sqrt(u*u+4*v)); printf(" Zeros: %.13g, %.13g\n", s1, s2) end if end proc; P := proc (x) options operator, arrow; x^5+(-1)*3.5*x^4+2.75*x^3+2.125*x^2+(-1)*3.875*x+1.25 end proc; n := degree(P(x)); a := Array(0 .. n); a[0] := P(0); for i to n do a[i] := coeff(P(x), x^i) end do; itmax := 10; TOL := 10^(-10); r := -1; s := -1; Bairstow(n, a, r, s, itmax, TOL)

  1   -0.6441699    0.1381090       0.3558        1.138
  2   -0.5111131    0.4697336       0.1331       0.3316
  3   -0.4996865    0.5002023      0.01143      0.03047
  4   -0.5000001    0.5000000   -0.0003136   -0.0002023
  5   -0.5000000    0.5000000    6.413e-08    9.268e-09
  6   -0.5000000    0.5000000            0            0

Q(x)=(1)x^3 + (-4)x^2 + (5.25)x^1 + (-2.5)
Remainder: 0(x-(-0.5))+(0)
Quadratic Factor: x^2-(-0.5)x-(0.5)
 Zeros: 0.4999999993, -0.9999999997

 

P := unapply(convert(simplify(evalc(P(x)/((x-1/2)*(x+1)))), rational), x); n := degree(P(x)); a := Array(0 .. n); a[0] := P(0); for i to n do a[i] := coeff(P(x), x^i) end do; itmax := 10; TOL := 10^(-10); r := -1; s := -1; Bairstow(n, a, r, s, itmax, TOL)

proc (x) options operator, arrow; -5/2+x^3-4*x^2+(21/4)*x end proc

 

  1    1.2413793    3.1982759        2.241        4.198
  2    0.6131527   -2.7786768      -0.6282       -5.977
  3    1.3134052   -1.2310587       0.7003        1.548
  4    1.7930398   -1.0627832       0.4796       0.1683
  5    1.9953201   -1.2091046       0.2023      -0.1463
  6    2.0002856   -1.2499754     0.004965     -0.04087
  7    1.9999999   -1.2499999   -0.0002857   -2.449e-05
  8    2.0000000   -1.2500000    1.419e-07   -8.164e-08
  9    2.0000000   -1.2500000    1.367e-14   -2.012e-14

Q(x)=(1)x^1 + (-2)
Remainder: 2.01221e-14(x-(2))+(1.70836e-14)
Quadratic Factor: x^2-(2)x-(-1.25)
 Zeros: 1 +-(0.5)i

 

P := unapply(convert(simplify(evalc(P(x)/((x-1+I*(1/2))*(x-1-I*(1/2))))), rational), x)

proc (x) options operator, arrow; x-2 end proc

(1)

NULL

NULL

NULL


 

Download Bairs.mw

 

the attached, maybe?


 

for i from 1 to 20 do:
    printf("new variable %2d\n", i)
end do;

new variable  1
new variable  2
new variable  3
new variable  4
new variable  5
new variable  6
new variable  7
new variable  8
new variable  9
new variable 10
new variable 11
new variable 12
new variable 13
new variable 14
new variable 15
new variable 16
new variable 17
new variable 18
new variable 19
new variable 20

 

 


 

Download prStuff.mw

A caveat - I'm using 64-bit Windows 7

Maple 18     Created file size 56kB.      only thing GIMP renders is the text associated with axes tickmarks.

Maple 2015 Created file size 8.6MB.    Plot looks fine in GIMP

Maple 2016 Created file size 7.3MB.    Plot looks fine in GIMP

Maple 2017 Created file size 7.3MB.    Plot looks fine in GIMP

Maple 2018 Created file size 7.3MB.    Plot looks fine in GIMP

Maple 2019 Created file size 0!.           Maple "hangs", has to be terminated with CTRL+ALT+Delete

Maple 2020 Created file size 0!.           Maple "hangs", has to be terminated with CTRL+ALT+Delete

 

but maybe something like the attached?


 

  restart;
  with(plots):
  with(plottools):

  K:= 48: R1:= 30: R2:= 20: OA:= R1+R2:
  Vodilo:= display( cuboid( [-3, -3, 0],
                            [OA+4, 4, 2]
                          )
                  ):
  Koleso1:= display( cylinder([0, 0, 5], R1, 18),
                     cylinder([0, 0, -2], 1, 16)
                   ):
  Koleso2:= display( cylinder([0, 0, 7], R2, 6),
                     cylinder([0, 0, -2], 1, 16)
                   ):
  for i from 0 to K do
      Amplituda:= .3;
      phi:= sin(6.28*i/K)*Amplituda;
      phi1:= phi*(R1+R2)/R2;
      P1:= rotate(Koleso2, 0, 0, -phi1);
      P2:= rotate(Vodilo, 0, 0, -phi);
      x:= OA*cos(phi);
      y:= OA*sin(phi);
      P[i]:= display( Koleso1,
                      translate(P1, x, y, 0),
                      P2,
                      translate(P1, -x, -y, 0)
                    );
  end do:
  display( seq(P[i], i = 0 .. K),
           orientation = [55, 101],
           insequence = true,
           scaling = constrained,
           axes = normal
         );

 

 


 

Download animplt2.mw

what you were looking for?


 

  restart;
  with(plots):
  with(plottools):

  K:= 48: R1:= 30: R2:= 20: OA:= R1+R2:
  Vodilo:= display( cuboid( [-3, -3, 0],
                            [OA+4, 4, 2]
                          )
                  ):
  Koleso1:= display( cylinder([0, 0, 5], R1, 18),
                     cylinder([0, 0, -2], 1, 16)
                   ):
  Koleso2:= display( cylinder([0, 0, 7], R2, 6),
                     cylinder([0, 0, -2], 1, 16)
                   ):
  for i from 0 to K do
      Amplituda:= .3;
      phi:= sin(6.28*i/K)*Amplituda;
      phi1:= phi*(R1+R2)/R2;
      P1:= rotate(Koleso2, 0, 0, -phi1);
      P2:= rotate(Vodilo, 0, 0, -phi);
      x:= OA*cos(phi);
      y:= OA*sin(phi);
      P[i]:= display( Koleso1,
                      translate(P1, x, y, 0),
                      P2
                    );
  end do:
  display( seq(P[i], i = 0 .. K),
           orientation = [55, 101],
           insequence = true,
           scaling = constrained,
           axes = normal
         );

 

 


 

Download animplt.mw

4 5 6 7 8 9 10 Last Page 6 of 132