acer

32358 Reputation

29 Badges

19 years, 331 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

restart;

eq16 := r(t) = d[vol]*V/(KUS*V^2+L*tau);

r(t) = d[vol]*V/(KUS*V^2+L*tau)

F := (u,v) -> sort(expand(u/v),order=tdeg(v),ascending):

lhs(eq16) = F(numer(rhs(eq16)),V) / F(denom(rhs(eq16)),V);

r(t) = d[vol]/(L*tau/V+KUS*V)

 

Download divid.mw

[edited] This kind of question has been asked several times before -- enough that people have written short re-usable procedures for it. So, you could use A.Jakubi's dividend procedure (and optionally pass it a sorting post-processor similar to what I used above).

restart;

dividend := proc(ex,f:=numer(ex),t:=[expand,expand])
  local n,d;
  n := t[1](numer(ex)/f):
  d := t[2](denom(ex)/f);
  n/d;
end proc:

eq16 := r(t) = d[vol]*V/(KUS*V^2+L*tau);

r(t) = d[vol]*V/(KUS*V^2+L*tau)

dividend(rhs(eq16), V);

d[vol]/(V*KUS+L*tau/V)

dividend(rhs(eq16), V,
         u->sort(expand(u),order=tdeg(V),ascending));

d[vol]/(L*tau/V+KUS*V)

 

Download dividend_ex.mw

You can check that the floats agree under rounding (to less digits), as Kitonum suggested.

Or you can increase the working precision, and then round down to current (default 10) Digits.

Or you can use testfloat to get details, and optionally use testfloat to do a coarser comparison.

restart;

a := (x^2+y^2)^3 = y^2:
b := (x^2+y^2)^(3/2) = y:

listA := evalf(sort([seq(evalf[20](solve({y >= 0, subs(x = i, a)}, {y})), i = 0.5e-1 .. .5, 0.5e-1)])):
listB := evalf(sort([seq(evalf[20](solve({y >= 0, subs(x = i, b)}, {y})), i = 0.5e-1 .. .5, 0.5e-1)])):

evalb( listA = listB );

true

restart;

a := (x^2+y^2)^3 = y^2:
b := (x^2+y^2)^(3/2) = y:

listA := sort([seq(evalf(solve({y >= 0, subs(x = i, a)}, {y})), i = 0.5e-1 .. .5, 0.5e-1)]):
listB := sort([seq(evalf(solve({y >= 0, subs(x = i, b)}, {y})), i = 0.5e-1 .. .5, 0.5e-1)]):

evalb( listA = listB );

false

zip( (a,b)->testfloat( eval(y,a), eval(y,b), 1 ), listA, listB );

[true, true, true, [false, 6., `ulps `], true, true, true, [false, 2., `ulps `], true, [false, 2., `ulps `], true, [false, 3., `ulps `], [false, 2., `ulps `], true, true, [false, 2., `ulps `], [false, 2., `ulps `], [false, 4., `ulps `], true, true]

zip( (a,b)->testfloat( eval(y,a), eval(y,b), 10 ), listA, listB );

[true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true]

evalb( `and`( zip( (a,b)->testfloat( eval(y,a), eval(y,b), 1 )=true, listA, listB )[] ) );

false

evalb( `and`( zip( (a,b)->testfloat( eval(y,a), eval(y,b), 10 )=true, listA, listB )[] ) );

true

 

Download float_comparison.mw

In your code to produce the "A" you assign to the name x, like x:=Alist[i][1] . Then you try and construct the "D" while using x as if it were an unassigned name. That won't work. There likely are other issues with the code, which is quite disorganized.

You can keep all your variables organized by using a procedure's local names for dummy variables.

You should be able to produce an unrotated "A" with just three polygons, and without any (white) occluding background, in one line of code. Even better, construct a procedure that returns the plot structure for an "A" (upright, with lower left corner at the origin).

Here's an example. (You could add a color like "gray" to the polygons, of course.) 

