Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 182 days
Contact tomleslie

MaplePrimes Activity


These are replies submitted by tomleslie

Am I the only one who has noticed?...

tomleslie 13876
April 18 2018
0 0

that the previous suggestion starts the x-axis at 0.5 rather than 0.0, Surely you want the latter, as in

restart;
DynamicSystems:-DiscretePlot( [1,2,3,4,5],
                                                        [1,6,2,3,4],
                                                        style=stem,
                                                        view=[0..5.5, default]
                                                      )

Actually I've just noticed...

tomleslie 13876
April 17 2018
0 0

@Carl Love 

This code only stops if the current sample is the same as the previous one. But for some start numbers, one gets a recurring cycle of length greater than 1. Try yours and mine for a seed of 12345

The concrete example...

tomleslie 13876
April 17 2018
0 0

The thing you have to remember is that the cardioid intercepts the y-axis in two places, one negative and one positive. The positive y-intercept is the only valid solution

The attached worksheet contains four execution groups

  1. The first one just sets up the problem
  2. The second group uses the two-stage process using fsolve(). There is no way I know to control whether the first-stage fsolve() returns the valid, or the invalid solution: so I check it. If it has returned the invalid solution, then I just re-execute with the "avoid" option set for the invalid solution. fsolve() then returns the valid solution. This way I get all the valid solutions for 33 separate angles in ~1sec. So let's call this 30mSec/angle
  3. The third execution group uses DirectSearch() with limited constraints and the option AllSolutions=true, solutions=2. This returns both solutions for one angle in 1.8sec. So about 60X slower than the fsolve() process in (2) above - and I still have the overhead of figuring out which of the two returned solution is valid. This overhead would be minimal, so I'm not even going to count it.
  4. The fourth execution group uses DirectSearch() with AllSolutions=true, solutions=1, but with the added constraint that the y-intercept be positive. This correctly returns the single valid solution for one angle , but takes ~48secs. In other words it is ~1500X slower per angle, than the two stage fsolve() process in (2) above
> 

  restart;
  with(plots):
  with(plottools):
#
# Generate the cardioid, and set up a rotation
# just to check verything is working properly
#
  spX:=  piecewise
         ( theta < 0, undefined,
           theta <= Pi,   cos(theta+(1/4)*Pi)*exp((theta+(1/4)*Pi)*(1/2)),
           theta <= 2*Pi, sin(theta-3*Pi/4)*exp((theta-3*Pi/4)*(1/2)),
           theta > 2*Pi, undefined
         ):
  spY:=  piecewise
         ( theta < 0, undefined,
           theta <= Pi,   sin(theta+(1/4)*Pi)*exp((theta+(1/4)*Pi)*(1/2)),
           theta <  2*Pi, cos(theta-3*Pi/4)*exp((theta-3*Pi/4)*(1/2)),
           theta > 2*Pi, undefined
         ):
  p1:=plot( [ spX, spY, theta=0..2*Pi ],
              thickness=3,
              color=red,
              discont=true
          ):
#
# Define a rotation matrix which allows the above
# plot to be rotated by any specified angle beta,
# and produce the "rotated" x and y-coordinates
#
  mRot:=Matrix(2,2, [ [ cos(beta), -sin(beta)],
                      [ sin(beta),  cos(beta)]
                    ]
              ):
  rS:= mRot.< spX, spY >:
#
# Set the number of solutions to be used for all
# calculations.
#
  numSols:= 0..32:
> 

t1:= time():
#
# Attempt the same calculation using fsolve(). The problem
# with using fsolve() is that there are always two answers
# for where the cardioid crosses the y-axis, one positive
# and one negative. As far as I can tell there is no way
# to control which is returned:-(
#
# The only way I could come up with to produce only the
# the positive y-axis intercept, was a two stage process.
#
# For any give value of beta (the rotation angle)
#
#  1) Just run fsolve() to compute the angle theta for
#     which the x-coordinate is zero, without worrying
#     whether the resulting y-axis intercept is positive
#     or negative
#  2) Evaluate the y-axis intercept at the value of theta
#     obtained in (1). If the intercept is negative, then
#     re-un the fsolve() command with the 'avoid' option.
#     This results in the "other" value for theta, and
#     hence gives a positive y-intercept
##########################################################
#
# Get the first stage of solutions
#      
  ans:= seq
        ( fsolve
          ( eval
            ( rS[1],
              beta=k*Pi/(op(2, numSols)/2)
            ),
            theta=0..2*Pi
          ),
          k=numSols
        ):
