Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 181 days
Contact tomleslie

MaplePrimes Activity


These are replies submitted by tomleslie

As many values of x[30] as you want...

tomleslie 13876
October 28 2018
0 0

@Zeineb 

An issue with extrapolating an "autoregressive" fit is that it depends on the "start-value" you select, for which there are 13 possibiliies - any of the values in your original 'data' vector. The attched shows all possibilities.

In addition it uses the Statistics package to produce linear and quadratic fits of data against time, which can also be used for extrapolation. You could probably get almost any number you want - whihc at least demonstrates the dangers of extrapolation!!! ("Lies, damned lies and statistics", as a wise man once said - and we don't even know who!!)

> 

  restart:
  with(plots):
  data := <.150, .175, .150, .175, .200, .300, .250, .250, .350, .400, .450, .550, .550>:
  fit:= minimize
        ( add
          ( (data[n+1]-a*data[n]-b)^2,
             n=1..numelems(data)-1
          ),
          location
        );
  p1:= pointplot( data[1..-2],
                  data[2..-1],
                  symbol=solidcircle,
                  symbolsize=20,
                  color=red
                ):
  p2:= pointplot( data,
                  eval( a*data+~b,
                        fit[2,1,1]
                      ),
                  symbol=solidcircle,
                  symbolsize=20,
                  color=blue
                ):
  display( [p1, p2],
           labels=[x[n], x[n+1]],
           labelfont=[times,bold,20]
         );
 

0.2642696629e-1, {[{a = 1.035955056, b = 0.2314606742e-1}, 0.2642696629e-1]}

  

 
> 

#
# Possible values for n=30 depend
# on which of data[1]..data[13] you
# tale as the "starting" value.
#
# That's an unavoidable issue with
# using an autoregression model
#
# The following gives all possibilities
#
  f:=x->rhs(fit[2,1,1,1])*x+rhs(fit[2,1,1,2]):
  printf( "\n\n\t\tStart Index  Start Value  n=30 Value\n\n"):
  seq
  ( printf
    ( "\t\t%6d%15.3f%15.8f\n",
      j,data[j],(f@@(30-j))(data[j])
    ),
    j=1..numelems(data)
  );;



                Start Index  Start Value  n=30 Value

                     1          0.150     1.56715588
                     2          0.175     1.55763962
                     3          0.150     1.41635059
                     4          0.175     1.40748343
                     5          0.200     1.39675021
                     6          0.300     1.55937383
                     7          0.250     1.37023896
                     8          0.250     1.30033912
                     9          0.350     1.44283626
                    10          0.400     1.47175850
                    11          0.450     1.49615966
                    12          0.550     1.61074794
                    13          0.550     1.53250072

  
> 

#
# Produce a vector containing the "years". Starts
# at 1978, goes in steps of 2
#
  X:= Vector[column]([ seq
                       ( 1978+(j-1)*2,
                         j=1..numelems(data)
                       )
                     ]
                    ):
#
# Linearly fir data to year
#
  fit1:=Statistics:-LinearFit( [1, t],
                               X,
                               data,
                               t
                             );
#
# Quadratically fit data to year
#
  fit2:=Statistics:-LinearFit( [1, t, t^2],
                               X,
                               data,
                               t
                             );
#
# Plot actual data against both fits
# Bot look plausible - more or less
  p3:=pointplot( X,
                 data,
                 symbol=solidcircle,
                 symbolsize=20,
                 color=red):
  p4:=plot( [fit1, fit2],
            t=1978..1978+13*2,
            color=[blue, green]
          ):
  display([p3, p4]);
#
# Now plot the fits over the 60-year
# interval and get the end-values
#
# Which do you believe?
#
  p5:=plot( [fit1, fit2],
            t=1978..1978+30*2,
            color=[blue, green]
          ):
  display([p3, p5]);
  eval(fit1, t=1978+30*2);
  eval(fit2, t=1978+30*2);

-HFloat(34.68516483516481)+HFloat(0.017582417582417572)*t

  

HFloat(2709.8508241769696)-HFloat(2.740784215785366)*t+HFloat(6.930569430572319e-4)*t^2

  

 

 

HFloat(1.1478021978022)

  

HFloat(2.7057942057945183)

 (1)
> 

 

NULL

Download oddFit2.mw

Kaprekar??!!...

tomleslie 13876
October 26 2018
0 0

@phbmc

I'm still confused: is your objective to generate Kaprekar sequences?

I found the attached which I wrote (for fun!) a while ago. This generates Kaprekar sequences in for 4-digit numbers, in base 10, with leading zeros retained.

Rereading this, I don't think it would be too difficult to generalise it to handle different number bases, and different digit counts - but I don't propose to attempt this until I find out if that is your objective!

> 

  restart;
  Kap:= proc( x::integer)
              local toL:= proc(p)
                               local q:= p:
                               sort
                               ( [ seq
                                   ( irem(q, 10,'q'),
                                     r=1..4
                                   )
                                 ]
                               ):
                          end proc,
                    toNum:= proc(z)
                                 add( z[j]*10^(j-1),
                                      j=numelems(z)..1,-1
                                    );
                            end proc,
                    K2:= proc( x::integer, numInt:=0 )
                               local n1, n2, p:
                               n1:= toL(x);
                               n2:= ListTools:-Reverse(n1);
                               p:= toNum(n1)-toNum(n2);
                               if   p=x
                               then printf( "\n\t\tNumber of Iteration=%d\n",
                                            numInt
                                          );
                               else printf( "\t\t%4d-%4d -> %4d\n",
                                            toNum(n1), toNum(n2), p
                                          ):
                                    K2(p, numInt+1);
                               fi;
                         end proc,
              n1:= toL(x);
              if   length(x)<4
              then error "Input must have 4 digits":
              elif numelems( convert(n1, set))=1
              then error cat("Number must contain more",
                             " than 1 distinct digit"
                            ):
              else printf("\n\n\t\tIteration Sequence\n\n");
                   K2(x);
              fi;
        end proc:
  Kap(1112);
  Kap(1467);



                Iteration Sequence

                2111-1112 ->  999
                9990- 999 -> 8991
                9981-1899 -> 8082
                8820- 288 -> 8532
                8532-2358 -> 6174

                Number of Iteration=5


                Iteration Sequence

                7641-1467 -> 6174

                Number of Iteration=1

  
> 

 

Download Kaprekar.mw

Fairly simple...

tomleslie 13876
October 25 2018
0 0

but could be lengthy!

The following makes the substitution you require, and tidies up the matrix a bit.

Each entry in the matrix is a quartic in b, with coefficients which are linear combinations of the symbols x__1, x__2, y__1, and y__2. While Maple can symbolically solve such quartics, you are going to get four solutions for each matrix entry and each of these solutions is quite 'long'.

On the other hand, if x__1, x__2, y__1, and y__2 are numbers, then one could map fsolve() or roots() across the matrix, which would produce a more manageable solution.

Which do you want??

restart;
R := Vector[column](4, [a__1*b^3+a__2*b^2+(b-a__2-2)*b+b-a__1, a__1*b^3+(b-a__2-2)*b^2+a__2*b+b-a__1, (b-a__1)*b^3+a__2*b^2+(b-a__2-2)*b+a__1, (b-a__1)*b^3+(b-a__2-2)*b^2+a__2*b+a__1]):

L := Vector[column](5, [a__1*b^3+a__2*b^2+(b-a__1)*b+b-a__2-2-(b-a__2-2)*b^3-(b-a__1)*b^2-a__2*b-a__1, a__1*b^3+a__2*b^2+(b-a__2-2)*b+b-2*a__1-(b-a__1)*b^3-(b-a__2-2)*b^2-a__2*b, a__1*b^3+(b-a__1)*b^2+a__2*b+b-a__2-2-(b-a__2-2)*b^3-a__2*b^2-(b-a__1)*b-a__1, a__2*b^3+a__1*b^2+(b-a__1)*b+b-2*a__2-2-(b-a__2-2)*b^3-(b-a__1)*b^2-a__1*b, a__2*b^3+a__1*b^2+(b-a__2-2)*b+b-a__1-(b-a__1)*b^3-(b-a__2-2)*b^2-a__1*b-a__2]):

M:=Matrix(op(1,L),op(1,R), (i,j)->L[i]=R[j]):
f:=z->collect(expand(subs( [ a__1 = x__1*b + y__1,  a__2 = x__2*b + y__2], lhs(z)-rhs(z))),b)=0:
M:=f~(M);

 

Clarification needed...

tomleslie 13876
October 25 2018
0 0

I can think of ways to do this when the expression contains exactly one operand, but what 'types' would you like to be returned for expressions such as a-b+c and a*b/c?

So that would be...

tomleslie 13876
October 24 2018
0 0

@phbmc 

R := Vector[column](4, [a__1*b^3+a__2*b^2+(b-a__2-2)*b+b-a__1, a__1*b^3+(b-a__2-2)*b^2+a__2*b+b-a__1, (b-a__1)*b^3+a__2*b^2+(b-a__2-2)*b+a__1, (b-a__1)*b^3+(b-a__2-2)*b^2+a__2*b+a__1]):

L := Vector[column](5, [a__1*b^3+a__2*b^2+(b-a__1)*b+b-a__2-2-(b-a__2-2)*b^3-(b-a__1)*b^2-a__2*b-a__1, a__1*b^3+a__2*b^2+(b-a__2-2)*b+b-2*a__1-(b-a__1)*b^3-(b-a__2-2)*b^2-a__2*b, a__1*b^3+(b-a__1)*b^2+a__2*b+b-a__2-2-(b-a__2-2)*b^3-a__2*b^2-(b-a__1)*b-a__1, a__2*b^3+a__1*b^2+(b-a__1)*b+b-2*a__2-2-(b-a__2-2)*b^3-(b-a__1)*b^2-a__1*b, a__2*b^3+a__1*b^2+(b-a__2-2)*b+b-a__1-(b-a__1)*b^3-(b-a__2-2)*b^2-a__1*b-a__2]):

M:=Matrix(op(1,L),op(1,R), (i,j)->L[i]=R[j]);

 

Do loops, indexing etc...

tomleslie 13876
October 24 2018
0 0

@Ronan 

You can use a 'do' loop and create indexed table entries 'on-the-fly', but Maple provides other alternatives some of which may be 'more efficient' in certain situations and you should certainly be aware of these. Take a look at the 'toy' example in the attached - which do you prefer?

> 

#
# Toy examples
#
  a:=[1,2,3,4]:
#
# Method 1: create a table 'on-the-fly'
#
  for k from 1 to numelems(a) do
      b[k]:=a[k]^3;
  od:
  b();
#
# Method 2: Use 'seq()', producing a list
#
  b:=[seq( k^3, k in a)];
#
# Method 3: Define a function and use it
# 'elementwise'.
#
  f:= z-> z^3:
  b:=f~(a);

[1, 8, 27, 64]

  

[1, 8, 27, 64]

  

[1, 8, 27, 64]

 (1)
> 

 

Download toyIter.mw

I have modified my original response so that it no longer uses programmatically constructed names, nor even a 'do' loop - just a couple of functions which are applied 'elementwise' to lists. This a *lot* closer to the way I would *normally* code such a problem.

> 

  restart;
  with(plots):
#
# Define some functions/procs which will be used
# later
#

  Trunc := proc(F::{`+`, procedure}, odr::posint := 2, v::(list(name)) := [x, y])
                 local f;
                 description " Truncates an algebraic equation to required degree";
                 f := `if`(F::procedure, F(v[]), F);
                 if not f::`+`
                 then error "Can't truncate 1-term expression";
                 else select(q->degree(q, v) <= odr, f);
                 end if;
           end proc:
  afunc:= z->rhs~( solve
                    ( eval
                      ( [ ltg1,lprp1 ],
                        s = z
                      ),
                      [x, y]
                    )[]
                  ):
  doPlot:=z-> plot( [z[], r = -20 .. 20],
                  colour = "MediumSeaGreen"
                ):
#########################################
# Start the actual calculation
#
  C := (x^2+y^2+12*x+9)^2-4*(2*x+3)^3:
  xp1 := 1:
  yp1 := 2:
  p21 := implicitplot( C,
                       x = -3 .. 3,
                       y = -5 .. 5,
                       colour = "Salmon",
                       gridrefine = 2,
                       size = [300, 300]
                     ):

  sol1 := [solve(C, y)]:
  f1 := expand
        ( eval
          ( C,
            [x = X+r, y = Y+s]
          )
        ):
  ltg1 := expand
          ( eval
            ( Trunc(f1, 1, [Y, X]),
              [X = x-r, Y = y-s]
            )
          ):
  lprp1:= -coeff(ltg1, y)*x+coeff(ltg1, x)*y+K1:
  lprp1:= eval
          ( lprp1,
            [ K1 = solve
              ( eval
                ( lprp1,
                  [x = xp1, y = yp1]
                ),
                K1
              )
            ]
          ):

  Pedal:=afunc~(eval~(sol1, x=r)):
  display( [ p21, doPlot~(Pedal)[]]);

 
> 

  restart;
  with(plots):
#
# Define a couple of procedures
#
  Trunc:= proc(F::{`+`, procedure}, odr::posint := 2, v::(list(name)) := [x, y])
               local f;
               description " Truncates an algebraic equation to required degree";
               f := `if`(F::procedure, F(v[]), F);
               if not f::`+`
               then error "Can't truncate 1-term expression";
               else select(q->degree(q, v) <= odr, f);
               end if;
           end proc:
  Pedal:=proc( Curv, xp, yp )
             #
             # Define/initialise local variables
             # and functions
             #
               local sol1:= [ solve(Curv, y) ],
                     f1:= expand
                          ( eval
                            ( Curv,
                              [x = X+r, y = Y+s]
                            )
                          ),
                     ltg1:= expand
                            ( eval
                              ( Trunc(f1, 1, [Y, X]),
                                [X = x-r, Y = y-s]
                              )
                            ),
                     lprp1:= -coeff(ltg1, y)*x+coeff(ltg1, x)*y+K1,
                     afunc:= z->rhs~( solve
                                      ( eval
                                        ( [ ltg1,lprp1 ],
                                          s = z
                                        ),
                                        [x, y]
                                      )[]
                                    ),
                     doPlot:=z-> plot( [z[], r = -20 .. 20],
                                       colour = "MediumSeaGreen"
                                     ),
                     Ped;
             #
             # Actually do the calculation
             #
               lprp1:= eval
                       ( lprp1,
                         [ K1 = solve
                                ( eval
                                  ( lprp1,
                                    [x = xp1, y = yp1]
                                  ),
                                  K1
                                )
                         ]
                       ):
               Ped:=afunc~(eval~(sol1, x=r)):
               return doPlot~(Ped)[];
          end proc:
#######################################
# Actually perform the calculations
#
  C:= (x^2+y^2+12*x+9)^2-4*(2*x+3)^3:
  xp1:= 1:
  yp1:= 2:
  p21:= implicitplot( C,
                      x = -3 .. 3,
                      y = -5 .. 5,
                      colour = "Salmon",
                      gridrefine = 2,
                      size = [300, 300]
                    ):

  Pedal2:=Pedal(C, xp1, yp1):
  display( [ p21, Pedal2 ],
           size = [600, 600],
           scaling = constrained,
           caption = " Curve and it's Pedal curve"
         );

 

NULL

> 

 

> 

``

NULL

``

NULL


 

Download pedCur2.mw

 

An answer or two...

tomleslie 13876
October 23 2018
0 0

@Ronan 

Use of []:

Sometimes you have a list [a,b,c], but you actually want the sequence a,b,c: the construct [a,b,c][] will return the sequence form the list. This probably strikes you as a bit pointless! Sometimes, however, a Maple command will return an answer as a 'list of lists' like [[a,b],[c,d]] or even [[a,b]] when you actually want [a,b],[c,d] or [a,b]. The inclusion of [] essentiallu 'unpeels' one set of brackets from the preceding expression, so [[a,b]][] will return [a,b]. It is almost always better to avoid the extra pair of [] brackets in this first place, and this is generally possible with really careful programming - but sometimes, just tacking on the extra '[]' is a convenient, quick and dirty fix to get returned values to desired type. If anyone wanted to describe this as sloppy programming - I would not attempt to defend it!.

Programming Style

Well there are a few issues here. For my own work I always use 'worksheet mode' with '1-D input'. I just find it easier to read code that way. There is a fundamental dichotomy when writing a Maple worksheet: do you see it as 'mathematics' or 'programming'. It is actually a mixture - so neither a 'mathematics' nor a 'programming' style feels completely correct. I prefer to stick with the 'programming' style for several reasons

  1. When I started using Maple, there was no other choice: you had either 'worksheet mode' with '1-D input' or nothing. So I got used to doing it this way (and I'm probably too old to change!)
  2. Use of 2-D input adds an extra 'layer' of parsing to code: most of the time this 'parsing layer' is foolproof - but every once in a while it can cause a problem which is really hard to identify. I just prefer to avoid this possibility.
  3. Whilst I will (often) use 2-D input when responding to questions on this site, I (sometimes) find it difficult to debug stuff written this way (See items (1) and (2) above) - so my default is to convert to 1-D input, which (to me) makes things easy to read/debug

The indentation style I use is particularly 'idiosyncratic' and I would not recommend it to anyone. My only excuse is that in a previous existence, I spent years writing code in LISP (don't ask why!). As everyone knows, LISP stands for Lots of Irrelevant, Superfluous Parentheses, and the real pain in LISP was keeping track of opening/closing parentheses and making sure they 'matched'. There were various 'styles' of doing this, but aligning matching parentheses 'vertically' was one of the most basic: so rather then writing

command( arguments )

one would write

command( arguments
                )

or even

command
( arguments
)

So I'm more-or-less stuck with a LISP-y style!

For my own work, I tend to develop Maple code in a 'programmer's editor' which I have tuned to handle the indentation and parenthesis-matching automatically (more-or-less). Now, my brain is just wired to read code in this style - but I would not recommend that anyone else adopt it!

Trivial...

tomleslie 13876
October 23 2018
0 0

@Gabriel samaila 

The attached file contains a couple of (commented) additions to the previous file “R_0_fix.mw code”, Thes plot f(eta) versus eta, and f'(eta) versus eta. Note that the structure of ode means that neither of these should be affected by parameter values in the list 'L' - and they are not, whihc is good!

> 

restart
> 

with*plots; ode1 := diff(f(eta), eta, eta, eta)+(1/2)*f(eta)*(diff(f(eta), eta, eta)) = 0

diff(diff(diff(f(eta), eta), eta), eta)+(1/2)*f(eta)*(diff(diff(f(eta), eta), eta)) = 0

 (1)
> 

ode2 := (diff(theta(eta), eta, eta))*(1+(4*R*(1/3))*(theta(eta)+C[T])^3)+4*R*(theta(eta)+C[T])^2*(diff(theta(eta), eta))^2+(1/2)*pr*f(eta)*(diff(theta(eta), eta)) = 0

(diff(diff(theta(eta), eta), eta))*(1+(4/3)*R*(theta(eta)+C[T])^3)+4*R*(theta(eta)+C[T])^2*(diff(theta(eta), eta))^2+(1/2)*pr*f(eta)*(diff(theta(eta), eta)) = 0

 (2)
> 

bcs1 := f(0) = 0, (D(f))(0) = 0, (D(f))(10) = 1

f(0) = 0, (D(f))(0) = 0, (D(f))(10) = 1

 (3)
> 

fixedparameter := [pr = .72, a = .6, C[T] = .2]

[pr = .72, a = .6, C[T] = .2]

 (4)
> 

ode3 := eval(ode2, fixedparameter)

(diff(diff(theta(eta), eta), eta))*(1+(4/3)*R*(theta(eta)+.2)^3)+4*R*(theta(eta)+.2)^2*(diff(theta(eta), eta))^2+.3600000000*f(eta)*(diff(theta(eta), eta)) = 0

 (5)
> 

bcs2 := theta(10) = 0, (D(theta))(0) = -a*(1-theta(0))

theta(10) = 0, (D(theta))(0) = -a*(1-theta(0))

 (6)

####### now i want to vary 'R' instead of   'a' that is kepping 'a' and varying 'R'

> 

L := [0., 0.5e-1, .1, .2, .4, .6, .8, 1, 5, 10, 20]

[0., 0.5e-1, .1, .2, .4, .6, .8, 1, 5, 10, 20]

 (7)
> 

#
# Check the indeterminates in the system. Note that the
# parameter 'a' is undefined. Since this parameter is ued
# in the definition of D(theta)(0), this boundary condition
# is also undefined, and hence the system cannot be solved
#
  indets(eval({bcs1, bcs2, ode1, ode3}, R = L[1]), 'name');

{a, eta}

 (8)
> 

for k to 10 do sol_All := dsolve(eval({bcs1, bcs2, ode1, ode3}, R = L[k]), [f(eta), theta(eta)], numeric, output = listprocedure); Y_sol || k := rhs(sol_All[5]); YP_sol || k := -rhs(sol_All[6]) end do

Error, (in dsolve/numeric/bvp/convertsys) too few boundary conditions: expected 6, got 5

  
> 

  #### do you dow why this errow massage?
> 

#
# Reason for error message is explained above.
#
# Rerun the above, but this time set the value of a=1, just
# as a check. I have no idea what the value you want!!
#
  for k to 10 do
      sol_All := dsolve
                 ( eval
                   ( {bcs1, bcs2, ode1, ode3},
                     [R = L[k], a=1]
                   ),
                   [f(eta), theta(eta)],
                   numeric,
                   output = listprocedure
                 );
      Y_sol || k := rhs(sol_All[5]);
      YP_sol || k := -rhs(sol_All[6]);
#
# Addition:- Extract f(eta) and f'(eta) from
# the returned results
#
      feta || k := rhs(sol_All[2]);
      fpeta || k := rhs(sol_All[3])
  end do:
> 

for k to 10 do L[k], [(Y_sol || k)(0), (YP_sol || k)(0)] end do

0., [HFloat(0.7718221287933378), HFloat(0.22817787120666155)]

  

0.5e-1, [HFloat(0.7794871400247606), HFloat(0.220512859975239)]

  

.1, [HFloat(0.7867583087402996), HFloat(0.21324169125970144)]

  

.2, [HFloat(0.8001174320465659), HFloat(0.19988256795343393)]

  

.4, [HFloat(0.8223999433940956), HFloat(0.17760005660590436)]

  

.6, [HFloat(0.8397214606309522), HFloat(0.16027853936904765)]

  

.8, [HFloat(0.8532960930254505), HFloat(0.14670390697454935)]

  

1, [HFloat(0.8641178449875935), HFloat(0.13588215501240694)]

  

5, [HFloat(0.9291481289213374), HFloat(0.0708518710786633)]

  

10, [HFloat(0.947214313258223), HFloat(0.052785686741776396)]

 (9)
> 

NULL
> 

#
# Plot the output obtained when the value of 'a'
# is set (arbitrarily) to 1
#
  plot( [ seq((Y_sol||j)(eta), j = 1..10)],
         eta = 0 .. 6,
         labels = [eta, theta(eta)],
         axes = boxed
      );
#
# Addition:- plot f(eta) vs eta and f'(eta) versus
# eta. Note that neither of these depend on the values
# the list L. This is to be expected because the ODE
# determining them (ode1 above) is independent of
# L-values
#
  plot( [ seq((feta||j)(eta), j = 1..10)],
         eta = 0 .. 6,
         labels = [eta, f(eta)],
         axes = boxed
      );
  plot( [ seq((fpeta||j)(eta), j = 1..10)],
         eta = 0 .. 6,
         labels = [eta, fprime(eta)],
         axes = boxed
      );

 

 

 
> 

 

> 

 

Download moreplots.mw

Observation...

tomleslie 13876
October 18 2018
0 0

@Carl Love 

Seriously doubtful population equations - I accept that methematically one can set these to absolutely anything, but (just for fun) let's assume that these equations bear some relation to reality.

The fox population is given by

k := k0*(-b*r0+a+1);

so if there are no rabbits (ie r0=0), and with a=0.15, then depending on how one handles hte "fractional fox" issue,either the fox population stays the same, or it grows slowly. Maybe they've become vegetarian!?

In a similar vein, the rabbit population is given by

r := r0*(b*e*k0-c+1);

so if there are no foxes (k0=0), and with c=0.02, then again depending on how one handles the "fractional rabbit" issue, either the rabbit population stays the same, or it grows very slowly. Trust me, that isn't how rabbits breed!!!

 

 

Well...

tomleslie 13876
October 18 2018
0 0

@brian bovril 

You can get a "general expression", but it involves summing over the roots of a cubic, and two of these roots are complex. Seems (to me) amazing that such a calculation can result in non-complex integers - but then (sometimes) Maple is magic!

See the attached

> 

restart;
a:=unapply(rsolve({{U(n)=U(n-1)+3*U(n-3), U(0) = 6, U(1) = 6, U(2) = 6, U(3)=24}}, U(n)), n);  # General formula
seq(expand(a(n)), n=0..9);  # The first 10 terms

Error, (in rsolve) First argument must be an equation or set of equations

  

a(0), a(1), a(2), a(3), a(4), a(5), a(6), a(7), a(8), a(9)

 (1)
> 

#
# Get the general formula
#
  U:= unapply
      ( rsolve
        ( { U(n)=U(n-1)+3*U(n-3), U(0) = 6, U(1) = 6, U(2) = 6, U(3)=24},
          U
        ),
        n
      );

proc (n) options operator, arrow; Sum(6*(1/_R)^n/((9*_R^2+1)*_R), _R = RootOf(3*_Z^3+_Z-1)) end proc

 (2)
> 

#
# The above "solution" is a bit WTF?! So just check
# what values it generates
#
  seq(evalf(U(j)), j=0..10);
#
# If one eliminates the imaginary parts - they are
# all zero, and then round() the real parts, then
# hey presto
#
   seq(round(Re(evalf(U(j)))), j=0..10);

6.000000000+0.*I, 6.000000000+0.*I, 6.000000000+0.*I, 24.00000000+0.*I, 42.00000000+0.*I, 60.00000000+0.*I, 132.0000000+0.*I, 258.0000000+0.*I, 438.0000000+0.*I, 834.0000000+0.*I, 1608.000000+0.*I

  

6, 6, 6, 24, 42, 60, 132, 258, 438, 834, 1608

 (3)
> 

#
# But what exactly does that sum() in the original
# solution actually mean?
#
# Well it is a sum over the roots, What are these
# roots??
#
  rootVals:=[allvalues(op([2,2],U(n)) )];
#
# Regenerate the series by explicitly summing over
# the rootVals. Note that no floating point arithmetic
# is used (or necessary)
#
  seq( simplify(evalc(add( op(1,U(n)), _R in rootVals))), n=0..10);

[(1/6)*(36+4*85^(1/2))^(1/3)-(2/3)/(36+4*85^(1/2))^(1/3), -(1/12)*(36+4*85^(1/2))^(1/3)+(1/3)/(36+4*85^(1/2))^(1/3)+((1/2)*I)*3^(1/2)*((1/6)*(36+4*85^(1/2))^(1/3)+(2/3)/(36+4*85^(1/2))^(1/3)), -(1/12)*(36+4*85^(1/2))^(1/3)+(1/3)/(36+4*85^(1/2))^(1/3)-((1/2)*I)*3^(1/2)*((1/6)*(36+4*85^(1/2))^(1/3)+(2/3)/(36+4*85^(1/2))^(1/3))]

  

6, 6, 6, 24, 42, 60, 132, 258, 438, 834, 1608

 (4)
> 

 

> 

 

Download recur.mw

 

Some answers...

tomleslie 13876
October 18 2018
0 0

@David1 

In all that follows remember that I am far from being an expert on optimization

  1. Nothing *wrong* with using the "solutions=1" option. AFAIK, this only changes the output display. So won't save any time. One argument for not using this option is that one can compare the various optima obtained and determine whether the "global" one is "better" than the next few on the output list. Think of a sine wave: all the local optima are equally good: which one should be reported as the global optima (and would you care?). Now in some multidimensional problem, would you want to know it multiple identical(-ish) are obtained? Since it costs nothing, I tend not to restrict the output to reporting one solution
  2. The final figure in the DirectSearch() output is the number of "major iterations" taken to produce this partcular line in the output. Two factors (at least) contribute to this
    1. Depending on the "method" being used there will be a different number of function evaluations on "test points" adjacent to the "current" point, in order to figure out which direction to move: possibly also by how far to move in each variable, since the objective function may be "slowly-varying" in one variable and "rapidly-varying" in other
    2. The dimensionality of the problem: the more variables there are, the more possibilities there are to check in the above
  3. Notice that each line (ie "solution") in the output produces a different answer for this final entry. My interpretation of this is as follows. In order perform a "global" search, what most(?) optimizers actually do, is to use many different "starting" points, which may produce many different "finishing" points (think - local minima). The "global" optimum is then selected as the best of the "local" optima. AFAIK, the final number in each line of the output is the number of major iterations (see above) from each of the associated "starting" points
  4. Time taken: this has a number of contributing factors, and it is generally "risky" to compare optimizers on this criterion alone (unless one is doing the comparison across a large number of problems. There are two major contributing factors to the total time
    1. How many different starting points are used
    2. How many iterations are attempted from each starting point before deciding to "give up" and try a new starting point

One can usually "tune" both of the above - but you have to be aware of what the optimizer does "by default"? If two optimizers have produced the same answer, but one has tested (say) 100 starting points, and the other 200 (taking twice as long), then the former is definitely "quicker" but the latter has done a better check that the obtained result is truly "global". Which is "better"?

If you think of "starting" points as a grid, then their total (default) number should probably(?) be linearly(?) dependent on the number of variables in the problem - but some optimizers (probably?) don't do this!

If the user doesn't define ranges for each of the variables, then what ranges of starting points should be covered for each variable - +/-10, +/-100, or +/-10000000000. In an "ideal" world, one would like to check the finest possible grid of starting points over the widest possible ranges - but that would get seriously time-consuming

Tuning the number of iterations from each starting point is similarly "subjective". Most optimizers work on the basis that if the difference between the "current" point and any feasible "next" point is less than epsilon then stop. A moment's thought should convince you that the value of epsilon should be "relative" to the magnitude of the objective value: but should it be 1%, 0.1%, 0.01%.......The larger you make this number then the fewer iterations you will perform (so faster) but going smaller (so slower) *may* produce a "better" solution. Again - which is better?

Other factors which can significantly contribute to time taken include

how hard is it to evaluate the objective function?

how hard is it to evaluate constraints - because at each function evaluation one has to check whether or not there is a constraint violation: if there is, then what should be done?

hardware/software floats. Pretty much all optimizers use hardware floats. So will DirectSearch(), unless you have set the Digits() variable to a quantity sufficiently high, in which case it will use software floats - this will slow things down (a lot!)

I have just done a bit of reading about yalmip. Interestingly this is not technically an optimiser. It is a consistent (Matlab) interface which can call one of a number of different, separately-installed optimizers: yalmip's authors refer to these as low-level solvers. Whether the authors of these optimizers would agree with this definition, is a different matter!

 

Syntax/logical issue...

tomleslie 13876
October 17 2018
0 0

related to the definition of variables in the function 'f()'. Do this correctly and you will get answer - see the attached.

Bearing in mind all the comments I (and others) have made on this thread regarding how much one can *trust* optimizers for globally solving non-convex problems, I wouldn't want to guarantee that this is a *the* global optimum, but it is probably a pretty good guess. It is also unfair to blame DirectSearch() when the user has incorrectly specified the problem!

restart

with(DirectSearch)

with(LinearAlgebra)

A := RandomMatrix(5, 3)

b := RandomVector(5)

``f := proc (x1, x2, x3) options operator, arrow; Norm(A.`<,>`(x1, x2, x3)-b, 1) end proc

proc (x1, x2, x3) options operator, arrow; Norm(Typesetting:-delayDotProduct(A, `<,>`(x1, x2, x3))-b, 1) end proc

 (1)

constr := x1^2+x2^2+x3^2 = 1

x1^2+x2^2+x3^2 = 1

 (2)

GlobalSearch(f, {constr})

Array(%id = 18446744074601136486)

 (3)

``

Download how_to_use_directSearch.mw

Before updating...

tomleslie 13876
October 17 2018
0 0

@Thomas Richard 

Output of Physics:-Version() was

"C:\Users\TomLeslie\maple\toolbox\2018\Physics Updates\lib\Physics Updates.maple", `2018, July 24, 12:24 hours, version in the MapleCloud: unable to determine, version installed in this computer: 78`

After updating, output of Physics:-Version() is

"C:\Users\TomLeslie\maple\toolbox\2018\Physics Updates\lib\Physics Updates.maple", `2018, October 17, 14:31 hours, version in the MapleCloud: unable to determine, version installed in this computer: 145`

And the matrix multiplication problem has "gone away"

  • How?
  • Why?
  • The Physics package can "break" 2D-input matrix multiplication !!!!

Can you...

tomleslie 13876
October 17 2018
0 0

specify OS/build version?

I think I found your issue...

tomleslie 13876
October 17 2018
0 0

@ebrahimina 

It may such a serious Maple bug that I wanted others to examine/verify, so I created a 'trivial' illustration and posted as a new question as below

https://www.mapleprimes.com/questions/225701-Matrix-Multiplication-With-2D-Input

First 82 83 84 85 86 87 88 Last Page 84 of 207