restart;

with(plottools):

makeA:=proc()
  uses plots;
  display(
    polygon([[0,0],[1,0],[4,8],[3,8]]),
    polygon([[3,8],[4,8],[7,0],[6,0]]),
    polygon([[2,4],[5,4],[5,3],[2,3]]),
  scaling=constrained);
end proc:

makeA();

P1:=makeA():

plots:-display(rotate(P1, 2*Pi/3));

 

 

plot_letterA.mw

interface( rtablesize = 11 ):  # or higher

w := (x, y) -> piecewise(y <= .5, -2*tanh(y-.25), .5 < y, 2*tanh(.75-y)):

Matrix( 11, 11, (i,j) -> w(i-1,j-1) );

First note that a sequence of items is different from a list of items.

restart;

f := proc()
  print(sprintf( "nargs = %a", nargs ));
  print(sprintf( "args[1] = %a", args[1] ));
  return NULL;
end proc:

f( 1,2,3,4,5 );
                          "nargs = 5"
                         "args[1] = 1"

f( [1,2,3,4,5] );
                          "nargs = 1"
                  "args[1] = [1, 2, 3, 4, 5]"

When f above is passed a sequence of 5 numbers then the number of arguments is 5, and the first argument is just 1. (nargs stands for number of arguments, and args stands for arguments.)

But when f above is passed a list with 5 entries then the number of arguments passed to f is just 1, and the first argument is the whole list.

Your maximum procedure and example call didn't work as intended because you passed it a single list (so nargs is just 1), whereas you've written it to work with a sequence.

Here are some versions that work (several other variations are possible). The first works when passed a sequence of numeric items, and the second works when passed a list of numeric items. I also give a third (which is a minor variation on the second) where the procedure's first parameter L is explicitly declared.

restart;

maximum := proc()
 local r, i;
 r := args[1];
 for i from 2 to nargs do
   if args[i] > r then
     r := args[i];
   end if;
 end do;
 r;
end proc:

maximum( 6,9,4,1 );
                               9

restart;

maximum := proc()
 local r, i, L;
 L := args[1];
 r := L[1];
 for i from 2 to nops(L) do
   if L[i] > r then
     r := L[i];
   end if;
 end do;
 r;
end proc:

maximum( [6,9,4,1] );
                               9

restart;

maximum := proc(L)
 local r, i;
 r := L[1];
 for i from 2 to nops(L) do
   if L[i] > r then
     r := L[i];
   end if;
 end do;
 r;
end proc:

maximum( [6,9,4,1] );
                               9

Secondly, note that the following two lists are NOT similar. One is a list of lists and the other is a list of scalars. The first is a list of 6 lists (each with 2 scalar entries). The second is a list with 12 scalar entries. You cannot operate on them in the exact same manner. BTW, if you changed your trailing colons to semicolons then your output wouldn't be suppressed and you'd be able to more easily see this for yourself.

LC_L:=[[0,0],[0,h],[w,h],[w,w],[b,w],[b,0]];

    LC_L := [[0, 0], [0, h], [w, h], [w, w], [b, w], [b, 0]]

LC_L2:=[(0,0),(0,h),(w,h),(w,w),(b,w),(b,0)];

         LC_L2 := [0, 0, 0, h, w, h, w, w, b, w, b, 0]

Note that round brackets do not produce a structure in any way similar to how square brackets or squiggly brackets do. Square brackets [...] construct a list. Squiggly brackets {...} construct a set. But round brackets serve as delimiters, eg. to mark the scope of arithmetic operations. Round brackets do not construct separate substructures which are subsequently accessible separately using the op command or indexed referencing.

note: round brackets also serve to denote a function call, eg f(x,y), but that's not really relevant to the points made here.

[edited]
A few more subtleties worth knowing: You can index into an expression sequence. But as soon as you pass the expression sequence to a procedure then within that procedure it (normally) gets treated as separate arguments instead of as a unified sequence.

restart;

