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!

Hello everyone,

I think that I accidentely deleted a file from inside of a Maple workbook. Is there a way to restore it inside of that workbook?

Regards,

Hi,

I have managed to create the following plot.  It won't plot on the site's plotter:

fieldplot([1, y^2 + x], x = -10 .. 10, y = -6 .. 6, fieldstrength = fixed, color = abs(y^2 + x + 1))

I am trying to assign a colour gradient to the different vectors based on vector magnitude.  In the color option, I entered color= expresssion for the magnitude of the vectors

This only half worked.  It currently scales the colours such that the largest and smallest vectors are the same colour.  How do I assign a gradient such that the small magnitude vectors are one colour, and they then transition to another colour as their magnitude gets larger?

Hi,

What's the steps to prevent simplify() from doing ugly stuff like this :

Ofc I expect something like (1 - R) instead of -(-1+R)

Thanks!

Is it possible to generate a XMLElement by using part of the code to be generated with a string, or are the only methods a strict use of AddAttribute and AddChild?

Hope the code below explains what I mean. The first line is a strict copy of a declaration in the help file.

The question is if I could use some part of the definition which is saved in the str variable and combine it into the first declaration?

Download xml_string.mw

Is there a command that can show me all multipermutations?

If I have a set (1,2,3) I want an output of all the multipermutations of length n.

For example for n=2:
[1,1],[1,2],[1,3],[2,1],[2,2],[2,3],[3,1],[3,2],[3,3].

I suspect I have to use the Iteration command MultiPartition, but I don't understand the examples given on maplesoft:
https://www.maplesoft.com/support/help/Maple/view.aspx?path=Iterator%2FMultiPartition

Dear all,

In some step of my program, Maple cannot understand that the two following vectors are equal:

V1 := Vector[column](8, [1, 2, 2, 1, 3, A, B, 1/(A + B)^2]);

V2 := Vector[column](8, [1, 2, 2, 1, 3, A, B, 1/(A^2 + 2*A*B + B^2)]);

I tried to use the following two commands:

LinearAlgebra:-Equal(V1[6 .. 8], V2[6 .. 8]);

verify(V1, V2, 'Vector(expand)');

but Maple still returns 'false' instead of 'true'

Could somebody help me please ?
Best regards,

I'd like to define a proc, which takes first argument to be either an ode (i.e. type `=`) or set of ode's, or list of ode's.

However, I do not know how to tell Maple that the list or set, if that is the type, to be a list of `=` and no other type.

Here is what I tried to make it more clear

If I do this

restart;

interface(warnlevel=4);
kernelopts('assertlevel'=2): 


foo:=proc(ode::{`=`, set,list},func::`function`,$)
   print("ode=",ode);
   print("func=",func);
end proc:

Then Maple will check that the first argument is ANY one of `=`, set or list. 

But does not check if the list or set contains only equations of type `=`.  So I am able to call the above like this

ode1:=diff(y(x),x)=1:
ode2:=diff(y(x),x)=x:
foo(ode1,y(x));  #this is OK
foo([ode1,a],y(x))  #this is wrong

Which is wrong, since `a` is not of type `=`.

Next I tried this (I also wanted to check that if first argument is list or set, that it is not empty, so added extra check)

restart;
interface(warnlevel=4);
kernelopts('assertlevel'=2):

foo:=proc(ode::{`=`, 
                   And(set,satisfies(x-> (numelems(x)<>0 and type(x,`=`))  )),
                   And(list,satisfies(x-> (numelems(x)<>0 and type(x,`=`))  ))
                   },
                   func::`function`,$                                   
                   )

  print("ode=",ode);
  print("func=",func);
end proc;

But the above gives an error

ode1:=diff(y(x),x)=1:
ode2:=diff(y(x),x)=x:

foo([ode1,a],y(x))  

I also tried

foo:=proc(ode::{`=`, 
                   And(set,satisfies(x-> type(x,`=`)  )),
                   And(list,satisfies(x-> type(x,`=`) ))
                   },
                   func::`function`,$                                   
                   )

   print("ode=",ode);
   print("func=",func);

