Carl Love

Carl Love

28025 Reputation

25 Badges

12 years, 313 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

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).

 

Although pasting works sometimes, your learning is much improved by explicitly typing examples even when pasting works. Perhaps as you typed, you'd realize that what's pasted above is nonsense because it's missing both the multiplication and exponentiation operators.

Obtaining non-real floats (i.e., with nonzero imaginary parts) when you expect only reals is not a sign of an error or a cause for concern. There are two reasons that you might get these: round-off error or ignoring symbolic constants.

Round-off error: The decimal approximation of some real constants often necessarily involves detours into the complex numbers.  For example, ((-1)^(1/3))^3 is clearly -1, but if I include a decimal point:

((-1)^(1/3.))^3;
      -.9999999999+3.865045973*10^(-10)*I

The small magnitude of that imaginary part is a sign that that imaginary part is simply due to round-off error and can be ignored. My rule of thumb for this is if you're expecting a real number and

Im(z) <> 0 and Re(z) <> 0 and ilog10(Im(z)) < 3-Digits and ilog10(Im(z)) - ilog10(Re(z)) < 3-Digits

and these relations remain true when you increase Digits, then the imaginary part can be ignored and/or discarded. (It's just a rule of thumb: there are a few situations where it doesn't work.) There are easy ways to remove these imaginary parts which we can go over when the need arises. At this developmental phase of your coding project, you should acknowledge their presence and mentally ignore them without explicitly removing them.