#
# Check whether the resulting y-intercept is positive
# If it isn't rerun the fsolve() command with the 'avoid'
# option
#
  ans:= seq
        ( `if`
          ( evalf
            ( eval
              ( rS[2],
                [ beta=k*Pi/(op(2, numSols)/2),
                  theta=ans[k+1]
                ]
              )<0
            ),
            fsolve
            ( eval
              ( rS[1],
                beta=k*Pi/(op(2, numSols)/2)
              ),
              theta=0..2*Pi,
              avoid={ theta=ans[k+1] }
            ),
            ans[k+1]
         ),
         k=numSols
       ):
  plt:= seq
        ( display
          ( [ pointplot
              ( [ 0,
                  evalf
                  ( eval
                    ( rS[2],
                      [ beta=k*Pi/(op(2, numSols)/2),
                        theta=ans[k+1]
                      ]
                    )
                  )
                ],
                symbol=solidcircle,
                symbolsize=40,
                color=blue
              ),
              rotate( p1, k*Pi/(op(2, numSols)/2))
            ],
            scaling=constrained
          ),
          k=numSols
       ):
  display(plt, insequence=true);
t2:=time()-t1;

 

.920

 (1)
> 

#
# Get both solutions for a given angle. Since the
# constraint for a positive intercept is removed
# no-one knows which of these solutions is the
# desired one
#
  with(DirectSearch):
  k:=0:
  CodeTools:-Usage
             ( SolveEquations
               ( eval
                 ( rS[1],
                   beta=k*Pi/(op(2, numSols)/2)
                 ),
                 [ #eval
                   #( rS[2],
                   #  beta=k*Pi/(op(2, numSols)/2)
                   #) >0,
                   theta>=0,
                   theta<=2*Pi
                 ],
                 AllSolutions=true, solutions=2
               )
             );

memory used=148.61MiB, alloc change=16.00MiB, cpu time=1.79s, real time=1.77s, gc time=109.20ms

  

Matrix(2, 4, {(1, 1) = 0.7382467866e-24, (1, 2) = Vector(1, {(1) = 0.8592128878e-12}), (1, 3) = [theta = .7853981633970569], (1, 4) = 21, (2, 1) = 0.1695803329e-20, (2, 2) = Vector(1, {(1) = 0.4118013269e-10}), (2, 3) = [theta = 5.497787143773578], (2, 4) = 22})

 (2)
> 

#
# Re-insert the constraint that rS[2] be positive.
# Now require only one solution. Note how the
# execution time increases by 30X
#
  with(DirectSearch):
  k:=0:
  CodeTools:-Usage
             ( SolveEquations
               ( eval
                 ( rS[1],
                   beta=k*Pi/(op(2, numSols)/2)
                 ),
                 [ eval
                   ( rS[2],
                     beta=k*Pi/(op(2, numSols)/2)
                   ) >0,
                   theta>=0,
                   theta<=2*Pi
                 ],
                 AllSolutions=true, solutions=1
               )
             );

memory used=4.26GiB, alloc change=344.00MiB, cpu time=47.47s, real time=45.47s, gc time=3.81s

  

Vector[row](4, {(1) = 0.3031554666e-22, (2) = Vector(1, {(1) = -0.5505955563e-11}), (3) = [theta = .785398163399959], (4) = 21})

 (3)
> 

 


 

Download DSprob2.mw

I agree...

tomleslie 13876
April 17 2018
0 0

with vv

We've all had a lot of "fun" with this problem, but essentially all that anyone has achieved is to find the positive y-intercept of the rotating cardioid. The "disk" representing the cam follower has then been offset from this intercept by a fixed amount, equal to the radius of the "disk". This may "look" OK on the plots which you have seen, but it absolutely does not guarantee that the disk and the cardioid never overlap.

This can be most easily visualised by changing the "disk" to a "circle" and increasing its size as in the attached. Use the drag handles on the plot to increase its size as much as your screen will allow, and then use the toolbar command to "step" through the animation. You will observe that the "disk" and the "cardioid" overlap in nearly every frame.

camFoll2.mw

This is not an "error" in the calculations provided, just the manifestation of the relative radii of curvature of the "disk" and the "cardioid", which means that "tangency" is never guaranteed

