Question: On seeing what Maple is taking so long for

Hi!

I have the following code to calculate the optimal quadrature rule with preassigned nodes (within a list), while the number of nodes that are to be added is n. The calculated quadrature rule is then used to approximate an integral.

restart;
with(LinearAlgebra):     
with(ListTools):
with(PolynomialTools):
Erweiterung:= proc(Unten, Oben, f,g,Liste,n)::real;
  #Unten:= Untere Intervallgrenze; Oben:= Obere Intervallgrenze; f:= zu integrierende Funktion;
  #G:= Gewicht; Liste:= Liste der alten Knoten, n:= Anzahl hinzuzufügender Knoten;
  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;
  uses LinearAlgebra, ListTools;
  Basenwechsel1:=proc(Dividend, m)::list; #stellt ein gegebenes Polynom über eine Linearkombination der orthogonalen Polynomen dar.
   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;
  Basenwechsel2:=proc(Dividend, m)::list; #stellt ein gegebenes Polynom über eine Linearkombination der orthogonalen Polynomen dar.
   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;
p[-1]:=0;
p[0]:=1;
for i from 1 to (numelems(Liste)+n)*2 do
  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)); #Berechnung einer Folge orthogonaler Polynome bezüglich der gegebenen Gewichtsfunktion und des gegebenes Intervalles
  pprint(p[i]);
c[i-1]:=coeff(p[i],x,i)/coeff(p[i-1],x,i-1); #Berechnung der dreigliedrigen Rekursion der errechneten orthogonalen Polynome
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));
if i <> 1 then
  e[i-1]:=coeff(p[i]-(c[i-1]*x+d[i-1])*p[i-1],x,i-2)/coeff(p[i-2],x,i-2);
else
  e[i-1]:=0;
end if;
end do;
pprint(Liste[1],numelems(Liste));
Hn:=mul(x-Liste[i],i=1..numelems(Liste));
pprint(Hn);
 Koeffizienten:=FromCoefficientList(Basenwechsel1(Hn,numelems(Liste)+1),x,termorder=reverse); #Die Koeffizienten der orthogonalen Polynome werden hier als Koeffizienten der Monome gespeichert.
pprint(Koeffizienten,HIER);

pprint(c,d,e);
a[0][0]:=1; #Beginn der Berechnung der orthogonalen Produkterweiterungen, die Koeffizienten der orthogonalen Polynome werden wieder über die Monome gespeichert (2*x^2+2 bedeutet bspw. [2,0,2,0,0...] für die Koeffizienten)
a[1][0]:=x;
a[1][1]:=-e[1]*c[0]/c[1]+(d[0]-d[1]*c[0]/c[1])*x+c[0]/c[1]*x^2;
for s from 2 to numelems(Liste)+n do
  a[s][0]:=x^s;
  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);
    pprint (coeff(a[s][1],x,s),as1s);
end do;
for s from 2 to numelems(Liste)+n do
  for j from 2 to s do
    
      pprint(c[j-1]*sum(coeff(a[s][j-1],x,k-1)/c[k-1]*x^k,k=abs(s-j)+2..s+j));  pprint(sum((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));  pprint(c[j-1]*sum(e[k+1]*coeff(a[s][j-1],x,k+1)/c[k+1]*x^k,k=abs(s-j)..s+j-2));pprint(e[j-1]*sum(coeff(a[s][j-2],x,k)*x^k,k=s-j+2..s+j-2));

     a[s][j]:=c[j-1]*sum(coeff(a[s][j-1],x,k-1)/c[k-1]*x^k,k=abs(s-j)+2..s+j)+sum((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]*sum(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]*sum(coeff(a[s][j-2],x,k)*x^k,k=abs(s-j)+2..s+j-2);

      
   
    
  end do;
end do;
for s from 0 to numelems(Liste)+n do
  for j from 0 to s do
    pprint(a[s][j], Polynom[s][j]);
  end do;
end do;
M:=Matrix(n,n); #Beginn der Erstellung eines linearen Gleichungssystems, dessen Lösung die Koeffizienten der orthogonalen Polynome sind, deren Summe Em die hinzuzufügenden Knoten als Nullstellen hat.
V:=Vector(n);
 
  for s from 0 to n-1 do
    for j from 0 to s do
      M(s+1,j+1):=sum(coeff(a[s][j],x,k)*coeff(Koeffizienten,x,k),k=0..n);
      if s<>j then
        M(j+1,s+1):=M(s+1,j+1);
      end if;
      pprint(M,1);
    end do;
    pprint(testb1);pprint(coeff(a[n][s],x,2));pprint(coeff(Koeffizienten,x,2));
    pprint(testb2); pprint(Koeffizienten);
    M(s+1,n+1):=sum(coeff(a[n][s],x,k)*coeff(Koeffizienten,x,k),k=0..n);
    
    pprint(M,V);
  end do;
pprint(M,V);
S:=LinearSolve(M,V);
K:=evalindets(S,name,()->2);
pprint(K,LinSolve);

Em:=add(p[i]*K[i+1],i=0..n); #Erstellen von Em, dessen Nullstellen die hinzuzufügenden Knoten sind
Hnm:=Hn*Em; #Erstellen von Hnm, welches alle Knoten als Nullstelle besitzt
KnotenHnm:=fsolve(Hnm,complex); #Knotenberechnung

   
pprint(Hn,alt,Em,neu,Hnm);
pprint(Testergebnis,nNeu);
pprint(fsolve(Hnm),fsolve(nNeu));
KoeffizientenHnm:=Reverse(Basenwechsel2(Hnm,n+numelems(Liste)));  #Das Polynom Hnm wird über die orthogonalen Polynome dargestellt.
pprint(KoeffizientenHnm);
h0:=int(g,x=Unten..Oben); #Beginn der Berechnung der Gewichte
 pprint(h0,HO);
b[n+numelems(Liste)+2]:=0;
b[n+numelems(Liste)+1]:=0;
  for i from 1 to nops([KnotenHnm]) do
    for j from n+numelems(Liste) by -1 to 1 do
      pprint(test21,KnotenHnm,Hnm);
      b[j]:=KoeffizientenHnm[j+1]+(d[j]+KnotenHnm[i]*c[j])*b[j+1]+e[j+1]*b[j+2];
      pprint(b[j]);
    end do;
    pprint(test23);
    gxi:=quo(Hnm,x-KnotenHnm[i],x);
   pprint(gxi);
    Gewichte[i]:=c[0]*b[1]*h0/eval(gxi,x=KnotenHnm[i]);
    pprint(b);
   
    Delta[i]:=c[0]*b[1];
  end do;
print(DieKnoten,KnotenHnm);
print(DieGewichte, Gewichte);
Ergebnis:=add(eval(f,x=KnotenHnm[k])*Gewichte[k],k=1..nops([KnotenHnm]));
Endergebnis:=evalf(Ergebnis)
end proc

The problem is that the code takes very, very long to run if the weight function is not a polynomial.

Erweiterung(-1,1,2*x^2+2,1,[-0.906179845938664],4)

for example, is done immediately (1, the 4th entry, being the weight function), while

Erweiterung(-1,1,2*x^2+2,2/sqrt(1-x^2),[-.8660254037, 0, .8660254037],4)

takes ages to finish. Is there a tool for me to see what exactly is taking Maple so long? Is there an easy fix, such as evalf()'ing key calculations (other than using (2*x^2+2)*2/sqrt(1-x^2) as the integrand and 1 as the weight function, since the quadrature rules I am looking at are supposed to be good with certain weird weight functions)? Thank you in advance!