tomleslie

13876 Reputation

20 Badges

15 years, 175 days

MaplePrimes Activity


These are answers submitted by tomleslie

Making a wild guess. When Maple says

"Error, (in Bits:-GetBits) argument 1, the number, must be a nonnegative integer"

maybe it means that the first argument should be a non-negative integer, Now I have actually tried this with positive and negative integer arguments; GetBits() handles either. However, if the first argument is a float, then I get the above Error Message. Since you can't be bothered to post the code which leads to the above error message, I cannot check further - but I'm betting that you are prociding a float as first argument.

If you are not, then post the code which leads to the above error so I can debug properly

is probably to use the dismantle() command, as in

f:=x^2*exp(3/2)*sin(Pi/3-1/x);
dismantle(f);

although you do have to read the output of the dsimantle() command quite carefully to appreciate the details

 

You have been unfortunate enough to generate somethiing called an ill-conditioned polynomial. The final stage of your calculation, the computation of the determinant, produces a polynomial in P of degree 18, for which you require the roots. Unfortunately, in an ill-conditioned polynomial the roots vary a lot if any of the coefficients vary by even a tiny amount. This (rather counter-intuitive) phenomenon was first(?) reported by a guy called Wilkinson, who illustrated it with what (at first sight) seems to be a really "innocuous" polynomial. See the article at

https://en.wikipedia.org/wiki/Wilkinson%27s_polynomial

This will probably put you off believing the results of any polynomial roots calculation forever! However all is not lost,yet. One way to get around this problem is to avoid floating point numbers at all costs. Just don't use floats! One can achieve this by converting all of your parameter values to exact rationals. Consider the simple case of 1/3: it is an exact rational, but this is not the same as 0.33333333 or 0.333333333333333333 or 0.333333333333333333333333333333. You may claim that you are not using anything parameters with recurring decimals, but when you write Ec := 0.380e12, whihc will be stored as a floating point number, are you expecting it to be stroed (and used) as

0.380000000000000000000000000000000e-12 or maybe
0.3800000000000000000000000000000000000000000000000000000000000000000000000000000e-12, or maybe even
0.3800000000000000000000000000000000000010000000000000000000000000000000000000000e-12

because when such a number is used in the computation of coefficients of an ill-conditioned polynomial this is (probably) going to make a difference

By coverting all parameters to exact rational fractions in Maple, all subsequent calculations will then be carried out using a combination of integers and exact rationals, completely avoiding the issue of their representation as floating point numbers. The Maple woiksheet for doing this is attached later: unsurprisingly, this gives consistent results, irrespective of the number of digits requested for the final calculation step.

Comparing this calculation with that returned by Matlab is difficult  (for me!). There are (at least) two symbolic calculators which can be used within Matlab. Once upon a time, Matlab used a "stripped-down" version of Maple as its symbolic toolbox: but Maple and Matlab had an argument, and Matlab swiched to using a different engine called muPad, I believe) for symbolic calculations. Your "Matlab" code only functions with the muPad engine. My Matlab installation is set up to use Maple as the symbolic calculation engine for the simple reason that having paid for Maple and Matlab, I could find no reason to pay even more for the "Matlab symbolic toolbox" (aka muPad). Since I can use Matlab with Maple as the symbolic engine - why bother.

For your specific Matlab code there are only a few 'symbolic' commands - 'solve', 'diff'', and 'subs' being the obvious ones. The functionality of these commands in MAple and MuPad is similar but their syntax is (slightly) different. Hence if I try to run your "Matlab" code, I just get numerous syntax errors. I could probably rewrite your "Matlab" code to do all of the symbolic calculations in Maple, but since this is the same as using Maple directly - what would be the point? Hence I cannot really comment on the best way to use muPad to deal with the roots of ill-conditioned polynomials. That is probably a question for a Matlba forum - not this one

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