end proc;

Also gives erorr when called 

foo([ode1,ode2],y(x))

Ofcourse, I can just leave the check as in first case above, and in the proc itself, do the check myself manually by going over each entry in the list or set to make sure each entry is of type `=`. 

But I wanted Maple to do this work for me, if possible.

What is the correct syntax for doing so?

Maple 2020.2

[Perhaps this should be a Post rather than a Question.]

As far as I can tell, RootOf is the only Maple command that takes a bound variable on input and returns an unevaluated form with the bound variable renamed (to _Z). Why is that? Either all should do it (definite integrals, limit, sum, etc.) or none should. 

Pros:

Making a substitution for a variable that occurs both free and bound in the same expression is a major source of programming error. Renaming bound variables helps to ameliorate that. 

Cons:

The unexpected appearance of _Z in their results seems to be a great source of confusion to new users. However, if more commands did this, perhaps it would be more expected.

Feel free to start a brainstorming discussion on this, as if it were a Post, even if you don't have a direct Answer for the Question. If that's the way that the thread heads, I'll change it to a Post.

Hi,

I'm experiencing something quite annoying.

When copying Maple input code from an executable block into a document block, just after a line of text was written.

So for example <CTRL+T><write some text><CLICK ENTER FOR \newline><CTRL+M><PASTE MAPLE INPUT 1D CODE><CLICK ENTER>

The block doesn't get executed. I don't see the blue Maple output. Why is that? Is this a bug or I'm I doing something wrong. Normally logic would dictate after pressing enter it would execute the 1D code and simply show the blue output just below the 1D one-liner?

Another question:

I want to write the values of function g(x,t) as scientific notation.  I mean I want to write 1.95*10^(-3)  instead of 0.001953125 etc.

Yes, by clicking right-click and clicking the numeric formatting, I can transform all columns of Array A, but I want to transform just the last column g(x,t) in Array A.

restart:
  f:=(x,t)->x*t;
  g:=(x,t)->x^2*t;interface(rtablesize=20):
  A:=Array( [ [`x`,`t`,`f(x,t)`, `g(x,t)`],
              seq
              ( seq
                ( [i, j, f(i,j), g(i,j)],
                  j=0.125..0.875, 0.25
                ),
                i=0.125..0.875, 0.25
              )
            ]
          );  

lambda:=0.1:N:=5:M:=sqrt(N(N+1)):omegap:=10:phi:=0:

var:={n(t),u(t)}:
dsys:={diff(n(t),t)=-2*(n(t)-N)+(u(t)-abs(M)*exp(I*phi))*exp(-2*I*omegap*t/lambda)+conjugate((u(t)-abs(M)*exp(I*phi))*exp(-2*I*omegap*t/lambda)),diff(u(t),t)=-2*(1-I*delta)*u(t)+2*(n(t)-N)*exp(-2*I*omegap*t/lambda)+2*abs(M)*exp(I*phi)}:
res1:=dsolve(dsys union {n(0)=0,u(0)=0},numeric,output=listprocedure):

# (this command need to be changed) P1:=plots[odeplot3](res1,[[t,(n(t))]],0..10,axes=boxed,tickmarks = [3, 2], color = black, thickness = 4, linestyle = solid, titlefont = [Helvetica, roman, 18], labeldirections = [horizontal, vertical], labelfont = [Helvetica, roman, 24]);

The member command allows testing, if a value is in a table.

But is it possible to check if the index is used, not the value?

Download list.mw

I am trying to plot generator reactive output (Q) over a range of generator power output (P) with field current (Ifldc) and terminal voltage (Et) constant using “solve”. The system of equations includes an interpolating fiction Ifld(el).

I’m looking for some help as to how to configure a solution. Attached are a few of my failed attempts.

The upload of my file deleats the "solve" fo my first try it is:

Download Qcalc1.mw