Amongst other things, you might want to consider the profile of the cam. As defined, the "cusp" (ie the "dip" in the cardioid) is a discontinuous point. It doesn't matter how small you make the "follower", it will never reach the "bottom" of this "dip", so why make the cam like this???

Constraint removal causes problems...

tomleslie 13876
April 17 2018
0 0

For any angle, the OP's problem has two solutions, only one of which is valid.

Removing the constraint rS[2]>0, and using "AllSolutions=true, solutions=1" speeds things up a bit, but means that only one solution is returned - but it may be the wrong one!

Removing the constraint rS[2]>0, and using "AllSolutions=true", returns both solutions, and is even faster, but then further processing is required, in order to find the valid solution

So the constraint rS[2]>0 is necessary in order to avoid extra processing - and removing it is not " Improving the last piece of your code "

Aaaah - floor()...

tomleslie 13876
April 17 2018
0 0

I went down a rabbit-hole trying to use frem() for more or less the same purpose, and I never reappeared.

Anyhow, for what it's worth, another version

> 

  restart;
  with(plots):
  with(plottools):
#
# Define the "Archimedean" cardioid
#
  r1:= t -> piecewise
            ( t<Pi,   sqrt(exp((t+(1/4)*Pi)*(1/2))),
              t<2*Pi, sqrt(exp((9*Pi/4-t)*(1/2)))
            ):
  r:= t-> r1(t - 2*Pi*floor(t/(2*Pi))):
#
# Set the radius of the "ball" and length of
# the line attached to the "ball" for the follower
#
  dRad:=0.2:
  ll:=1.0:
#
# Combined plot for one "frame" of the animation
#
  F:= u-> display
          ( [ disk
              ( [ 0, dRad+r(Pi/2-u)*sin(Pi/2) ],
                dRad,
                color=blue
              ),
              line
              ( [ 0, r(Pi/2-u)*sin(Pi/2)+dRad ],
                [ 0, ll+2*dRad+r(Pi/2-u)*sin(Pi/2) ],
                color=blue,
                thickness =5
              ),
              plot
              ( [ r(t-u)*cos(t),
                  r(t-u)*sin(t),
                  t=0..2*Pi
                ],
                discont=true
              )
            ]
          ):
#
# Set number of frames and produce animation
#
  nframes:=64:
  display
  ( seq( F(p),
         p = 0..evalf(2*Pi), evalf(2*Pi/nframes)
       ),
    insequence=true,
    scaling=constrained,
    axes=none
  );

 
> 

 

Download camFoll.mw

"Archimedean Spiral" in polars...

tomleslie 13876
April 16 2018
0 0

@Earl 

I'm pretty sure that the following is what you need for your specific cardioid. However when I tried to patch this into vv's solution (which has got to be the way to go) something broke, and I couldn't immediately fix it. It's past my bedtime, I may get back to this tomorrow

> 

restart;
with(plots):
with(plottools):
r:=t->piecewise(t<Pi, sqrt(exp((t+(1/4)*Pi)*(1/2))),
                t<2*Pi, sqrt(exp((9*Pi/4-t)*(1/2)))
               ):
p1:= plot( [r(t)*cos(t), r(t)*sin(t), t=0..2*Pi], discont=true):
p2:= [seq(rotate(p1, u), u=0..evalf(2*Pi), evalf(Pi/32))]:
display(p2, insequence=true);

 
> 

 

Download archSP.mw

The code file you requested is attached...

tomleslie 13876
April 16 2018
0 0

@Markiyan Hirnyk 

  1. The first execution group just sets up the problem
  2. The second execution group uses the two-stage fsolve() approach, just so that one can see the "correct" solution. Takes about 1 second to compute the solution for 33 angles
  3. The third execution group attempts to solve the same problem using DirectSearch:-SolveEquations()  without the AllSolutions=true option, just to demonstrate that some incorrect solution are obtained. This takes about 1 second to execute for 33 angles
  4. The fourth execution group attempts to solve the same problem using DirectSearch:-SolveEquations()  with the AllSolutions=true option. Because I was unsure of the format of the output, I removed all post-processing of the results returned. I killed this execution group after about 5 minutes, because I wasn't sure if I was getting anywhere
  5. The fifth execution group uses DirectSearch:-SolveEquations()  with the AllSolutions=true option, but for only one angle ( rather than the 33 angles used in bullets 2-4 above). This took 10secs to return two solutions. Based on the residues, one of these is obviously not a "solution". Now I could select the desired solution, by searching for the minimum residue, but given the time taken (~10secs/angle, so ~330secs for 33 angles?), I'll stick with the approach in (2) above which does the whole thing in 1 sec

