Spirithaunter

75 Reputation

4 Badges

6 years, 247 days

MaplePrimes Activity


These are replies submitted by Spirithaunter

@acer Thank you!

Speed and precision is what is mostly needed, the matrix is symmetrical, other than that anything can appear.

 

evalindets(S,name,()->1);

has done it for me, though. Thank you very much! :)

@acer Perfect, that did it! Thank you very, very much!

 

@Spirithaunter  Basenwechsel.mw seems to be broken, this should work instead:

Basenwechsel_.mw

@acer Apologies, when I type version(), the result is this:

version()
 User Interface: 1194701
         Kernel: 1194701
        Library: 1194701
                            1194701

Checking for updates yields "No updates are available", therefore I thought I was up to date. I got Maple in the year 2016.

"Fixing" the unwanted local declaration meant removing it, there was not much else I thought I could do to get around that issue. "Correcting" is not the term I was looking for, I just fixed the problem of having a local declaration which Maple did not like. That probably leads to the old problem, yes, I still thought I could try it, the procname fix potentially could have done it alone. Thus, I did not include the code, which is just your proposal minus the declaration. I can, however, post it if you like.

Basenwechsel.mw

The issue is here:

Basenwechsel_Error.mw

The first use of the entire procedure is compiled with my Maple, the second is your computation of your original code.

You seemed confused of how I tried to program a recursive procedure, as in "(perhaps an attempt at a recursive call to itself)". Therefore I thought the entire structure was unnecessarily complicated. Apparently, your intended correction is the way of usually programming it, and not just your attempt at preserving my code (which I did not know), therefore I asked how it is done properly.

What do you mean by "The code I posted can be adjusted to declare locals at the start of the procedures, if that was an issue for your version." ? The declaration is right after the declaration of the subprocedure, so where should I put it?

 

Thank you again for your effort :)

Daniel

@acer Hi acer!

 

Thank you very much for your time in trying to answer my question, unfortunately, your solution does not work for my version of Maple, however. It complains about Basenwechsel having a local variable "Error, unexpected `local` declaration in procedure body" and even after fixing that the result is still [0,0,0,1,0].

Basenwechsel was indeed designed to be a recursive procedure. Given a set of orthogonal polynomials, it divides a polynomial by the orthogonal polynomial of the same degree d and then divides the remainder (of degree up to d-1) by the orthogonal polynomial of the degree d-1, and so on, down to 0. The result of these divisions are the coefficients needed to express the polynomial entered into the recursive procedure in terms of the orthogonal polynomials.

Declaring the locals is something I wanted to do after my code is finished, I think this is easier :)

 

Is there a better way to do recursive algorithms then?

 

Greetings, Daniel

@Carl Love Thank you :D

@Carl Love Thank you sir, that problem in the code is fixed now :)

@Carl Love That seems to be even easier, thank you :D

@Kitonum Thank you very much :) I think the seq with i=1..n was too much, but without it it works perfectly :D

@acer 15035

Thank you! :D

The problem remained after the name changes, however evaluating the eigenvectors helped. The computer still mentions

"Error, (in Matrix) triangular[] storage expects name parameters" , however it does calculate the eigenvectors and -values.  :)

 

@Kitonum 

Thank you for your answers so far :)

I have done that now, it now works for the matrix in the example, not for mine however.

Does the matrix have to be generated in a specific way for it to work?

Here is my code. The matrix is generated at the bottom, its values are correct, I checked that. 
with(LinearAlgebra): is in the 6th line.

nodes, vac := Eigenvectors(M)  still just doesn't do anything.

And what does "Error, (in Matrix) triangular[] storage expects name parameters" even mean?

 

lower:= -1;
upper:= 1;
G:= x;
f:=x;
n:=2;
   
 
with(LinearAlgebra):
 
 
A := [seq(a[i], i = 0 .. n)];
B := [seq(b[i], i = 0 .. n)];
P := [seq(p[i], i = -1 .. 2*n+1)];  #?
S := [seq(s[i], i = -1 .. floor(n/2))];
T := [seq(t[i], i = -1 .. floor(n/2))];
p[-1]:= 0;
p[0]:=1;
for v from -1 to floor(n/2) do
  s[v]:=0;
  t[v]:=0
end do;
for j from 1 to 2*n+1 do
  RekursivesZwischenergebnis:= x^j;
  for i from 0 to j-1 do
    RekursivesZwischenergebnis:= RekursivesZwischenergebnis -
    (int(x^j*p[i]*diff(G,x),x=lower..upper)/int(p[i]*p[i]*diff(G,x),x=lower..upper))*p[i]                  #Gram-Schmidt algorithm
  end do;
  p[j]:=RekursivesZwischenergebnis;