Code for Maple calculation

  restart:
  with(LinearAlgebra):
  beta:=convert(1, rational, exact):          nu:=convert(0.3, rational, exact):
  lambda:=convert(2, rational, exact):        G:=convert(5, rational, exact):  
  ko:=convert(.5*0.1e7, rational, exact):     Ec:=convert(0.380e12, rational, exact):
  Em:=convert(0.70e11, rational, exact):      ri:=convert(0, rational, exact):
  ro:=convert(0.5, rational, exact):          ti:=convert(0.1e-1, rational, exact):
  K:=convert(0.4e-1*0.1e10, rational, exact): n:=convert(1, rational, exact):
  m:=convert(20, rational, exact):            alpha:=convert(0.1, rational, exact):
  t:= ti*(1+alpha*r/ro)^beta:
  Er:= (Ec-Em)*(r/ro)^n+Em:
  W:= simplify
      ( add
        ( a[n]*ChebyshevT(n, r),
          n = 0 .. m
        )
      ):
  sys:= { eval
          ( W,
            r = ro
          ),
          eval
          ( diff(W, r),
            r = ri
          )
        }:
  W:= subs
      ( solve
        ( sys,
          {a[0], a[1]}
        ),
        W
      ):
  d1:= diff(W, r):
  d2:= diff(d1, r):
  W:= subs
      ( solve
        ( { eval
            ( ko*d1-Ec*ti^3*(1+alpha)^(3*beta)*(d1*nu+d2*r)/(12*r*(-nu^2+1)),
              r = ro
            )
            = 0
          },
          a[2]
        ),
        W
      ):
  d1:= diff(W, r):
  d2:= diff(d1, r):
  Uf:= int
       ( K*(1+lambda*r/ro+G*(r/ro)^2)*W^2*r,
         r = ri .. ro
       ):
  NUM:= int
        ( ( ( d2+d1/r)^2-(2*(1-nu))*d2*d1/r)*Er*t^3*r/(12*(-nu^2+1)),
          r = ri .. ro
        )
        +
        ko*ro*eval
              ( d1,
                r = ro
              )^2
        +
        Uf:
  DEN:= int
        ( d1^2*r,
          r = ri .. ro
        ):
  PT:= NUM-Em*ti^3*P*DEN/ro^2:
  for c from 3 to m do
      eq[c]:= diff(PT, a[c])
  end do:
  MT:= GenerateMatrix
       ( [ seq
           ( eq[j],
             j = 3 .. m
           )
         ],
         [ a[j]$j = 3 .. m ]
       )[1]:
  ans:= solve
        ( Determinant(MT)
        ):
  for k from 10 by 10 to 100 do
      evalf[k](ans[1]):
  od;

10.80811852

 

10.808118523654378125

 

10.8081185236543781253355008657

 

10.80811852365437812533550086568014055867

 

10.808118523654378125335500865680140558670033046240

 

10.8081185236543781253355008656801405586700330462399584643526

 

10.80811852365437812533550086568014055867003304623995846435264865883824

 

10.808118523654378125335500865680140558670033046239958464352648658838238354076898

 

10.8081185236543781253355008656801405586700330462399584643526486588382383540768981114494012

 

10.80811852365437812533550086568014055867003304623995846435264865883823835407689811144940120067749174

(1)

 

Download illCon.mw

 

Assumes that "sensitivity" wrt each parameter is as defined in Rouben Rostamian's response to the OP' previous post

restart;

#
# Set up numerical values for all problem parameters
#
  params:=[       k=.5,         tau=.95,      mu=0.1e-1,
                 pi=116.1, vartheta=0.8e-2,  phi=0.25e-2,
            epsilon=0.2e-2,     rho=0.5e-1, beta=0.115e-1,
                chi=0.598e-2,     q=.5,      eta=.2,
              delta=.1,       alpha=0.57e-1,   p=.2,
            Upsilon=1.2
          ]:

#
# Define main function
#
  R:= k*tau*(rho*(Upsilon*(mu+alpha+eta)+chi)+(1-rho)*(Upsilon*(1-q)*eta+mu+beta+chi))
      *
      (pi*(-mu*p+mu+phi)/(mu*(vartheta+mu+phi))+epsilon*pi*(mu*p+vartheta)
      /
      (mu*(vartheta+mu+phi)))/((mu+beta+chi)*(mu+alpha+eta)-chi*eta*(1-q));

k*tau*(rho*(Upsilon*(mu+alpha+eta)+chi)+(1-rho)*(Upsilon*(1-q)*eta+mu+beta+chi))*(pi*(-mu*p+mu+phi)/(mu*(vartheta+mu+phi))+epsilon*pi*(mu*p+vartheta)/(mu*(vartheta+mu+phi)))/((mu+beta+chi)*(mu+alpha+eta)-chi*eta*(1-q))

(1)

#
# Compute "all" derivatives and evaluate numerically.
#
# For the purposes of this calculation "all"
# derivatives, means the derivatives with respect to
# every variable returned by indets(R, name)
#
# Output a list of two element lists where each of
# the latter is
#
# [ varName,
#   eval( diff( R, varName), params )
# ]
#
 [ seq( [j, eval( diff( R, j), params )],j in indets(R, name))];

[[Upsilon, 45499.94339], [alpha, -242658.4807], [beta, -2203579.010], [chi, -1207856.747], [epsilon, 62441.36156], [eta, 25848.6895], [k, 131376.6247], [mu, -7221570.068], [p, -62316.47885], [phi, 3039828.235], [pi, 565.7908043], [q, -107402.8680], [rho, 75126.34862], [tau, 69145.59197], [vartheta, -3191819.648]]

(2)

#
# Compute all "sensitivities" (where the sensitivity
# is as defined in Rouben Rostamian response to the
# OP's earlier post) and evaluate numerically.
#
# For the purposes of this calculation "all" sensitivities
# means the sensitivity with respect to every variable
# returned by indets(R, name)
#
# Output a list of two element lists where each of
# the latter is
#
# [ varName,
#   eval( varName*diff( R, varName)/R, params )
# ]
#
  seq( [j, eval( j*diff( R, j)/R, params )],j in indets(R, name));

[Upsilon, .8311970593], [alpha, -.2105630804], [beta, -.3857788046], [chi, -.1099584247], [epsilon, 0.1901140685e-2], [eta, 0.7870103093e-1], [k, 1], [mu, -1.099369098], [p, -.1897338403], [phi, .1156913660], [pi, 1.000000000], [q, -.8175188560], [rho, 0.5718395397e-1], [tau, 1], [vartheta, -.3887229899]

(3)

``

Download sens.mw

It isn't clear what "data" you want to write. You have two main options

  1. Generate a plot of the solution surface (over ranges you specify) and then use Maple's exportplot() command
  2. Generate a Matrix of the solution at specified x,y values, and output this using Maple's ExportMatrix() command

It seems as if you may want the latter. See the attached

NULL

  restart;
  eq := diff(Tf(x, y), x) = diff(Tf(x, y), y, y):
  Cl := {1 = Tf(0, y), (D[2](Tf))(x, 0) = 1, (D[2](Tf))(x, 1) = 1}:
  pds := pdsolve(eq, Cl, numeric):
#
# generate a plot of the solution over the specified
# xRange yRange, Can then use Maple's exportplot()
# command to export in a variety of formats. See
# help page at ?exportplot
#
  p1:= pds:-plot3d( x = 0 .. 1,
                    y = 0 .. 1
                  );
#
# Generate a matrix of solution values at specified [x,y]
# pairs. Can then use Maple's ExportMatrix() command to
# export in a variety of formats. See help page at
# ?ExportMatrix
#
  M:= Matrix
      ( [ seq
          ( [ seq
              ( rhs(pds:-value(x = k)(j)[3]),
                j=0..1,0.2
              )
            ],
            k=0..1,0.2
          )
        ]
      );

 

_rtable[18446744074332911126]

(1)

 

 

``

Download pdExport.mw