Code for this exercise is attached below

> 

  restart;
  with(plots):
  with(plottools):
#
# Generate the cardioid, and set up a rotation
# just to check verything is working properly
#
  spX:=  piecewise
         ( theta < 0, undefined,
           theta <= Pi,   cos(theta+(1/4)*Pi)*exp((theta+(1/4)*Pi)*(1/2)),
           theta <= 2*Pi, sin(theta-3*Pi/4)*exp((theta-3*Pi/4)*(1/2)),
           theta > 2*Pi, undefined
         ):
  spY:=  piecewise
         ( theta < 0, undefined,
           theta <= Pi,   sin(theta+(1/4)*Pi)*exp((theta+(1/4)*Pi)*(1/2)),
           theta <  2*Pi, cos(theta-3*Pi/4)*exp((theta-3*Pi/4)*(1/2)),
           theta > 2*Pi, undefined
         ):
  p1:=plot( [ spX, spY, theta=0..2*Pi ],
              thickness=3,
              color=red,
              discont=true
          ):
#
# Define a rotation matrix which allows the above
# plot to be rotated by any specified angle beta,
# and produce the "rotated" x and y-coordinates
#
  mRot:=Matrix(2,2, [ [ cos(beta), -sin(beta)],
                      [ sin(beta),  cos(beta)]
                    ]
              ):
  rS:= mRot.< spX, spY >:
#
# Set the number of solutions to be used for all
# calculations.
#
  numSols:= 0..32:
> 

#
# Attempt the same calculation using fsolve(). The problem
# with using fsolve() is that there are always two answers
# for where the cardioid crosses the y-axis, one positive
# and one negative. As far as I can tell there is no way
# to control which is returned:-(
#
# The only way I could come up with to produce only the
# the positive y-axis intercept, was a two stage process.
#
# For any give value of beta (the rotation angle)
#
#  1) Just run fsolve() to compute the angle theta for
#     which the x-coordinate is zero, without worrying
#     whether the resulting y-axis intercept is positive
#     or negative
#  2) Evaluate the y-axis intercept at the value of theta
#     obtained in (1). If the intercept is negative, then
#     re-un the fsolve() command with the 'avoid' option.
#     This results in the "other" value for theta, and
#     hence gives a positive y-intercept
##########################################################
#
# Get the first stage of solutions
#      
  ans:= seq
        ( fsolve
          ( eval
            ( rS[1],
              beta=k*Pi/(op(2, numSols)/2)
            ),
            theta=0..2*Pi
          ),
          k=numSols
        ):
#
# Check whether the resulting y-intercept is positive
# If it isn't rerun the fsolve() command with the 'avoid'
# option
#
  ans:= seq
        ( `if`
          ( evalf
            ( eval
              ( rS[2],
                [ beta=k*Pi/(op(2, numSols)/2),
                  theta=ans[k+1]
                ]
              )<0
            ),
            fsolve
            ( eval
              ( rS[1],
                beta=k*Pi/(op(2, numSols)/2)
              ),
              theta=0..2*Pi,
              avoid={ theta=ans[k+1] }
            ),
            ans[k+1]
         ),
         k=numSols
       ):
  plt:= seq
        ( display
          ( [ pointplot
              ( [ 0,
                  evalf
                  ( eval
                    ( rS[2],
                      [ beta=k*Pi/(op(2, numSols)/2),
                        theta=ans[k+1]
                      ]
                    )
                  )
                ],
                symbol=solidcircle,
                symbolsize=40,
                color=blue
              ),
              rotate( p1, k*Pi/(op(2, numSols)/2))
            ],
            scaling=constrained
          ),
          k=numSols
       ):
  display(plt, insequence=true);

 
> 