# A list is stored as a container which could hold an
# expression sequence (amogst other things).
expr := a,b,c,d,e:

expr;
                         a, b, c, d, e

L := [expr]:

L;
                        [a, b, c, d, e]

op(L);
                         a, b, c, d, e

# `op` sees six arguments here, the 3, followed by `a`, and `b`, etc.
op(3, expr);
Error, invalid input: op expects 1 or 2 arguments, but received 6

op(3, L);
                               c

# We can unambiguously index (reference) into the sequence.
expr[2];
                               b

L[2];
                               b

# But once passed to procedure `G`, the sequence is seen as multiple
# arguments, rather than as a single argument.
G := proc(U)
  U[2];
end proc:

# Parameter `U` of `G` takes on the value of the first argument passed.
# Here, that first argument is the list `L`.
G( L );
                               b

# Parameter `U` of `G` takes on the value of the first argument passed.
# Here, that first argument is just `a` and not the whole sequence.
G( expr );
                              a[2]

The largest known Mersenne primes have been found using a program named Prime95. It has been highly optimized over the years, with multicore support and use of AVX2 and FMA3 chipset extensions.

The 50th known Mersenne prime (2^77232917-1) took the Prime95 software 6 days to find on a quad-core Intel i5-6600 @3.3Ghz.

I suspect that one important question is whether isprime knows that your trying to do a primality test of a number with the particular form (ie, 2^p-1), and whether it would use a special approach if it did. So Maple might not end up doing something like this and, even if it did, the implementation would very likely not be nearly as well optimized as Prime95. So your machine might well take orders of magnitude more resources than did the machine which actually found the 50th known Mersenne prime.

It looks like you're using Maple 16.00.
 

Here's something...

restart

with(plots):

c := 10^0:

L := -10^2:

V := -(1/2)*signum(1+2*c*x)*abs(x)*(c*x^2+x-2*L)/(1+2*c*x)^2:

Student:-Calculus1:-ExtremePoints(V, x = 0 .. 1, numeric);

[0., .5050829867, 1.]

simplify(diff(V,x)) assuming x>=0, x<=1;
select(u->is(u>=0 and u<=1), evalc([solve(%)]));
evalf(%);

-(1/2)*(3*x^2+2*x^3-398*x+200)/(1+2*x)^3

[-(1/6)*2397^(1/2)*sin((1/3)*arctan((2/9)*579^(1/2))+(1/6)*Pi)-1/2+(1/6)*3^(1/2)*2397^(1/2)*sin(-(1/3)*arctan((2/9)*579^(1/2))+(1/3)*Pi)]

[.505082985]

kernelopts(version);

`Maple 16.00, X86 64 LINUX, Mar 3 2012, Build ID 732982`

 

Download 2_2.mw

Since you're interested in performance you should not be calling diff each time you call w2.

Try this,

w2:=unapply(diff(v2(x,y),x)-diff(u2(x,y),y),
            [x,y]);

 

If either diff(v2(x,y),x) or diff(u2(x,y),y) return something which contains something other than what you called an "explicit function" or an unevaluated diff call, etc, then show us the full example.

[edited] Below is an example where both u2(x,y) and diff(u2(x,y),y) return unevaluated for nonnumeric arguments x and y. In this case the above approach will not work. So the question becomes, how to evaluate the derivative of this "black box" (non-explicit) function numerically? One way is to use the fdiff command, and variant on that idea is to use evalf(D[..](...)) which has a kernel hook into fdiff.

Since you want the whole call to w2 to run under fast evalhf mode (since v2 might still run under it) then I'll use eval to allow the call to evalf(D[..](...)) to be called from within evalhf.

restart;

v2 := proc(x, y) sin(x)*cos(y) end proc:

u2 := proc(x, y)
        if not (x::numeric and y::numeric) then
          return 'procname'(args);
        end if;
        3*sin(x)-cos(y);
      end proc:

diff(u2(x,y),y); # returns unevaluated, so cannot use original approach

diff(u2(x, y), y)

w2:=unapply(diff(v2(x,y),x)-'eval'('evalf'(D[2](u2)(x,y))),
            [x,y]):

CodeTools:-Usage(evalhf(add(add(abs(w2(i, j)), i = 0 .. 10), j = 0 .. 10)));

memory used=8.34MiB, alloc change=1.00MiB, cpu time=45.00ms, real time=45.00ms, gc time=0ns

93.6605340977501584

 

Download evalf_D.mw

If either of your u2 or v2 is the solution to a numeric pde problem (ie. coming from pdsolve numeric) then you might want to consider that the derivatives in question might be obtainable from the original differential equations (as some function of their terms, all supplied by the pdsolve solution). If so then that would likely be more numerically robust than using fdiff which is a numeric differentiation command and brings with it yet another layer of rough discretization.

The curves meet where the two expressions attain the same value, ie. where the difference of the two expressions is zero. So you're looking for the roots of BesselJ(0,x) - BesselJ(1,x) .

Student:-Calculus1:-Roots(BesselJ(0, x)-BesselJ(1, x), x=0..10, numeric);

            [1.434695651, 4.680102554, 7.836002335]

Here's a slightly more detailed look at it, considering just the first of the three roots:

ans := Student:-Calculus1:-Roots(BesselJ(0, x)-BesselJ(1, x), x=0..10, numeric);

         ans := [1.434695651, 4.680102554, 7.836002335]

x=ans[1];
                        x = 1.434695651

# Evaluate both expressions at that point. They are effectively equal.

eval( [ BesselJ(0, x), BesselJ(1, x)] , x=ans[1] );

                  [0.5479464494, 0.5479464495]

# Evaluate their difference at that point. It is effectively zero.

eval( BesselJ(0, x)-BesselJ(1, x), x=ans[1] );

                                 -10
                            -1 10   

That is not a bug.

By default LinearAlgebra arithmetic commands operate in hardware double precision. The computation of the 2nd entry in your result is similar to this:

restart;

HFloat(1.99987);

                        1.99987000000000
HFloat(-2.);

                              -2.

evalhf(HFloat(1.99987) + HFloat(-2.));

                    -0.000129999999999963478

Vector([%], datatype=float[8]);

                [-0.000129999999999963]

Note that neither 1.99987 nor -0.00013 is expressible as a hardware double precision float with "infinite" precision.

lprint( HFloat(1.99987) );
  HFloat(1.99987000000000004)

lprint( HFloat(-0.00013) );
  HFloat(-0.129999999999999989e-3)

If you want numeric LinearAlgebra operating with base-10 floats and "machine epsilon" like 10^(-Digits), rather with base-2 hardware floats with "machine epsilon" like evalhf(DBL_EPSILON), then you can force LinearAlgebra to operate in software floats rather than its default of hardware double precision floats.

You can do that by setting the UseHardwareFloats Maple envionment variable to false, or by raising Digits higher than 15 and leaving UseHardwareFloats at its default setting of deduced.

restart:

UseHardwareFloats; # default is `deduced`

deduced

with(LinearAlgebra):
Ax := <1.00004, 1.99987, -0.12e-3>:
b := <1., 2., 0.>:

Add(Ax, b, 1, -1);

_rtable[18446884781111771974]

UseHardwareFloats:=false:

Add(Ax, b, 1, -1);

_rtable[18446884781111774014]

 

Download LA_Add_hw_sw.mw

What do you mean by "search interval"?

restart;
solve( {sin(sqrt(x))}, {x}, allsolutions );

                         /      2    2\ 
                        { x = Pi  _Z1  }
                         \            / 

about(_Z1);
  Originally _Z1, renamed _Z1~:
    is assumed to be: integer

restart;
Student:-Calculus1:-Roots( sin(sqrt(x)), x=0 .. 200 );

                 [     2      2      2       2]
                 [0, Pi , 4 Pi , 9 Pi , 16 Pi ]