Can confirm the problem, but I'm not sure why it happens. Obvious suggestion is that it is something to do with the "CONCATENATE" formula. Maple usually doesn't mind cell values defined by formulae in Excel - at least when they produce numeric values.

Workarounds

1. Import only the first couple of columns of data and do the concatenation in Maple, as in

restart;
L:= ExcelTools:-Import("C:/Users/TomLeslie/Desktop/Book1.xls",1, "A3:B8"):
M:= Matrix( op([1,1],L),
                    3,
                    (i,j)->`if`( j<3,
                                   L[i,j],
                                   cat(L[i,1], ",", L[i,2])
                                 )
                 );
V:=Vector[column]( op([1,1],L), i->cat(L[i,1], ",", L[i,2]) );

depending on whether you want an n*3 matrix, or just a vector of the concatenated data.

2. Export the data from the Excel file as "plain vanilla csv", then use Maple's standard "import" function, with target = matrix, then select the relevant data, as in

L:=Import("C:/Users/TomLeslie/Desktop/Book1.csv", output=Matrix);
M:=Matrix( L[2..6,1...3]);
M:=Vector[column]( L[2..6,3]);

again depending on whether you want an n*3 matrix, or just a vector of the concatenated data.

Read comments in the file.

You will have to change the 'destination' string to something appropriate for your machine

#
# Define the destination directory for the
# Excel files. OP will have to change this
# to something appropriate for his/her
# machine
#
  destination:="C:/Users/TomLeslie/Desktop/";
#
# Generate a random 10*10 matrix for test
# purposes
#
  myMat:=LinearAlgebra:-RandomMatrix(10,10);
#
# Define a coupel of headers to be included
# in the output Excel file
#
  headers:=[ "data1", "data2"]:
#
# Output the data, two columns at a time to
# each of five different files
#
  for j from 1 to 9 by 2 do
    #
    # Generate the filename. Will result in
    # names
    #
    #  myXL1.xlsx,
    #  myXL2.xlsx,
    #      .
    #      .
    #  myXL5.xlsx,
    #
      fname:= cat( destination,
                   "myXL",
                   convert((j+1)/2, string),
                   ".xlsx"
                 ):
    #
    # Export the headers to the location
    #
    #  worksheet=sheet1,
    #  starting cell = B2
    #
    # in the specified file
    #
      ExcelTools:-Export( headers,
                          fname,
                          "sheet1",
                          "B2"
                        );
    #
    # Export the data to the location
    #
    #  worksheet=sheet1,
    #  starting cell = B4
    #
    # in the specified file
    #
      ExcelTools:-Export( myMat[.., j..j+1],
                          fname,
                          "sheet1",
                          "B4"
                        ):
  od:
 

 


Download exportXL.mw

it depends on the nature of the calculations being performed.

So far as I can tell the intel i7-6500U has two cores, with two threads per core, so if you execute the Maple command

kernelopts(numcpus);

the answer ought to be 4 (=numberOfCores*numberOfThreads), For comparison, I'm running an intel i7-3770K which has 4 cores and two threads/core, so in my case

kernelopts(numcpus);

returns 8.

As a general rule, when executng a Maple worksheet, only one thread will be used - so (in your case) you will only be using a maximum of 25% of your processor's capabilities (or for my machine 12.5%).

So can more of of the processor capabilities be utilised? Yes! But it depends on whether or not the calculation you are performing can be parallelized. Some calculations can be, and some can't. Maple will not decide this for you - you just have to know and program appropriately!

If calculations can be parallelized then your best option is probably the Threads:-Task() sub-package. The Parallel Programming chapter in Maple's ProgrammingGuide.pdf is a useful introduction. If you don't have the pdf then ?ProgrammingGuide,Chapter15 at the Maple prompt will give more or less the same info. Careful reading is recommended

The following works (ie no errors) in Maple 2017

restart;
int(int(sin(x^2*y), x = 1 .. Pi), y = 1 .. Pi);
int(int(sin(x^2*y), x = 1 .. evalf(Pi)), y = 1 .. evalf(Pi));