Ignoring symbolic constants: I don't know where you're at in your mathematical education, but you're probably at least vaguely aware of the following (each of these is either covered in a precalculus course  or a first course in differential equations).

  1. The general solution of a homogenoous linear ODE with constant coefficients involves finding the roots of a univariate polynomial with the same coefficients and with the degree of each term being the differential order of the corresponding term in the ODE.
  2. That if the roots r[k] of that polynomial are distinct then the general solution is y(t) = Sum(C[k]*exp(r[k]*t), k= 1..n), where n is the differential order of the ODE and the C[k] are constants which can't be specified until initial or boundary conditions are provided. (Only minor adjustments to that are needed when the roots aren't distinct.)
  3. The roots of a polynomial are in general non-real even if the coefficients are real. 
  4. If the coefficients are real, then any non-real roots occur in complex-conjugate pairs (e.g., x+I*y and x-I*y).
  5. For real texp(I*t) = cos(t) + I*sin(t). (That's still true if t is non-real, but we only need to consider the real case here.)
  6. The sum of any two solutions of a linear homogenous ODE is also a solution.
  7. There is no formula or algorithm for finding the exact roots of a general polynomial of degree greater than 4. (This is one reason why you often see polynomial roots expressed in Maple's RootOf form.) Indeed, it's been proven that no such formula or algorithm will ever be found. There are algorithms to find decimal approximations of the roots to arbitrary accuracy.

If you put these facts together, the result is that the real solution to a real-coefficient ODE without initial/boundary values might be expressed as y(t) = sum(C[k]*exp(r[k]*t), k= 1..n) where some of the r[k] are manifestly non-real (not of the "round-off error" variety described above), and in complex-conjugate pairs. When the values of the C[k] are determined, some of them will also be non-real, but they will be such that the overall solution is real valued, upto the limitations of round-off error.

Numeric BVP solver: I think that your problem can be solved by Maple's numeric BVP solver, which is accessed by 

Sol:= dsolve(
   {
some ODEs,
    some boundary conditions (BCs) expressed at the two endpoints of a closed real interval of the independent variable
   },
   numeric,
   
possibly some other numeric options (as needed)
);

You should at least try. Symbolic parameters can be included as long as you have exactly 1 BC for every symbolic parameter plus the total differential order of the ODEs. The BCs can be quite complicated as long as the only explicit values of the independent variable that are used are the two endpoints of the interval. You can use a parameter for the eigenvalue. The solution of the system will include real numeric values for all parameters. This technique completely avoids non-real values. See help page ?dsolve,numeric,bvp.

Lengthy output messages: The message "Length of output exceeds ..." is not an error message nor even an indication that there was any snag in computing the solution that you requested. The message is merely there to warn you that the computed output is so long that you may not want it displayed on your screen, or at least you may not want it in prettyprinted format. The easiest way to look at the solution is

lprint(%);

which'll give it in plaintext form. However, this is always a good time to ask yourself "Is it reallly worth looking at? Can I analyze it in some other way (numeric evaluation, plot, breaking it into parts, etc.)?"

 

For close to 0, is a close approximation of sin(x). You can see that from its Maclaurin series. So all you're seeing are shifted approximations of sqrt(10).

My workaround is similar to Tom's, and additionally it'll continue to work if and when Maplesoft changes the parameters, regardless of how they're changed (even if it's more complicated than just switching their order):

MyGammaD:= proc(k, theta)
uses St= Statistics;
local a,b, GD:= 'GammaDistribution'(a,b);
   eval(GD, solve({St:-Mean(GD) = k*theta, St:-Variance(GD) = k*theta^2}, {a,b}))
end proc:

I'm not claiming to know whether Maplesoft has ever changed the parameters. The above works regardless of changes.

Even more bothersome to me than this thing with GammaDistribution is that there doesn't seem to be standardization across mathematical software and literature as to whether the second parameter of Normal is sigma or sigma^2.

I don't know what you're trying to achieve with the unevaluation quotes, but they have no effect whatsoever when used on a procedure parameter. Other than that, it looks like you're trying to match expressions that are invariant under op. Any conditional can be made a type by using type satisfies. So, the type is satifies(x-> x = op(x)). This won't error out if op(x) is multi-term, so there's no need for square brackets.

With tools like unapply and ->, there's no practical limit to the level of "meta differential operators" that you can create, each capable of being iterated with @@. In this example, I define an operator Newton which takes an "operator template" (f,c) and returns the differential operator that constructs the Newton's method iteration operator for solving f(x) = c. This final operator can be iterated with @@. Each level of these operators is only one line of code.

Newton:= proc(f,c) local x; subs(_c= c, _c-> unapply(x - f(x,_c)/D[1](f)(x), x)) end proc:
Sqrt:= Newton(((x,c)-> x^2-c), c): 
(Sqrt(2)@@4)(1.5); #1.5 is arbitrary starting point.
                          1.414213562
%^2;
                          1.999999999

The output produced by Newton itself is very ugly to look at. But the result of Sqrt(2) is what you'd expect.

Regardless of whether its number of elements is even or odd, the median[1] of a nonempty[2] sorted list X equals

add(X[[floor,ceil]((nops(X)+1)/2)])/2

That idea is used in the following, where I also address an issue that you'll encounter when your input has symbolic real constants that you may not want to clobber with evalf in the output:

#Note that the default sort order doesn't respect the usual order
#of real numbers when the numbers are symbolic:
sin~([$1..4]): (sort, evalf)(%);
     [sin(1), sin(2), sin(3), sin(4)], [0.8414709848, 0.9092974268, 0.1411200081, -0.7568024953]

#We can get the desired order by sorting with evalf:
(x::list)-> ((S,n)-> (S, add(S[[floor,ceil]((n+1)/2)])/2, add(S)/n))(sort(evalf(x)), nops(x)):
   
%(sin~([$1..4]));
     [-0.7568024953, 0.1411200081, 0.8414709848, 0.9092974268], 0.4912954964, 0.2837714810

#We can do the same thing plus not clobber the input
#by using sort(..., output= permutation):
(x::list)-> 
   ((P,n)-> (x[P], add(x[P[[floor,ceil]((n+1)/2)]])/2, add(x)/n))
      (sort(evalf(x), 'output'= 'permutation'), nops(x))
:   
%(sin~([1,2,3,4]));
                                      1          1         
    [sin(4), sin(3), sin(1), sin(2)], - sin(3) + - sin(1), 
                                      2          2         

      1          1          1          1       
      - sin(1) + - sin(2) + - sin(3) + - sin(4)
      4          4          4          4 
      
evalf([%]);
     [[-0.7568024953, 0.1411200081, 0.8414709848, 0.9092974268], 0.4912954964, 0.2837714811]

[1]I mean median in the sense that it's taught in elementary math, not the more-technical sense discussed by @sand15 .

[2]And do we really want to give serious consideration to defining the median of an empty list?

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