end do;
a[0]:=-coeff(p[1],x,0);

  #p[0+1]=(x-a[0])*p[0]-b[0]*p[0-1] -> p[1]=x*p[0]-a[0]*p[0]-b[0]*p[-1] ->
  #p[1]=x*1-a[0]*1-0 -> a[0]=x-p[1] -> a[0]= -coeff(p[1],x,0) since p[1] is monic and of degree 1
 
b[0]:=int(p[0]^2, x=lower..upper); #by definition
for j from 1 to ceil(3*n/2) do
 
  #to be exact, a only needs to be initialised up to floor(3/(2*n)). It does not hurt, however, to   #do it further. The value for a[ceil(3/(2*n))] is needed for a[ceil(3/(2*n)]anyway. Not cutting
  #the initialisation for a[ceil(3/(2*n))] if ceil(3/(2*n) is unequal to floor(3/(2*n)) saves        #code and the evaluation of the comparison.
    
                                     
  a[j]:= coeff(p[j],x,j-1)- coeff(p[j+1],x,j);
    
    #p[j+1]=(x-a[j])*p[j]-b[j]*p[j-1] -> p[j+1]=x*p[j]-a[j]*p[j]-b[j]p[j-1] ->
    #coeff(p[j+1],x,j)=coeff(x*p[j],x,j)-coeff(a[j]*p[j],x,j)
      #(since b[j]*p[j-1] is of degree j-1) ->
    #coeff(p[j+1],x,j)=coeff(x*p[j],x,j)-a[j], since p[j] is monic ->
    #coeff(p[j+1],x,j)=coeff(p[j],x,j-1)-a[j]->
    #a[j]=coeff(p[j],x,j-1)-coeff(p[j+1],x,j)
 
  b[j]:=  quo((x-a[j])*p[j]-p[j+1],p[j-1],x);
    

     #p[j+1]=(x-a[j])*p[j]-b[j]*p[j-1] -> -p[j+1]+(x-a[j])*p[j]= b[j]*p[j-1]
     #b[j]=((x-a[j])*p[j]-p[j+1])/p[j-1]

  end do;    
t[0]:=b[n+1];
for m from 0 to n-2 do
  u:=0;
  for k from floor((m+1)/2) to 0 by -1 do
    l:=m-k;
    u:=u+(a[k+n+1]-a[l])*t[k]+b[k+n+1]*s[k-1]-b[l]*s[k];
    s[k]:=u
  end do;
  for v from -1 to floor(n/2) do
    Hilfsvariable:=s[v];
    s[v]:=t[v];
    t[v]:=Hilfsvariable
  end do;
end do;
for j from floor(n/2) to 0 by -1 do
    s[j]:=s[j-1]
end do;
for m from n-1 to 2*n-3 do
  u:=0;
  for k from m+1-n to floor((m-1)/2) do
    l:=m-k;
    j:=n-1-l;
    u:=u-(a[k+n+1]-a[l])*t[j]-b[k+n+1]*s[j]+b[l]*s[j+1];
    s[j]:=u
  end do;
  if m mod 2 = 0 then
    k:= m/2;
    a[k+n+1]:=a[k]+(s[j]-b[k+n+1]*s[j+1])/t[j+1]
  else
    k:=(m+1)/2;
    b[k+n+1]:=s[j]/s[j+1]
  end if;
  for v from -1 to floor(n/2) do
    Hilfsvariable:=s[v];
    s[v]:=t[v];
    t[v]:=Hilfsvariable
  end do;
end do;
a[2*n]:=a[n-1]-b[2*n]*s[0]/t[0];
print(A);
print(B);
M:=Matrix(2*n+1);
M(1,1):=a[0];
for m from 2 to (2*n+1) do #generate Gauss-Kronrod-Matrix
  M(m-1,m):= sqrt(b[m-1]);
  M(m,m-1):= sqrt(b[m-1]);
  M(m,m):= a[m-1];
end do;
nodes, vac := Eigenvectors(M);
weights:=Vector(2*n+1);

print(nodes);
print(vac);
for m from 1 to 2*n+1 do
  CurrentNormalizedVector:= Normalize(Column(vac,m),Euclidean);
  print(CurrentNormalizedVector);
  weights[m]=CurrentNormalizedVector[1]^2*B[m];
  print(weights[m])
end do;
RekursivesZwischenergebnis:=0;
for i from 1 to 2*n+1 do
  RekursivesZwischenergebnis:=RekursivesZwischenergebnis + weights[i]*eval(f,x=nodes[i])
end do;
FinalResult:=RekursivesZwischenergebnis

 

I forgot to write it down here, apologies for that. In my actual program I have written "with(LinearAlgebra)"

 

1 2 Page 2 of 2