#
# Using DirectSearch in a braindead manner - note
# that some (but not all)of the solutions are
# incorrect
#
  with(DirectSearch):
  ans:= [ seq
          ( rhs( SolveEquations
                 ( eval
                   ( rS[1],
                     beta=k*Pi/(op(2, numSols)/2)
                   ),
                   [ eval
                     ( rS[2],
                       beta=k*Pi/(op(2, numSols)/2)
                     ) >0,
                     theta>=0,
                     theta<=2*Pi
                   ],
                   usewarning=false
                 )[3][]),
            k=numSols
          )
        ]:
  plt:=seq( display
            ( [ pointplot
                ( [ 0,
                    evalf
                    ( eval
                      ( rS[2],
                        [beta=k*Pi/8, theta=ans[k+1] ]
                      )
                    )
                  ],
                  symbol=solidcircle,
                  symbolsize=40,
                  color=blue
                ),
                rotate( p1, k*Pi/(op(2, numSols)/2))
              ],
              scaling=constrained
            ),
            k=numSols
         ):
#
# NB The following animationa is OK for some values,
# but wrong for others
#
  display(plt, insequence=true);

 
> 

#
# Using DirectSearch with the AllSolution=true option
# Don't know what the output format will be so just
# generate the solutions, and worry about what to do
# with them later
#
# Unforunately despite waiting for about 10 mins, no
# output appeared - see next execution group
#
  ans:= [ seq
          ( SolveEquations
            ( eval
              ( rS[1],
                beta=k*Pi/(op(2, numSols)/2)
              ),
              [ eval
                ( rS[2],
                  beta=k*Pi/(op(2, numSols)/2)
                ) >0,
                theta>=0,
                theta<=2*Pi
               ],
               AllSolutions=true
            ),
            k=numSols
          )
        ]:

Warning,  computation interrupted

  
> 

#
# using DirecSearch() with AllSolotions=true again, but
# instead of trying to find the sequence of solutions
# which took waaaaay too long, just find the first
# solution
#
  k:=0:
  CodeTools:-Usage
             ( SolveEquations
               ( eval
                 ( rS[1],
                   beta=k*Pi/(op(2, numSols)/2)
                 ),
                 [ eval
                   ( rS[2],
                     beta=k*Pi/(op(2, numSols)/2)
                   ) >0,
                   theta>=0,
                   theta<=2*Pi
                 ],
                 AllSolutions=true
               )
             );

memory used=4.21GiB, alloc change=32.00MiB, cpu time=55.27s, real time=49.82s, gc time=10.02s

  

Matrix(2, 4, {(1, 1) = 0.2631193646e-24, (1, 2) = Vector(1, {(1) = 0.5129516202e-12}), (1, 3) = [theta = .7853981633972148], (1, 4) = 27, (2, 1) = 1.096640025369003, (2, 2) = Vector(1, {(1) = 1.0472058180553634}), (2, 3) = [theta = 3.1415926535897905], (2, 4) = 55})

 (1)
> 

 

Download DSprob.mw

Don't use sum()...

tomleslie 13876
April 16 2018
0 0

@tsunamiBTP 

use add() or Sum() - the trouble with sum() is that it spends ages trying to find a closed-form solution before realising that it can't. So having spent ages getting nowhere, it tries add() abd if that doesn't work it just returns sum() ie unevaluated

If you want to avoid the long-winded and pointless process of trying sum() where it will never work, then stick to add() or Sum()

Check out Acer's solution...

tomleslie 13876
April 16 2018
0 0

@tsunamiBTP 

I tried a variety of methods, which can be characterised as

evalf + add()
evalf + Sum()
evalhf + add() <- suggested by Acer, with whom one should never argue
evalhf + Sum()

Results based on my machine


                                 Memory      Cpu
 evalf()+add()        121.51MiB   983.00ms
 evalf()+Sum()        145.84MiB     1.39s
 evalhf()+add()          3.74KiB    15.00ms
 evalhf()+Sum()          2.57KiB    16.00ms

 Conclusions

  1. Difference between add() and Sum() small:
    1. Using software floats add() is slightly  faster than Sum() and uses slightly less memory
    2. Using hardware floats add() is slightly faster than Sum() but uses slightly more memory
  2. Difference between evalhf() and evalf() - huge
    1. For add(), memory footprint is improved by about 30000X, speed improvent is about 70X
    2. For Sum, memory footprint is improved by 55000X, speed improvement is about 85X

Hence - definitely use hardware floats

