Carl Love

Carl Love

28085 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

The immediate first thing that I noticed was the same thing discussed by acer. But there are several other issues, and I decided to show you how to do a more thorough analysis of your code.

These are some steps to take that will quickly show you which line of code is taking so long. But first I did some preliminary steps to avoid some obviously unnecessary computations:

  • I replaced all int(...) with evalf(Int(...)). Note the use of capital I!
  • I replaced all sum with add.

Next, I dealt with your calls to pprint. I assume that this is a procedure that you use to print supplementary information perhaps as an aid in debugging. The problem with that is that even if the procedure does nothing at all (indeed, even if the procedure doesn't even exist), it still takes a significant amount of time simply to evaluate its arguments. There is a built-in command made just for this situation: userinfo. I replaced all pprint(...with userinfo(1, Erweiterung, ...). This call will evaluate and print the arguments if and only if infolevel[Erweiterung] >= 1 where infolevel is a global table that's already predefined explicitly for use with userinfo. (It's traditional to use a procedure's name as the 2nd argument to userinfo, but actually any name can be used. The name must match the infolevel entry.)

Next, I used the command CodeTools:-Profiling:-Profile to "profile" Erweiterung. "Profiling" means to create a line-by-line analysis showing the number of times each line is executed and how much time and memory it used. This is cumulative over all calls to the procedure in the current session; so, if you want to analyze a single call to the procedure, start with restart. Thus:

#The procedure code and (very importantly!) the restart are above. 
#I just didn't want to clutter up this Answer with that code.

Digits:= 15: #experimental value

(*
#The following line of code is commented out because I *don't* 
#want the userinfo calls to take any time.
infolevel[Erweiterung]:= 1:
*)

CodeTools:-Profiling:-Profile(Erweiterung):

#This was your troublesome call:
Erweiterung(-1,1,2*x^2+2,2/sqrt(1-x^2),[-.8660254037, 0, .8660254037], 4);

DieKnoten, -0.999999999950535, -0.866025403700000, 
  -0.499999999951372, 0., 0.499999999951372, 0.866025403700000, 
  0.999999999950535

DieGewichte, TABLE([1 = 0.523598776221123, 2 = 1.04719755075381, 
  3 = 1.04719755105594, 4 = 1.04719755111737, 
  5 = 1.04719755105594, 6 = 1.04719755075381, 
  7 = 0.523598776221123])

                        18.8495559215366


To view the "profile", use CodeTools:-Profiling:-PrintProfiles:

CodeTools:-Profiling:-PrintProfiles(Erweiterung);
Erweiterung
Erweiterung := proc(Unten, Oben, f, g, Liste, n)
local Basenwechsel, p, i, c, d, e, Hn, Koeffizienten, a, s, j, M, V, S, K, nNeu, Em,
Hnm, KnotenHnm, KoeffizientenHnm, h0, b, gxi, Gewichte, Delta, Ergebnis, 
  Endergebnis, Koeffizient, Rest, Basenwechsel1, Basenwechsel2, k;
     |Calls Seconds  Words|
PROC |    1   1.438 17286788|
   1 |    1   0.000     37| Basenwechsel1 := proc( Dividend, m ) :: list;
                                local Basenwechsel;
                                Koeffizient := quo(Dividend,p[m],x);
                                Rest := rem(Dividend,p[m],x);
                                if m = 0 then
                                    Basenwechsel := [Koeffizient*evalf(Int(g*p
                                      [m]^2,x = Unten .. Oben))];
                                else
                                    Basenwechsel := [Koeffizient*evalf(Int(g*p
                                      [m]^2,x = Unten .. Oben)), op(procname(
                                      Rest,m-1))];
                                end if;
                            end proc;
   2 |    1   0.000     31| Basenwechsel2 := proc( Dividend, m ) :: list;
                                local Basenwechsel;
                                Koeffizient := quo(Dividend,p[m],x);
                                Rest := rem(Dividend,p[m],x);
                                if m = 0 then
                                    Basenwechsel := [Koeffizient];
                                else
                                    Basenwechsel := [Koeffizient, op(procname(
                                      Rest,m-1))];
                                end if;
                            end proc;
   3 |    1   0.000    273| p[-1] := 0;
   4 |    1   0.000     10| p[0] := 1;
   5 |    1   0.000      6| for i to 2*numelems(Liste)+2*n do
   6 |   14   1.266 15477600|     p[i] := x^i-add(evalf(Int(x^i*p[j]*g,x = Unten
                                  .. Oben))*p[j]/evalf(Int(p[j]^2*g,x = Unten
                                  .. Oben)),j = 0 .. i-1);
   7 |   14   0.000      0|     userinfo(1,Erweiterung,p[i]);
   8 |   14   0.000    778|     c[i-1] := coeff(p[i],x,i)/coeff(p[i-1],x,i-1);
   9 |   14   0.000   1149|     d[i-1] := (coeff(p[i],x,i-1)-c[i-1]*coeff(p[i-\
                                  1],x,i-2))/coeff(p[i-1],x,i-1);
  10 |   14   0.000     42|     if i <> 1 then
  11 |   13   0.000   2504|         e[i-1] := coeff(p[i]-(x*c[i-1]+d[i-1])*p[i
                                      -1],x,i-2)/coeff(p[i-2],x,i-2)
                                else
  12 |    1   0.000    275|         e[i-1] := 0
                                end if
                            end do;
  13 |    1   0.000      0| userinfo(1,Erweiterung,Liste[1],numelems(Liste));
  14 |    1   0.000     52| Hn := mul(x-Liste[i],i = 1 .. numelems(Liste));
  15 |    1   0.000      0| userinfo(1,Erweiterung,Hn);
  16 |    1   0.015  34905| Koeffizienten := PolynomialTools:-
                              FromCoefficientList(Basenwechsel1(Hn,numelems(
                              Liste)+1),x,termorder = reverse);
  17 |    1   0.000      0| userinfo(1,Erweiterung,Koeffizienten,HIER);
  18 |    1   0.000      0| userinfo(1,Erweiterung,c,d,e);
  19 |    1   0.000    549| a[0][0] := 1;
  20 |    1   0.000    286| a[1][0] := x;
  21 |    1   0.000     68| a[1][1] := -e[1]*c[0]/c[1]+(d[0]-d[1]*c[0]/c[1])*x
                              +c[0]/c[1]*x^2;
  22 |    1   0.000      2| for s from 2 to numelems(Liste)+n do
  23 |    6   0.000   1746|     a[s][0] := x^s;
  24 |    6   0.000    620|     a[s][1] := -e[s]*c[0]/c[s]*x^(s-1)+(d[0]-d[s]*
                                  c[0]/c[s])*x^s+c[0]/c[s]*x^(s+1);
  25 |    6   0.000    194|     pprint(coeff(a[s][1],x,s),as1s)
                            end do;
  26 |    1   0.000      2| for s from 2 to numelems(Liste)+n do
  27 |    6   0.000      0|     for j from 2 to s do
  28 |   21   0.000      0|         userinfo(1,Erweiterung,c[j-1]*add(coeff(a[
                                      s][j-1],x,k-1)/c[k-1]*x^k,k = abs(s-j)+2
                                      .. s+j));
  29 |   21   0.000      0|         userinfo(1,Erweiterung,add((d[j-1]-c[j-1]*
                                      d[k]/c[k])*coeff(a[s][j-1],x,k)*x^k,k =
                                      abs(s-j)+1 .. s+j-1));
  30 |   21   0.000      0|         userinfo(1,Erweiterung,c[j-1]*add(e[k+1]*
                                      coeff(a[s][j-1],x,k+1)/c[k+1]*x^k,k = 
                                      abs(s-j) .. s+j-2));
  31 |   21   0.000      0|         userinfo(1,Erweiterung,e[j-1]*add(coeff(a[
                                      s][j-2],x,k)*x^k,k = s-j+2 .. s+j-2));
  32 |   21   0.000  31718|         a[s][j] := c[j-1]*add(coeff(a[s][j-1],x,k-\
                                      1)/c[k-1]*x^k,k = abs(s-j)+2 .. s+j)+add
                                      ((d[j-1]-c[j-1]*d[k]/c[k])*coeff(a[s][j-\
                                      1],x,k)*x^k,k = abs(s-j)+1 .. s+j-1)-c[j
                                      -1]*add(e[k+1]*coeff(a[s][j-1],x,k+1)/c[
                                      k+1]*x^k,k = abs(s-j) .. s+j-2)+e[j-1]*
                                      add(coeff(a[s][j-2],x,k)*x^k,k = abs(s-j
                                      )+2 .. s+j-2)
                                end do
                            end do;
  33 |    1   0.000      2| for s from 0 to numelems(Liste)+n do
  34 |    8   0.000      0|     for j from 0 to s do
  35 |   36   0.000      0|         userinfo(1,Erweiterung,a[s][j],Polynom[s][
                                      j])
                                end do
                            end do;
  36 |    1   0.000   6724| M := Matrix(n,n);
  37 |    1   0.000   3919| V := Vector(n);
  38 |    1   0.000      0| for s from 0 to n-1 do
  39 |    4   0.000      0|     for j from 0 to s do
  40 |   10   0.000   1674|         M(s+1,j+1) := add(coeff(a[s][j],x,k)*coeff
                                      (Koeffizienten,x,k),k = 0 .. n);
  41 |   10   0.000     30|         if s <> j then
  42 |    6   0.000     54|             M(j+1,s+1) := M(s+1,j+1)
                                    end if;
  43 |   10   0.000      0|         userinfo(1,Erweiterung,M,1)
                                end do;
  44 |    4   0.000      0|     userinfo(1,Erweiterung,testb1);
  45 |    4   0.000      0|     userinfo(1,Erweiterung,coeff(a[n][s],x,2));
  46 |    4   0.000      0|     userinfo(1,Erweiterung,coeff(Koeffizienten,x,2
                                  ));
  47 |    4   0.000      0|     userinfo(1,Erweiterung,testb2);
  48 |    4   0.000      0|     userinfo(1,Erweiterung,Koeffizienten);
  49 |    4   0.000    955|     M(s+1,n+1) := add(coeff(a[n][s],x,k)*coeff(
                                  Koeffizienten,x,k),k = 0 .. n);
  50 |    4   0.000      0|     userinfo(1,Erweiterung,M,V)
                            end do;
  51 |    1   0.000      0| userinfo(1,Erweiterung,M,V);
  52 |    1   0.016 193237| S := LinearAlgebra:-LinearSolve(M,V);
  53 |    1   0.000     48| K := evalindets(S,name,() -> 2);
  54 |    1   0.000      0| userinfo(1,Erweiterung,K,LinSolve);
  55 |    1   0.000    260| Em := add(p[i]*K[i+1],i = 0 .. n);
  56 |    1   0.000     14| Hnm := Hn*Em;
  57 |    1   0.016 161926| KnotenHnm := fsolve(Hnm,complex);
  58 |    1   0.000      0| userinfo(1,Erweiterung,Hn,alt,Em,neu,Hnm);
  59 |    1   0.000      0| userinfo(1,Erweiterung,Testergebnis,nNeu);
  60 |    1   0.000      0| userinfo(1,Erweiterung,fsolve(Hnm),fsolve(nNeu));
  61 |    1   0.000  11230| KoeffizientenHnm := ListTools:-Reverse(
                              Basenwechsel2(Hnm,numelems(Liste)+n));
  62 |    1   0.000      0| userinfo(1,Erweiterung,KoeffizientenHnm);
  63 |    1   0.000     12| h0 := evalf(Int(g,x = Unten .. Oben));
  64 |    1   0.000      0| userinfo(1,Erweiterung,h0,HO);
  65 |    1   0.000    277| b[n+numelems(Liste)+2] := 0;
  66 |    1   0.000     14| b[n+numelems(Liste)+1] := 0;
  67 |    1   0.000      4| for i to nops([KnotenHnm]) do
  68 |    7   0.000     14|     for j from numelems(Liste)+n by -1 to 1 do
  69 |   49   0.000      0|         userinfo(1,Erweiterung,test21,KnotenHnm,
                                      Hnm);
  70 |   49   0.000   2108|         b[j] := KoeffizientenHnm[j+1]+(c[j]*
                                      KnotenHnm[i]+d[j])*b[j+1]+e[j+1]*b[j+2];
  71 |   49   0.000      0|         userinfo(1,Erweiterung,b[j])
                                end do;
  72 |    7   0.000      0|     userinfo(1,Erweiterung,test23);
  73 |    7   0.000  21442|     gxi := quo(Hnm,x-KnotenHnm[i],x);
  74 |    7   0.000      0|     userinfo(1,Erweiterung,gxi);
  75 |    7   0.000   1230|     Gewichte[i] := c[0]*b[1]*h0/eval(gxi,x = 
                                  KnotenHnm[i]);
  76 |    7   0.000      0|     userinfo(1,Erweiterung,b);
  77 |    7   0.000    361|     Delta[i] := c[0]*b[1]
                            end do;
  78 |    1   0.125 1293248| print(DieKnoten,KnotenHnm);
  79 |    1   0.000  34818| print(DieGewichte,Gewichte);
  80 |    1   0.000    340| Ergebnis := add(eval(f,x = KnotenHnm[k])*Gewichte[
                              k],k = 1 .. nops([KnotenHnm]));
  81 |    1   0.000      0| Endergebnis := evalf(Ergebnis)
end proc

And so we see that the line of code labelled 6 in the left column is where the vast majority of the the time is being used. Line 78 is also a problem, and likely should be replaced with userinfo. Interpretting the memory numbers is more of an expert-level topic and you shouldn't concern yourself with that at the moment. Regarding the times reported as 0.000: Maple can't measure execution times less than 0.015 seconds; they just get reported as 0. If you need them measured (which doesn't seem to be the case at the moment), then you need to accumulate several runs of the code in one profile.

Next, I decided to measure the effect of precision (as set via Digits) on performance. I set Digits to 30 and redid the above. (Extremely important: This is done starting with restart. Without that, the measurements are difficult or impossible to interpret. I simply don't show the whole thing here because other than Digits, the code is exactly the same.) Here's what I got:

restart:
#omitted procedure code ....
Digits:= 30:
#omitted code that sets up profiling ....
Erweiterung(-1,1,2*x^2+2,2/sqrt(1-x^2),[-.8660254037, 0, .8660254037], 4);

Error, (in Erweiterung) unable to compute coeff

So, you'll need to investigate that error further; and I'll look into it also. Nonetheless, the Profiling information collected up until the point of the error has still been collected, and in this case it's highly relevant: 

Erweiterung
Erweiterung := proc(Unten, Oben, f, g, Liste, n)
local Basenwechsel, p, i, c, d, e, Hn, Koeffizienten, a, s, j, M, V, S, K, nNeu, Em,
Hnm, KnotenHnm, KoeffizientenHnm, h0, b, gxi, Gewichte, Delta, Ergebnis, 
  Endergebnis, Koeffizient, Rest, Basenwechsel1, Basenwechsel2, k;
     |Calls Seconds  Words|
PROC |    1  70.578 932086456|
   1 |    1   0.000     37| Basenwechsel1 := proc( Dividend, m ) :: list;
                                local Basenwechsel;
                                Koeffizient := quo(Dividend,p[m],x);
                                Rest := rem(Dividend,p[m],x);
                                if m = 0 then
                                    Basenwechsel := [Koeffizient*evalf(Int(g*p
                                      [m]^2,x = Unten .. Oben))];
                                else
                                    Basenwechsel := [Koeffizient*evalf(Int(g*p
                                      [m]^2,x = Unten .. Oben)), op(procname(
                                      Rest,m-1))];
                                end if;
                            end proc;
   2 |    1   0.000     31| Basenwechsel2 := proc( Dividend, m ) :: list;
                                local Basenwechsel;
                                Koeffizient := quo(Dividend,p[m],x);
                                Rest := rem(Dividend,p[m],x);
                                if m = 0 then
                                    Basenwechsel := [Koeffizient];
                                else
                                    Basenwechsel := [Koeffizient, op(procname(
                                      Rest,m-1))];
                                end if;
                            end proc;
   3 |    1   0.000    273| p[-1] := 0;
   4 |    1   0.000     10| p[0] := 1;
   5 |    1   0.000      6| for i to 2*numelems(Liste)+2*n do
   6 |   12  70.578 932082225|     p[i] := x^i-add(evalf(Int(x^i*p[j]*g,x = Unten
                                  .. Oben))*p[j]/evalf(Int(p[j]^2*g,x = Unten
                                  .. Oben)),j = 0 .. i-1);
   7 |   12   0.000      0|     userinfo(1,Erweiterung,p[i]);
   8 |   12   0.000    687|     c[i-1] := coeff(p[i],x,i)/coeff(p[i-1],x,i-1);
   9 |   11   0.000    993|     d[i-1] := (coeff(p[i],x,i-1)-c[i-1]*coeff(p[i-\
                                  1],x,i-2))/coeff(p[i-1],x,i-1);
  10 |   11   0.000     33|     if i <> 1 then
  11 |   10   0.000   1886|         e[i-1] := coeff(p[i]-(x*c[i-1]+d[i-1])*p[i
                                      -1],x,i-2)/coeff(p[i-2],x,i-2)
                                else
  12 |    1   0.000    275|         e[i-1] := 0
                                end if
                            end do;

Every line after that has an execution count of 0 due to the error. We see that line 6 is still the time hog. Now look at lines 8 and 9. There is no conditional-statement branching between these lines, yet line has an execution count (11) one less than line ​​​​​​8 (12). How can that be? It's because the error must've occured on line 9.

Like this:

seq([seq](p), p= Iterator:-CartesianProduct([1,2,3] $ 2));

Let's say that you've already created the basic plot and assigned it to the name P:

P:= plot(...);

And let's say that your special point is (x,y), Then do

plots:-display(
    P, 
    plot(
        [[x,y]], 
        style= point, symbol= solidcircle, symbolsize= 20,
        color= red
    )
);

 

Here's a small modification to acer's final version that uses the very old but little-known type &under:

foo:= proc(
    ode::{`=`, And({set,list}(`=`), Not(0) &under nops)},
    func::function(name), 
    $
)  
    #...
end proc:

type(e, T &under F) is equivalent to type(F(e), T). I'm not posting this because I think this is any better than acer's; rather, I just like to promote the use of this type. "Not zero under nops" is very understandable to me.

By the way, the quotes that you had on `function` do absolutely nothing, neither good nor bad. So, I hope that you weren't expecting them to do something.

The following will get you most of the way to your final form. To finish, you need to apply the well-knowm "half-angle identity" for sqrt(1-cos(x)). Maple has a strong preference for cos instead of sin, and I can't persuade it to change its mind about that.

R:= 1/n*sum(exp(I*(i-1)*k*varepsilon), i= 1..n):
simplify(abs(R)) assuming real, n::posint;
                                           (1/2)
               (-2 cos(k n varepsilon) + 2)     
               ---------------------------------
                                          (1/2) 
               n (2 - 2 cos(k varepsilon))      

 

The command print never produces programmatically accessible output; the only way to access it is copy and paste.

I can't see any reason why you'd want your output as a set, although it's possible that you have a good reason.

Here are two syntax alternatives for the creation of a column Vector. The first two lines require relatively recent versions of Maple:

f:= x-> local y; fsolve(x*y=1):
R1:= <(for x from 1 by 0.2 to 2 do f(x) od)>;
R2:= <seq(f(x), x= 1..2, 0.2)>;

 

You need to correct the syntax in ode2. After that your dsolve with D(u)(0)=0 works perfectly:

ode2:= -v + mu*diff(r*diff(u(r), r), r)/r = 0:
dsolve({ode2, u(R)=0, D(u)(0)=0});
                                 2      2
                              v r    v R 
                       u(r) = ---- - ----
                              4 mu   4 mu

My corrections are shown in green. And, just to improve readability, I removed three levels of unnecessary parentheses.

It must be noted that while D(u)(0) = 0 implies u(0) is finite, the converse is far from true. D(u)(0) = 0 is a quite strong condition. 

The syntax for a multiple integral requires square brackets for the 2nd argument (the variables and ranges), like either of these:

simplify(Int(f(x,y,z)*Norm(n), [u= 1..2, v= 1..4]));
int(f(x,y,z)*Norm(n), [u= 1..2, v= 1..4], numeric);

However, your surface is undefined along the line y=3/4*x, which unfortunately passes through the rectangle (x= 1..2) x (y= 1..4). By some rearrangements of the variables, you can persuade int to give an answer, but that answer is Float(undefined).
 ​

This seems to me to be a bug in solve. The command eliminate, already shown by acer, is an algebraic generalization of solve, although it has far fewer options.

eqs:= {x*y^2-x, x*sin(Pi*y)}:
map2(eliminate, eqs, {x,y}); #Solve for each variable independently
         {[{x = 0}, {}], [{y = -1}, {}], [{y = 1}, {}]}

#Any return from eliminate whose 2nd element is not the empty set
#is not a complete solution (in the sense of solve). The possible
#presence of these incomplete solutions is what I mean by "algebraic
#generalization". There are no such solutions in this case; I include
#the next line only for the sake of completeness.
(C,InC):= selectremove(s-> op(2,s)={}, %);
      C, InC:= {[{x = 0}, {}], [{y = -1}, {}], [{y = 1}, {}]}, {}

#Remove the superfluous empty sets:
op~(..-2, C);
                  {{x = 0}, {y = -1}, {y = 1}}

#If you want to include the free variable in each solution:
map(solve, %, {x,y});
       {{x = 0, y = y}, {x = x, y = -1}, {x = x, y = 1}}


 

 

We must assume that discont returns all of the discontinuities, although it's pretty easy to have one's faith in this shaken. IMO computer algebra is not possible (at this time at least) without tacitly assuming that functions are piecewise continuous. 

The "flavors" used in RandomTools:-Generate can be nested arbitrarily. So, using the Vector flavor,

V:= RandomTools:-Generate('Vector'(variable(length = 1), 4));

The mul command will work for explicitly specified values of n. But to express the product for arbitrary n, you need to use product. Note how I iterate through the odd denominators only by using 2*a+1. Unlike (recent versions of) mul, the product command will not accept a "step" such as 2 as a 3rd argument. 

h:= x-> sin(x)/x:
J:= Int(product(h(x/(2*a+1)), a= 0..n), x= 0..infinity) 
    assuming n::nonnegint; 

Maple won't evaluate that integral; or at least it wouldn't do it in an amount of time that I was willing to wait. But it will do it for specified nonnegative integer values of n:

value(eval(J, n= 7));

For every n that I tried, it returned Pi/2. I suspect that there is a fairly easy proof that it's always Pi/2.

This seems too obvious---you must have considered it yourself: Is there any reason to not do simply

max(abs~(F-G))

where is the matrix for f and G is the matrix for g?

Aposthropes (ASCII code 39) are not special in strings. The only special ASCII characters are double quotes (ASCII 92) and backslash (ASCII 34). These can also be included in strings; they simply need to be escaped with backslash if they're hard-coded.

The error message that you show is typical when using sum, and you mentioned sum in your title. So which did you actually use? The correct command is add, not sum. Furthermore, the eval and indexing are not needed. You should be able to replace the entire command with

evalDG(add(L));

First 81 82 83 84 85 86 87 Last Page 83 of 395