restart;
Student:-Calculus1:-Roots( sin(sqrt(x)), x=80 .. 700 );

        [    2       2       2       2       2       2]
        [9 Pi , 16 Pi , 25 Pi , 36 Pi , 49 Pi , 64 Pi ]

Using an image for a plotted 3D surface can now be done using an option (ie. the image option).

But before that was possible, it could still be done more "manually".  See this old Post on mine, which has some examples with World map projections and even a few links for images of other planets surfaces. The technique was subsequently incorprated as new visualization functionality in Maple 18. Eg, with a user-friendly syntax,

plot3d(-1, t = -Pi .. Pi, p = 0 .. Pi,
       image = "http://www.maplesoft.com/data/examples/jpg/world.200412.400x200.jpg",
       lightmodel = none, coords = spherical, axes = none);

It is sad that the Help pages for the WorldMap functionality (new in Maple 2017) don't include any clear examples of how to display efficiently the DataSets,Builtin,WorldMap information on a 3D globe or sphere.

It is also sad that the Online Help doesn't have the new material from Maple 2016 and 2017 (it's 2018, folks).

 

In your Document I get right-click context menus to work OK on 2D Output, whether in a paragraph (Document Block) or an Execution Group.

But what doesn't work in your original is right-click context menus on 2D Input, whether in a paragraph (Document Block) or an Execution Group.

If I use the mouse to toggle "Show Command" on an offending Document Block, to reveal what the context-menu system inserted, then I see the expected command but with the placeholder %EXPR instead of the argument. Eg, solve(%EXPR, [b]) instead of solve(a=b+5, [b]) .

If I do all of the following then the right-click context menus work OK on 2D Input in your Document, for both new paragraphs and your two original paragraphs (Document Blocks):

1) Get rid of all the User Profile stuff in the Startup code region. I mean delete it all, not comment it out.
2) Get rid of the displayunit and outputmask attributes in the actual .mw file (manually).
3) Get rid of the UnitFormat-DefaultUnitFormat and NumericFormat-OutputMask in the View-Properties blob of the actual .mw file (manually).

So it looks like the combination of custom unit-formatting and custom numeric-formatting, loaded by the User Profile, are breaking context-menu actions for 2D Input (only).

You may have to give up on unit- and numeric formatting customizations from a User Profile, to get right-click context menus to work on 2D Input in new Documents with your Maple 2017.3. But context-menus might work regardless, on 2D Output.

note: You last saved your Document in Maple 2017.3, and I was working in Maple 2017.2 but it sounds like I was seeing the same issues as you were with your Document.

[edit] This may work as an alternative for a so-called User Profile. Take the key lines you want from your earlier Startup Code region, such as calls to the interface command for numeric formatting or the imaginary unit, and put them in the Startup Code region of a new blank Document. Don't copy over any of the lines mentioning User Profiles. Now save that blank Document somewhere permanent on your machine, and set it as the default "home" file, eg from the menubar Tools -> Options -> Interface -> Open worksheet at startup . If you do this then every time you open a New Document you should get a fresh copy of that customized template, with the entries in its Startup Code Region.

The f is not superfluous. You need to specify some numeric format (d, o, x, e, f, or g) in order to use z or Z.  The z or Z is just a modifier of that numeric format.

The reason that the "c" is in red in the Help page, where it describes "zc or Z", is that it is a placeholder (ie. it stands for the character you choose to type there).

One difference between "zc" and "Z" is that when modifying by Z you won't get the real and imaginary parts separated by anything other than the +/- sign of the imaginary component. The "zc" modifier allows you to separate them. Another difference it that using Z puts the imaginary unit at the end, whereas with zc it doesn't (but you can add it by hand).

If you only want to print both real and imaginary components together (separate by a single space, say), then it's as convenient to use the zc modifier as it is pass two values using Re and Im.