I rewote the worksheet in my original response to perform the above comparisons, and then compute all of the desired cases using hardware floats (as per Acer's suggestion)

> 

restart;
S4 :=x->(2*((2*k*Pi*tau*(-8*Pi^2*k^2*tau^3+4*Pi^2*k^2*tau^2*w+2*T^2*tau+T^2*w)*cos(w*Pi*k/T)+T*sin(w*Pi*k/T)*(-16*Pi^2*k^2*tau^3+4*Pi^2*k^2*tau^2*w+T^2*w))*exp(-(1/2)*w/tau)+16*Pi^3*k^3*tau^4-4*k*Pi*tau^2*T^2))*sin(2*k*Pi*x/T)/(4*Pi^2*k^2*tau^2+T^2)^2:
> 

a[0] := 0: Kappa := 4: tau := 1: N := 2:
T := M*tau: w := N*tau: M := Kappa*N:
> 

S5add:=unapply(a[0]+'add'(S4(x), k = 1 .. m), [m,x]):
S5Sum:=unapply(a[0]+Sum(S4(x), k = 1 .. m), [m,x]):
> 

#
# test combinations of evalf(), evalhf(), add() and sum()
# Conclusion - use hardware floats; doesn't seem to make
# much difference whether one uses Sum() or add()
#
  forget(sin): forget(cos): forget(evalf): forget(evalhf):
  peak:= CodeTools:-Usage( evalf(S5add(10000, 0.3998201722e-3)));
  (peak-1)*(1/2);
  forget(sin): forget(cos): forget(evalf): forget(evalhf):
  peak:= CodeTools:-Usage( evalf(S5Sum(10000, 0.3998201722e-3)));
  (peak-1)*(1/2);
  forget(sin): forget(cos): forget(evalf): forget(evalhf):
  peak:= CodeTools:-Usage( evalhf(S5add(10000, 0.3998201722e-3)));
  (peak-1)*(1/2);
  forget(sin): forget(cos): forget(evalf): forget(evalhf):
  peak:= CodeTools:-Usage( evalhf(S5Sum(10000, 0.3998201722e-3)));
  (peak-1)*(1/2);

memory used=121.51MiB, alloc change=74.01MiB, cpu time=999.00ms, real time=994.00ms, gc time=46.80ms

  

1.178180172

  

0.890900860e-1

  

memory used=146.00MiB, alloc change=-2.00MiB, cpu time=1.31s, real time=1.30s, gc time=31.20ms

  

1.178180179

  

0.890900895e-1

  

memory used=3.74KiB, alloc change=0 bytes, cpu time=16.00ms, real time=13.00ms, gc time=0ns

  

1.17818017931755792

  

0.890900895e-1

  

memory used=2.57KiB, alloc change=0 bytes, cpu time=15.00ms, real time=15.00ms, gc time=0ns

  

1.17818017931755792

  

0.890900895e-1

 (1)
> 

#
# Compute all of OP's cases using add() with hardware floats
#
peak:= CodeTools:-Usage( evalhf(S5add(10000, 0.3998201722e-3)));
(peak-1)*(1/2);
peak:= CodeTools:-Usage( evalhf(S5add(50000, 0.7999280138e-4)));
(peak-1)*(1/2);
peak := CodeTools:-Usage( evalhf(S5add(100000, 0.3998201722e-4)));
(peak-1)*(1/2);
peak := CodeTools:-Usage( evalhf(S5add(500000, 0.7999280138e-5)));
(peak-1)*(1/2);
peak := CodeTools:-Usage( evalhf(S5add(1000000, 0.3998201722e-5)));
(peak-1)*(1/2);
peak := CodeTools:-Usage( evalhf(S5add(5000000, 0.7999280138e-6)));
(peak-1)*(1/2);
peak := CodeTools:-Usage( evalhf(S5add(10000000, 0.3998201722e-6)));
(peak-1)*(1/2);
 

memory used=2.57KiB, alloc change=0 bytes, cpu time=16.00ms, real time=14.00ms, gc time=0ns

  

1.17818017931755792

  

0.890900895e-1

  

memory used=2.59KiB, alloc change=0 bytes, cpu time=62.00ms, real time=68.00ms, gc time=0ns

  

1.17881976187587423

  

0.894098810e-1

  

memory used=2.59KiB, alloc change=0 bytes, cpu time=140.00ms, real time=134.00ms, gc time=0ns

  

1.17889958508458759

  

0.894497925e-1

  

memory used=2.59KiB, alloc change=0 bytes, cpu time=687.00ms, real time=691.00ms, gc time=0ns

  

1.17896373808754462

  

0.894818690e-1

  

memory used=2.59KiB, alloc change=0 bytes, cpu time=1.39s, real time=1.39s, gc time=0ns

  

1.17897154636834256

  

0.894857730e-1

  

memory used=2.59KiB, alloc change=0 bytes, cpu time=6.94s, real time=6.93s, gc time=0ns

  

1.17897813653769612

  

0.894890685e-1

  

memory used=2.59KiB, alloc change=0 bytes, cpu time=13.81s, real time=13.81s, gc time=0ns

  

1.17897874270361847

  

0.894893715e-1

 (2)
> 

 

NULL


 

Download someSums2.mw

Close, but still no coconut...

tomleslie 13876
April 15 2018
0 0

Suggest you supply maple code for your equations, because if you think I am going to retype them in Maple - you are wrong.

Setting Q(x,t)=0 will help because this will simplify your equations

You supply boundary conditions for q(x,t), but you will also need boundary conditions for T(x,t)

If you have any questions...

tomleslie 13876
April 15 2018
0 0

@Earl 

Then post here

The worksheet I eventually posted is unnecessarily complicated because

  1. I left in a complete execution group dependent on DirectSearch(), and basically this doesn't work for all angles. You can play with various DirectSearch options here; you can swap from SolveEquations() to Search() to GlobalSearch() etc with appropriate minor reformulations of the objective function. You can supply different options to any/all of these commands: believe me I tried, and I couldn't come up with anything which "worked" for all angles - and I really don't understand why. Prior to your question I had always been very impressed with the capabilities whic DirectSearch() commands provide. But basically, this execution group doesn't provide what you want, so can be deleted
  2. The final two execution groups using the fsolve() approach are pretty much duplicates. The only real difference, is that the first just draws a "blob" at the positive y-intercept of the cardioid. Having correctly obtained the position of this y-intercept, it is relatively simple (using plot commands) to construct a "cam-follower" rather than "blob" at this y-intercept. So you only really need one of the two final execution groups - probably the second

 

tried many option...

tomleslie 13876
April 15 2018
0 0

@Markiyan Hirnyk 

I tried many different approaches using DirectSearch() - amongst these were

  1. Attempting to supply useful initial points
  2. Tightening solution tolerances (as well as increasing DIgits settings)
  3. Instead of using SolveEquations( f(beta, theta)= 0 ), I also tried using Search() and GlobalSearch() with minimizing (f(beta, theta)^2
  4. In all cases, I could get the "correct" solutions for many/most of the required angles (ie beta-values) - but for some angles (ie beta-values) , I would get obviously invalid solutions, and I have no idea why:-(
  5. Having spent a couple of hours trying to persuade DirectSearch() to provide the required correct solutions, basically I failed miserably.
  6. If I had left all the DirectSolutions() variants I tried, in the worksheet I eventually posted, the latter would have been about ten times as long, and filled with attempts which didn't work. Posting all the stuff which doesn't work seems a little pointless to me. I left in one execution group dependent on DirectSearch(), mainly so that others could try it with modifications if necessary (and also maybe if I wanted to come back to it???).
  7. Feel free to modify the DirectSearch() execution group any way you want, in order to come up with the correct solutions - just don't expect me to do it, because
  8. Having achieved the "correct" solutions using the two-stage fsolve() calculation, there seems little point
  9. If you wish to claim that DirectSearch() can provide all of the required solutions, fine then do it: because I tried and failed, and I'm not planning on trying any more

Hmmm, granted...

tomleslie 13876
April 12 2018
0 0

But which of the following  is "better"?

  1. solve the ode numerically
  2. solve the ode "exactly" knowing that the only way to evaluate the "exact" solution for a given argument is numerically

The code I posted does it both ways - and they disagree in the last couple of decimal places

Hmmm...

tomleslie 13876
April 12 2018
0 0

I've tested my code in Maple 2016.1 on windows7 x64, and it works perfectly.

One other thing which might be worth checking. Enter

?Fractals,LSystem,LSystemExamples;

at the Maple prompt. This should bring up the help window. In the help window toolbar, use

View->Open Page as Worksheet

This should open a new worksheet in your Maple session. Execute this worksheet using the !!! toolbar button.

If everything works, this suggests (but doess not prove!) that your problem is related to the specific worksheet you are using to generate the dragon curve

If you get any errors, this suggests (but doess not prove!) that your problem is more general - reinstalling Maple may be worth a try

First 99 100 101 102 103 104 105 Last Page 101 of 207