Carl Love

Carl Love

28065 Reputation

25 Badges

13 years, 25 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

It is done by Typesetting:-mn(gcd)(x,y). This command is undocumented; however, if you enter exports(Typesetting), then in addition to the documented commands, you'll see many whose names begin with m. These are all commands for doing this sort of thing.

As Acer says, the extra prompts make no difference to the meaning of the code. On the other hand, with respect to editing the code, I find them quite a nuisance (you can't backspace over them) as well as ugly. When they appear in my code, I try to delete them as soon as possible. So, yes, I'd say that there's an advantage to not having extra prompts.

What purpose do they serve? Is there some reason that Maplesoft doesn't remove the prompt from the bottom group when execution groups are merged (aka joined)? Or is it just an oversight?

And this seems as good a place as any to remind people once again that restart commands should always be placed in their own execution group. Not doing do only occasionally causes problems, but those problems can be quite weird.

The error message indicates that your variable O[L] (or O__L) does not have units. The absence of units is considered the same thing as having 1 as the unit, as indicated by the message. If you replace O[L] by O[L]*Unit(kN), it should work. I can't test it unless you supply your actual code.

Your command print(whattype(df)) is leading you to believe that df has been converted from a Matrix into a procedure. In fact, it's still a Matrix. It's an idiosyncracy of print combined with the fact that the symbol Matrix is both a type and a procedure. If you use lprint instead of print, it'll say "Matrix".

I have many more suggestions for how you can improve your code, which we can go into. 

Any of the following commands will return ln(a+b*c), which is, I think, what you're asking for.

indets(expr, function(Not(function)))[];
indets(expr, anyfunc(Not(function)))[];
indets(expr, function({`+`, `*`}))[];
indets(expr, anyfunc({`+`, `*`}))[];

You need to give deeper (i.e., more complex, more profound, more nested) examples for me to refine the above commands. You say that you want the first, but there is no second in your example, so I cannot show you how to distinguish those.

Recalling your previous Question, I need to add: If you are processing the input through InertForm:-Parse, then the above commands need to be changed because then operators will be turned into functions, so the above commands will not distinguish between sin and `%+` for example.

Among the mathematical functions, special attention needs to be given to sqrt with respect to how it behaves with respect to your type of work: sqrt(x) is either simplified (assuming that x is some expression for which that is possible) or converted to x^(1/2) (or some combination of simplification and conversion to exponent form, as with sqrt(12)). So, if you search for sqrt in some expression that you just entered, you won't find it, it becomes an operator rather than a function. This can be prevented by using the inert form %sqrt rather than sqrt. As you may have guessed by now, every function f has an inert form %f. The command value is used to convert the inert forms to active forms (as I did in my Answer to your previous Question).
 

You can ask Maple:

is(cosh(x) <= exp(x^2/2)) assuming x::real;
      true

This can mostly be handled by InertForm:-Parse followed by some filtering with indets and a special structured type (that I composed just for this situation). Here's your example done:

value(
   indets(
      InertForm:-Parse("y = s + (a/b)*log((a+c/b))"), #string!
      #selects operators with exactly two operands, both of which are names,
      #integers, or floats:
      And(specfunc({name,numeric}, {`%+`, `%*`, `%/`}), 2 &under nops) 
   )
);
                             /a  c\ 
                            { -, - }
                             \b  b/ 

Note that the argument to Parse is a string. This prevents any automatic simplifications, thereby avoiding the problem mentioned by Preben.

I said mostly because the minus sign is treated differently in Maple, even by InertForm. It is always (effectively) an unary operator when applied to a symbolic expression. So, for example:

-1 is atomic (any negative explicit integer or float is atomic).
-a is stored as (-1)*aeven by InertForm.
a-b is stored as `%+`(a, (-1)*b), where `%+` is the inert operator in prefix form.
-a-b is stored as `%+`((-1)*a, (-1)*b).
-a-2*b is stored as `%+`((-1)*a, (-1)*`%*`(2, b)).
-2*a-2*b is stored as `%+`(`%*`(-2, a), (-1)*`%*`(2, b)).

As far as I can tell so far, there is no weirdness with the handling of fractions or complex numbers: They are deconstructed into `%operator`(operands) form in the way that you probably expect. You didn't mention exponentiation, so I didn't investigate that.

If you consider the last 2 of my 6 examples above, I think that you can imagine some of the problems and ambiguities that can arise. So, if we are to proceed, then you need to show many examples with minus signs and say for each how you would like it deconstructed.

I see that you've loaded the command codegen[makeproc], but you haven't used it in the code that you're asking about. Using it is the way to solve your problem. Here's my example: I generate two long, random polynomials in variables x, y, z (to serve as examples in place of your two expressions, which I can't directly see). Then I translate them to optimized Fortran two ways: the first being the way you used, and the second by first applying codegen[makeproc]. Execute this example and observe the differences:

(p1,p2):= 'randpoly([x,y,z], degree= 32)' $ 2;

CodeGeneration:-Fortran(p1, optimize, resultname= dp1);
CodeGeneration:-Fortran(p2, optimize, resultname= dp2);

CodeGeneration:-Fortran(codegen[makeproc](p1, [x,y,z]), optimize, resultname= dp1);
CodeGeneration:-Fortran(codegen[makeproc](p2, [x,y,z]), optimize, resultname= dp2);

The second pair of Fortran commands might re-use some of the variable names, but, if it does, it won't cause any problems because those names will be local to the Fortran functions in which they occur.

It is not fixed.

A crude way that I work around such bugs is to replace arbitrary symbols, such as your a, with arbitrary (or randomly chosen) complex transcendental constants. For example, I replaced a with sin(23+I) (as well as several other such values). In all cases, the answer 2*Pi is obtained. Of course, this result cannot be accepted with mathematical certainty, but it's much safer and more certain than replacing a with a rational.

Regarding the OP's followup problem of finding integer solutions to x = a + sqrt(b*x+c), I here present two algorithms that are much faster than that posted by Kitonum. The first of these is about 35 times faster. It works with the symbolic solution of the equation and avoids repeated calls to solve. The procedure for this method can almost be automatically generated from the symbolic solution. This method could also be applied to the OP's first problem. My second algorithm is about 700 times faster than Kitonum's. To generate it, I analyzed the symbolic solution (by hand) and reduced the sufficient conditions to a few simple inequalities.
 

Author: Carl Love <carl.j.love@gmail.com>

17 Sep 2018

restart:

The Problem: Find triples of distinct integers (A, B, C) such that the following algebraic equation has two distinct integer solutions for x, and B and C are relatively prime.

Eq:= x = A + %sqrt(B*x+C);

x = A+%sqrt(B*x+C)

Search space:
(In order to get a reasonably measurable execution time for my fastest method, I need to expand the search space.)

Range:= [(-20..30) $ 3]:

Kitonum's method:

In order to compare solutions, I made the solutions for each triple a set rather than a list, and the overall list a set.

Eq1:= value(eval(Eq, [A,B,C]=~ [a,b,c])):

gc(): st:= time():
k := 0:
for a from op([1,1], Range) to op([1,2], Range) do
for b from op([2,1], Range) to op([2,2], Range) do
for c from op([3,1], Range) to op([3,2], Range) do
if igcd(b, c) = 1 and nops({a, b, c}) = 3 then X := {solve(Eq1)};
if nops(X)=2 and type(X[1], integer) and type(X[2], integer) then k := k+1; L[k] := [a, b, c, X] end if end if end do end do end do;
L1:= convert(L, set):
gc(): T1:= time()-st;

261.141

Carl's 1st method:

The key to the speed of this method is that solve is only used once, symbolically.  

SymSol:= subsindets(
   {solve(value(Eq), x, 'explicit')},
   'anything'^'identical'(1/2),
   s-> %sqrt(op(1,s))
);

{A+(1/2)*B-(1/2)*%sqrt(4*A*B+B^2+4*C), A+(1/2)*B+(1/2)*%sqrt(4*A*B+B^2+4*C)}

Using this symbolic solution, the procedure that checks whether an integer triple is a solution can be almost automatically generated.

#This procedure uses syntactic features introduced in Maple 2018.
#It won't parse in earlier Maple.
#There's an equivalent procedure for earlier Maple immediately below.
CheckSym:= proc(abc)
local a,b,c,S,V;
   if
      nops({((a,b,c):= seq(abc))}) = 3 and igcd(b,c) = 1                      and  
      (V:= value((S:= eval([Eq, SymSol], [A,B,C]=~ [a,b,c]))[2]))[1]::integer and
      nops(V) = 2 and andmap(evalb, value(map2(eval, S[1], x=~ V)))
   then
      [a, b, c, V]
   else #return NULL
   fi
end proc:

#equivalent procedure for Maple 2017 or older:
CheckSym2017:= proc(abc)
local a,b,c,S,V;
   (a,b,c):= seq(abc);
   if nops({a,b,c}) = 3 and igcd(b,c) = 1 then
      S:= eval([Eq, SymSol], [A,B,C]=~ [a,b,c]);
      V:= value(S[2]);
      if
         V[1]::integer and nops(V) = 2                 and
         andmap(evalb, value(map2(eval, S[1], x=~ V)))
      then
         [a, b, c, V]
      else
      fi
   else
   fi
end proc:

GenerateSols:= proc(Check::procedure, Range::list(range))
local abc;
   {seq(Check(abc), abc= Iterator:-CartesianProduct([`$`]~(Range)[]))}
end proc:
          

#This is just a procedure to run time tests.
#It's not related to the algebra problem at hand.
UsageWithGC:= proc(e::uneval)
   gc();
   CodeTools:-Usage(proc() local r:= eval(e); gc(); r end proc(), args[2..])
end proc:   
   

#Tests:
(L2A,T2A,L2B,T2B):= map(
   P-> UsageWithGC(GenerateSols(P, Range), 'output'= ['output', 'cputime']),
   [CheckSym, CheckSym2017]
)[]:

memory used=0.70GiB, alloc change=1.04MiB, cpu time=7.11s, real time=6.48s, gc time=875.00ms

memory used=0.70GiB, alloc change=1.04MiB, cpu time=6.95s, real time=6.36s, gc time=828.12ms

Carl's 2nd method:

A careful analysis of the symbolic solution and its substitution back into the original equation shows that most of the conditions that guarantee the type of solution required can be reduced to simple inequalities. This gives another big speed improvement.

#This procedure uses syntactic features introduced in Maple 2018.
#It won't parse in older Maple.
#There's an equivalent procedure for older Maple below.
CheckFast:= proc(abc)
local a,b,c,d,s;
   (a,b,c):= seq(abc);
   if
      b > 0 and a*b <= -c              and #Both roots satisfy x = a+sqrt(b*x+c).
      a <> b and a <> c and b <> c     and #user-specified requirement
      (d:= 4*(a*b+c)+b^2) > 0          and #There are two real roots.
      igcd(b,c) = 1                    and #user-specified requirement
      (s:= isqrt(d))^2 = d                 #Both roots are integer.
   then
      [a, b, c, {a+Scale2(b+s,-1), a+Scale2(b-s,-1)}] #return input and roots
   else #return NULL
   fi
end proc:

#Equivalent procedure for Maple 2017 or older:
CheckFast2017:= proc(abc)
local a,b,c,d,s;
   (a,b,c):= seq(abc);
   if b > 0 and a*b <= -c and a <> b and a <> c and b <> c then
      d:= 4*(a*b+c)+b^2;
      if d > 0 and igcd(b,c) = 1 then
         s:= isqrt(d);
         if s^2 = d then
            [a, b, c, {a+Scale2(b+s,-1), a+Scale2(b-s,-1)}]
         else
         fi
      else
      fi
   else
   fi
end proc:

#Tests:
(L3A,T3A,L3B,T3B):= map(
   P-> UsageWithGC(
      GenerateSols(P, subsop([2,1]= max(1, op([2,1], Range)), Range)),
      'output'= ['output', 'cputime']
   ),
   [CheckFast, CheckFast2017]
)[]:

memory used=4.68KiB, alloc change=0 bytes, cpu time=359.00ms, real time=237.00ms, gc time=296.88ms

memory used=4.68KiB, alloc change=0 bytes, cpu time=344.00ms, real time=235.00ms, gc time=296.88ms

#Verify that all the solutions are the same:
nops({L1, L2A, L2B, L3A, L3B});

1

nops(L1);

454

#computation speed increase factors:
T1 /~ [T2A, T2B, T3A, T3B], [T2A,T2B] /~ [T3A,T3B];

[36.73385848, 37.55803250, 727.4122562, 759.1308139], [19.80222841, 20.21220930]

#random sample of solutions:
combinat:-randcomb(L3B, min(32, nops(L3B)));  

{[-12, 1, 12, {-12, -11}], [-8, 29, 22, {6, 7}], [-6, 7, 30, {-3, -2}], [-6, 19, 24, {3, 4}], [-6, 25, 24, {1, 12}], [-6, 26, -9, {5, 9}], [-5, 19, 7, {3, 6}], [-5, 23, 3, {2, 11}], [-5, 23, 13, {1, 12}], [-5, 28, 25, {0, 18}], [-5, 30, -11, {2, 18}], [-4, 7, 18, {-2, 1}], [-4, 13, 30, {-2, 7}], [-4, 16, 1, {3, 5}], [-4, 19, 16, {0, 11}], [-4, 20, -11, {3, 9}], [-4, 23, -20, {3, 12}], [-4, 24, 1, {1, 15}], [-3, 11, 5, {1, 4}], [-3, 13, 9, {0, 7}], [-3, 20, 9, {0, 14}], [-3, 23, -7, {1, 16}], [-3, 25, -9, {1, 18}], [-2, 12, 13, {-1, 9}], [-2, 14, -5, {1, 9}], [-2, 29, 30, {-1, 26}], [-1, 7, -5, {2, 3}], [-1, 11, -17, {3, 6}], [-1, 15, -11, {1, 12}], [0, 8, -15, {3, 5}], [1, 10, -19, {2, 10}], [13, 1, -13, {13, 14}]}

 

 


 

Download IntegerSols.mw

We can generalize the idea to any container (Matrix, Vector, Array (of any number of dimensions), table, list, set) that contains a single element. The following procedure will either

  • unpack its argument if its argument is a container with a single element, or
  • simply return its argument otherwise.

And it applies itself recursively so that nested singleton containers are unpacked.

UnpackSingletonContainers:= proc(x, $)
   if x::'indexable' and not x::'string' and numelems(x) = 1 then
      thisproc(`if`(x::{table, rtable}, entries, op)(x))
   else
      x
   fi
end proc:

 

The warning says that the system is inconsistent, not that it's consistent.

In your dsolve line, you wrote u1(x,t) twice. You meant for one of those to be u2(x,t).

If by i you mean sqrt(-1), then you should begin your code with interface(imaginaryunit= i);

Your code works for me (using Maple 2018 2D Input, same as you). But I think that Acer figured out what your problem is.

Nonetheless, the following is much more robust, seems more scientific to me, and will be better in the long run:

E__n:= unapply(
   combine(
      evalindets(
         'h'*'c'*'R'[infinity],
         'name',
         x-> eval(
            value*'Unit'(units), 
            [ScientificConstants:-GetConstant(x)][2..]
         )
      )/(-n^2),
      units
   ), n
 );

eV:= unapply(x*convert(Unit('J'), units, Unit('eV')), x);

Now the code contains no "magic numbers"; they've all been replaced by internationally recognized symbols which are automatically replaced by the appropriate numbers as needed. Also, all the unit simplification is handled.

[Edit: See much-simplified procedure below.]

I have four suggestions which unfortunately I can't test because you haven't uploaded your code. (I see that you've posted it now.) Any one of these suggestions may be enough to get the integration working.

1. It is never necessary to use and (or three-part inequalities) in the conditions of piecewise in the way that you have done because the conditions are processed left-to-right and so every condition can tacitly use the assumption that all conditions to the left of it are false. (I am sure of what I've written so far.) I have a suspicion (although I'm far from sure about it) that these ands complicate the processing of piecewise expressions by int.

2. Maple's numeric integration by default will not return a result unless it can guarantee that that result is accurate to Digits significant digits. Sometimes it struggles to figure out how many extra digits it needs to use in internal computations (sometimes called guard digits) to guarantee that. You can guide it (and override the default) by using the options epsilon and digits. See ?evalf,Int.

3. Your expressions are essentially polynomial (from what I can see from your "code picture"). I wouldn't be surprised if symbolic integration would work.

4. piecewise expressions can usually be converted to Heaviside expressions. Sometimes (maybe always) these are handled better by int.

 

diff has no problem with the name sol__i__j; however the __i and __j are not indices, and they have no connection to the loop variables and j. They do display as subscripts, but subscripts are not necessarily indices.

One way to double-index sol and sys would be sol[i,j] and sys[i,j]. There are several other ways.

Before you use assign on dsolve's output, you should check that there's actually some output. It might return NULL if it can't solve the equations. In order for me to provide more details, you'll need to post your modified code (either as plaintext or an uploaded worksheet).

 

First 163 164 165 166 167 168 169 Last Page 165 of 395