I mentioned that the Z modifier prints the current symbol for the imaginary unit. If you use the zc modifier and want to insert the imaginary unit in manually then you could even make it robust enough to also respect that setting.

restart;

interface(imaginaryunit=I): # default

ee:=1234.5678-7733.6644*I;

1234.5678-7733.6644*I

printf("%z gI", ee);

1234.57 -7733.66I

printf("%Zg", ee);

1234.57-7733.66I

 

 

printf("%z fI", ee);

1234.567800 -7733.664400I

printf("%Zf", ee);

1234.567800-7733.664400I

 

 

printf("%z eI", ee);

1.234568e+03 -7.733664e+03I

printf("%Ze", ee);

1.234568e+03-7.733664e+03I

 

 

printf("%6.5z f\n", ee);

1234.56780 -7733.66440

printf("%6.5f %6.5f\n", Re(ee), Im(ee));

1234.56780 -7733.66440

 

interface(imaginaryunit=j):

ee;

1234.5678-7733.6644*I

printf("%Zf", ee);

1234.567800-7733.664400j

 

printf("%6.5z f%a\n", ee, sqrt(-1));

1234.56780 -7733.66440j

printf("%6.5z f%a\n", ee, interface(imaginaryunit));

1234.56780 -7733.66440j

 

interface(imaginaryunit=I):

printf("%6.5z f%a\n", ee, interface(imaginaryunit));

1234.56780 -7733.66440I

printf("%Zf", ee);

1234.567800-7733.664400I

 

Download printf_zcZ.mw

 

 

 

How about moving the line print('J'[n]=J[n]) to after the call to Jacobian, which constructs the Matrix and assigns it to J[n]?

restart

" iter:=5;  f[1](x,y):=3 x^(2)-y^(2);  f[2](x,y):=3 x^()*y^(2)-x^(3)-1;"

5

proc (x, y) options operator, arrow; 3*x^2-y^2 end proc

proc (x, y) options operator, arrow; 3*x*y^2-x^3-1 end proc

var := x, y

x, y

pointt := [x[n], y[n]]

[x[n], y[n]]

NULL

NULL

``

x[0] := 1; y[0] := 1

1

1

for n from 0 to iter do print('f1' = f[1](x[n], y[n]), 'f[2]' = f[2](x[n], y[n])); J[n] := Student[MultivariateCalculus][Jacobian]([f[1](x, y), f[2](x, y)], [var] = pointt, output = matrix); print('J'*[n] = J[n]); sol[n] := eval((Vector(2, {(1) = x[n], (2) = y[n]}))-1/J[n].(Vector(2, {(1) = f[1](x[n], y[n]), (2) = f[2](x[n], y[n])}))); x[n+1] := evalf(sol[n][1]); y[n+1] := evalf(sol[n][2]); print(x[n+1], y[n+1]) end do

f1 = 2, f[2] = 1

J*[0] = _rtable[18446884754823685054]

.6111111111, .8333333333

f1 = .4259259256, f[2] = 0.44924554e-1

J*[1] = _rtable[18446884754823680958]

.503659080870043, .852494422128773

f1 = 0.342706694679006e-1, f[2] = -0.296666580332421e-1

J*[2] = _rtable[18446884754823678422]

.499964121072352, .866045636385908

f1 = -0.142677224123089e-3, f[2] = -1.25763981939642*10^(-6)

J*[3] = _rtable[18446884754823674566]

.500000000014920, .866025401817003

f1 = 3.45245787514159*10^(-9), f[2] = -5.08916708774620*10^(-9)

J*[4] = _rtable[18446884754823670710]

.500000000000000, .866025403784439

f1 = 1.11022302462516*10^(-16), f[2] = -2.22044604925031*10^(-16)

J*[5] = Matrix(%id = 18446884754823666854)

HFloat(0.5), HFloat(0.8660254037844387)

``

``

Download q1nwtnnonlinearsys_1.mw

First 183 184 185 186 187 188 189 Last Page 185 of 336