However the second command above fails in Maple 2015 with the error message

Error, (in type/laurent) too many levels of recursion

A simple change to

restart;
int(int(sin(x^2*y), x = 1 .. Pi), y = 1 .. Pi);
evalf(%);

and it will work in Maple 2015 also

 

 

that what you want is to obtain an expression of the form A+B*t, where A and B are not dependent on the variable't'. In order to achieve this, the function psi1(t) has to be defined. As a 'toy' example. if I make psi1() a simple quadratic in 't', and assume that you want the Taylor series expansion about t=0, then the following code will achieve what you want

restart:
psi1:= z->3*z^2+2*z+1:
fun:= sin(phi2)*sin(psi1(t))*lk*m1*diff(psi1(t), t, t)+sin(phi2)*diff(psi1(t), t)^2*cos(psi1(t))*lk*m1;
taylor( fun, t=0, 2);
convert(%, polynom);

You can redefine the function psi1() to be more or less anything well-behaved, and expand the Taylor series about values other than t=0, should you so desire with trivial amendments to the above

You *ought* to be able to get this plot using the syntax in your original post, but as you have discovered -you can't, and I don't understand why!

I suppose (to be fair), I should point that the help page 'Plot with Units' only refers to plotting with functions or procedures - not points.

I tried several variations, none of which "worked", so I eventually decided that the simplest (?!) way is to extract values and units from your data, then plot the values, and add the units to the axes. This is painful, but for what it is worth, the following ought to work

restart;
my_scaled_points := [ [  1.0*Unit('kg'), 2.0*Unit('hour') ],
                                       [ 21.0*Unit('kg'), 7.0*Unit('hour') ]
                                    ];
getV:= z -> [ seq( map( j -> op(1,j), k ), k in z ) ]:
getU:= z -> map( j -> op(2,j), z[1] ):
plot( getV( my_scaled_points ), # get values
         style = point,
         symbol = solidcircle,
         symbolsize = 24,
         useunits = getU( my_scaled_points ) # get units
     );

There really ought to be a better way, and I hope someone can provide it.

 

You may want to do z-transforms from the basic definition for teacjing/learning purposes, which is fine. On the other hand if you just wnat to know the answers, then the optimum strategey is probably to use Maple's in-built ztrans() function. For y0ur two examples the code would be

   restart;
   ztrans( 1, n, z );
   ztrans( exp(a*k*Ta), k, z );
#
# Get the above in "standard" form
#
   expr:= simplify(normal(%));

 

becuase with Maples built-in routines

Optimization:-Maximize(x*y+y*z, [ x*y=1, y^2+z^2=1]);

returns

[1.49999999999997780, [x = 1.41421356249703, y = .707106781124551, z = .707106781248575]]

and

Optimization:-Minimize(x*y+y*z, [ x*y=1, y^2+z^2=1]);

returns

[.499999999989945598, [x = 1.41421339391716, y = .707106865409386, z = -.707106696967656]]

 

Below has fixes for your original code

# Letter L rotates in 2D only
restart:
with(plots):
with(plottools):
# 2D rotation only
thet:=Pi/3:
# lh := (x, y) -> [x*cos(thet),y]:  ???
#lh(1,2);
x0:=0:y0:=0:w:=2:h:=12:base:=7:ang:=1:
L1r:=[[x0*ang,y0],[x0*ang,y0+h],[(x0+w)*ang,y0+h],[(x0+w)*ang,y0]]:
L2r:=[[x0*ang,y0],[x0*ang,y0+w],[(x0+base)*ang,y0+w],[(x0+base)*ang,y0]]:
#NB The above rectangles overlap
q[0] := polygon(L1r):
r[0] := polygon(L2r):
print(`Single letter L`);
plots[display](q[0],r[0], style=patchnogrid, scaling=constrained,color=green);  

ang := 0: c1 := 1:
while evalf(ang-2*Pi) < 0 do

Lu:=[[x0*cos(ang),y0],[x0*cos(ang),y0+h],[(x0+w)*cos(ang),y0+h],[(x0+w)*cos(ang),y0]]:
q[c1]:=polygon(Lu, color=grey): #changed index to c1
Lb:=[[x0*cos(ang),y0],[x0*cos(ang),y0+w],[(x0+base)*cos(ang),y0+w],[(x0+base)*cos(ang),y0]]:
r[c1]:=polygon(Lb, color=green):#changed index to c1
#rotate
  ang := ang+Pi/4:
#
# swap order of the following two statements
#
pl[c1]:=plots[display](q[c1],r[c1], style=patchnogrid, scaling=constrained,color=green):
c1 := c1+1:
end do:
print(`Animated seq of 2D L`);
plots[display](seq(  pl[i], i=1..c1-1), insequence=true,style=patchnogrid, scaling=constrained);



 

 


Download 2DLfixed.mws

Still think that the following (slightly revised version of the one I posted for your earlier question) is better

  restart;
  with(plots):
  with(plottools):
  x0:=0: y0:=0: z0:=0: w:=2: h:=12:base:=7:
#
# 'Draw' the letter L
#
  L:=polygon( [ [    x0,    y0,  z0 ], [ x0+base, y0,  z0  ],
                  [x0+base, y0, z0+w], [   x0+w,  y0, z0+w ],
                  [ x0+w,   y0, z0+h], [   x0,    y0, z0+h ]
                ]
             ):
#
# 'Draw' the letter H in the x-z plane
#
  H:=polygon( [ [    x0,     y0,     z0     ], [   x0+w,    y0,     z0     ],
                [   x0+w,    y0, z0+h/2-w/2 ], [ x0+base-w, y0, z0+h/2-w/2 ],
                [ x0+base-w, y0,    z0      ], [  x0+base,  y0,     z0     ],
                [  x0+base,  y0,   z0+h     ], [ x0+base-w, y0,    z0+h    ],
                [ x0+base-w, y0, z0+h/2+w/2 ], [   x0+w,    y0, z0+h/2+w/2 ],
                [   x0+w,    y0,   z0+h     ], [    x0,     y0,    z0+h    ]
                ]
             ):
#
# Define rotation function, accepts letter polygon,
# (or a combination of them), rotation axis and
# number of frames
#
  doRot:= (let, ax, nf) -> display
                           ( [ seq
                               ( rotate
                                 ( let,
                                   ax[1]*j*2*Pi/nf,
                                   ax[2]*j*2*Pi/nf,
                                   ax[3]*j*2*Pi/nf
                                 ),
                                 j=0..nf
                               )
                             ],
                             insequence=true,
                             scaling=constrained
                           ):
###################################################
#             Examples
#
# Rotate 'L' about x-axis with 24 frames
#
  doRot(L, [1,0,0], 24);
#
# Rotate 'H' about y-axis with 24 frames
#
  doRot(H, [0,1,0], 24);
#
# Combine 'L' and 'H' in a single plot
# Rotate the combination about each of
# three axes with 24 frames
#
  comb:= display( [L, translate(H, base+w, 0, 0)],
                  scaling=constrained,
                  view=[ -2*base..2*base,
                         -2*base..2*base,
                         -h..h
                       ]
                ):
  doRot( comb, [1,0,0], 24);
  doRot( comb, [0,1,0], 24);
  doRot( comb, [0,0,1], 24);

 


Download letrot2.mw

The attached shows how I *might* do this if I was dealing only with individual (and combinations of) non-curved letters.

Can't (at the moment) think of a 'simple' way to extend this approach to include 'curvy' letters such as 'D', 'O', etc:-(

This was created (and works) in Maple2017, but I can't think of any reason why it wouldn't work in Maple7; it only uses pretty basic commands from the plottools() and plots() package

letrot.mw

First 142 143 144 145 146 147 148 Last Page 144 of 207