Good evening,

is it possible to construct a math table with Maple just with Code? Something like

      test_1    test_2

a       9,6        5,6

b     pending  8,4     

to be able to print it? There is

https://de.maplesoft.com/support/help/Maple/view.aspx?path=worksheet/documenting/table            and

https://de.maplesoft.com/support/help/Maple/view.aspx?path=worksheet%2freference%2ftableproperties   

, however, I do not understand how to just get any table going and assign calculated values within a procedure as is possible with a matrix for example. Thank you in advance! :)

Hi there!

I have a procedure that compares the (2n+1)-point Gauß-Kronrod-Quadrature to the (2n+1)-point Patterson-Quadratures for a range of n. I have plotted the results (the absolute and relative error if they "exist", meaning they need to posess certain features) in a graph, however they do not look very insightful. For the lowest n, the reader gets an impression for the different accuracy of the quadrature rules, however, for higher n's, the resulting points are basically just on the x-axis with no difference to see. Is it possible to also print a math table with Maple? Something like:

      GKQ   PZ+   PZ-   PY 

1     1,04   1,03  1,02  1,02

2     1,09   1,04     -       -

3     1,02   1,01  1,01  1,01

4     1,03   1,02  1,01  1,01

with - meaning no existance for that particular n? I havent found anything about that on the internet, it's all about plotting.

My (long) code is this:

restart:
with(LinearAlgebra):     
with(ListTools):
with(PolynomialTools):
with(CurveFitting):
with(plots):
Plotting:=proc(Unten,Oben,f,g,nUnten,nOben)::plot;

local SpeicherlisteX, SpeichervektorX, #speichert die Stützstellen
SpeichervektorXGekürzt, #streicht nicht existierende Quaraturformeln.
SpeicherlisteYAbs, SpeichervektorYAbs,  #speichert die Stützwerte des späteren Splines aus dem absoluten Quadraturfehler
SpeicherlisteYRel, SpeichervektorYRel,  #speichert die Stützwerte des späteren Splines aus den relativen Quadraturfehler
î, #Laufvariable  
InterpolationsfunktionAbs, #speichert den Spline aus dem absoluten Interpolationsfehler                      
InterpolationsfunktionRel, #speichert den Spline aus den relativen Fehlern von f
GraphAbsGK, GraphAbsPY, GraphAbsPZP, GraphAbsPZM, #speichert den Graphen aus dem Spline aus dem absoluten Interpolationsfehler          GraphRelGK, GraphRelPY, GraphRelPZP, GraphRelPZM, #speichert den Graphen aus dem Spline aus den relativen Fehlern  von f  
PunkteAbsGK, PunkteAbsPY, PunkteAbsPZP, PunkteAbsPZM,#speichert den Punktgraphen aus dem absoluten Interpolationsfehler
PunkteRelGK, PunkteRelPY, PunkteRelPZP, PunkteRelPZM, #speichert den Punktgraphen aus dem absoluten Interpolationsfehler
NichtexistenzGK, NichtexistenzPY, NichtexistenzPZP, NichtexistenzPZM, #speichert die Häufigkeit der Nichtexistenz

p,i,c,d,e,Hn,Koeffizienten,s,j,M,V,S,K,nNeu,Em,Hnm,KnotenHnm,KoeffizientenHnm,h0,b,gxi,Gewichte,Delta,Ergebnis,
Endergebnis,Koeffizient,Rest,a,VorgegebeneKnoten,TatsächlicherWert, DoppelterKnoten, KomplexerKnoten,

Text:= proc() #Prozedur zum Schreiben der Ausgabe
uses T= Typesetting;
     T:-mrow(seq(`if`(e::string, T:-mn(e), T:-Typeset(T:-EV(e))), e= [args]))
end proc,
OrtPol:= proc(G,N)::list; #Prozedur zum Berechnen der benötigten orthogonalen Polynome
  local q,r,R;
  q[-1]:=0;
  q[0]:=1;
 
  for r from 1 to N do
  q[r]:=(x^r-add(evalf(Int(x^r*q[R]*G,x=(-1)..1))*q[R]/evalf(Int(q[R]^2*G,x=(-1)..1)),R=0..r-1));
  end do;
  return(fsolve(q[N]));
end proc,
BasenwechselNormiert:=proc(Dividend, m)::list; #stellt ein gegebenes Polynom über eine Linearkombination der orthogonalen Polynome #dar.
   local BasenwechselNormiert;
 
  Koeffizient:=quo(Dividend, p[m],x);

  Rest:=rem(Dividend, p[m],x);
 
  if m=0 then
    BasenwechselNormiert:=[Koeffizient*evalf(Int(g*p[m]^2,x=Unten..Oben))];
  else

    BasenwechselNormiert:=[Koeffizient*evalf(Int(g*p[m]^2,x=Unten..Oben)),op(procname(Rest,m-1))];
   
  end if;
 
  end proc,
    Basenwechsel:=proc(Dividend, m)::list; #stellt ein gegebenes Polynom über eine Linearkombination der orthogonalen Polynome 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,
Erweiterung:= proc(Unten, Oben, f,g,Liste,n)::real; #Prozedur zur Berechnung der optimalen Erweiterung nach Knotenvorgabe
  #Unten:= Untere Intervallgrenze; Oben:= Obere Intervallgrenze; f:= zu integrierende Funktion;
  #g:= Gewicht; Liste:= Liste der alten Knoten, n:= Anzahl hinzuzufügender Knoten;
 
 
 
Hn:=mul(x-Liste[i],i=1..numelems(Liste));

 Koeffizienten:=FromCoefficientList(BasenwechselNormiert(Hn,numelems(Liste)+1),x,termorder=reverse); #Die Koeffizienten der orthogonalen Polynome werden hier als Koeffizienten der Monome gespeichert.

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):=add(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;
    end do;
    
    M(s+1,n+1):=add(coeff(a[n][s],x,k)*coeff(Koeffizienten,x,k),k=0..n);
    
    
  end do;

S:=LinearSolve(M,V);
K:=evalindets(S,name,()->2);

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

if (KnotenHnm[1]<-1-10^(-10)) or (KnotenHnm[n+numelems(Liste)]>1+10^(-10)) then
  return(false)
else
KomplexerKnoten:=false;
for i from 1 to n+numelems(Liste) do

 if(Im(KnotenHnm[i])>10^(-10)) then
  KomplexerKnoten:=true
 end if;
end do;
if KomplexerKnoten=true then
  return(false)
else
DoppelterKnoten:=false;
for i from 1 to n+numelems(Liste)-1 do
 
 if (KnotenHnm[i+1]-KnotenHnm[i]<10^(-10)) then
   DoppelterKnoten:=true
 end if;
end do;
if DoppelterKnoten=true then
 return(false)
else

KoeffizientenHnm:=Reverse(Basenwechsel(Hnm,n+numelems(Liste)));  #Das Polynom Hnm wird über die orthogonalen Polynome dargestellt.

h0:=evalf(Int(g,x=Unten..Oben)); #Beginn der Berechnung der Gewichte
 
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
      
      b[j]:=KoeffizientenHnm[j+1]+(d[j]+KnotenHnm[i]*c[j])*b[j+1]+e[j+1]*b[j+2];
      
    end do;
    
    gxi:=quo(Hnm,x-KnotenHnm[i],x);
   
    Gewichte[i]:=c[0]*b[1]*h0/eval(gxi,x=KnotenHnm[i]);
    
   
    Delta[i]:=c[0]*b[1];
  end do;

Ergebnis:=add(eval(f,x=KnotenHnm[k])*Gewichte[k],k=1..nops([KnotenHnm]));

Endergebnis:=Re(evalf(Ergebnis))
end if;
end if;
end if;
end proc:

p[-1]:=0;
p[0]:=1;
for i from 1 to (2*nOben+1)*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
 
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;

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 2*nOben+1 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 2*nOben+1 do
  for j from 2 to s do

     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;
for î from nUnten to nOben do
  VorgegebeneKnoten[î]:=OrtPol(g,î);
end do;
TatsächlicherWert:=evalf(Int(f*g,x= Unten..Oben));
GraphAbsGK:=plot([]); PunkteAbsGK:=plot([]); GraphAbsPZP:=plot([]); PunkteAbsPZP:=plot([]); GraphAbsPZM:=plot([]); PunkteAbsPZM:=plot([]); GraphAbsPY:=plot([]); PunkteAbsPY:=plot([]);
GraphRelGK:=plot([]); PunkteRelGK:=plot([]); GraphRelPZP:=plot([]); PunkteRelPZP:=plot([]); GraphRelPZM:=plot([]); PunkteRelPZM:=plot([]); GraphRelPY:=plot([]); PunkteRelPY:=plot([]);
SpeicherlisteX:=[];
SpeicherlisteYAbs:=[];
SpeicherlisteYRel:=[];
for î from nUnten to nOben do
  if Erweiterung(Unten,Oben,f,g,[VorgegebeneKnoten[î]],î+1) <> false then
    SpeicherlisteX:=[op(SpeicherlisteX),î]; #Stützstellen definieren                                   
    SpeicherlisteYAbs:=[op(SpeicherlisteYAbs),Erweiterung(Unten,Oben,f,g,[VorgegebeneKnoten[î]],î+1)-evalf(Int(f*g,   x=Unten..Oben))]; #Bestimmen des absoluten Fehlers von f für n=î
      if abs(TatsächlicherWert) > 10^(-10) then #Bestimmen des relativen Fehlers von f1 falls                                                          #dieser definiert ist
      SpeicherlisteYRel:=[op(SpeicherlisteYRel),abs(SpeicherlisteYAbs[-1]/TatsächlicherWert)];
    end if;
  end if;
end do;
if numelems(SpeicherlisteX)>0 then
  SpeichervektorX:=Vector[row](numelems(SpeicherlisteX),SpeicherlisteX);
  SpeichervektorYAbs:=Vector[row](numelems(SpeicherlisteYAbs),SpeicherlisteYAbs);
  PunkteAbsGK:= plot(SpeichervektorX,SpeichervektorYAbs,style = point, color=red, legend = ["GK"]);
    #  Generierung des Punktgraphen, der sich aus den absoluten Fehlern von f ergibt
  if numelems(SpeicherlisteX)>1 then
    InterpolationsfunktionAbs:=Spline(SpeichervektorX,SpeichervektorYAbs,n);
      #  Splines aus Stützpunkten, die sich aus den absoluten Fehlern von f ergeben
    GraphAbsGK:= plot(InterpolationsfunktionAbs, n=nUnten..nOben, color=red);
      #  Generierung des Graphen, der sich aus dem Spline aus den absoluten Fehlern von f ergibt
  end if;
end if;

if abs(TatsächlicherWert) > 10^(-10) then
  # falls der relative Fehler definiert ist analoges Vorgehen für die relativen Fehler
  if numelems(SpeicherlisteX)>0 then
    SpeichervektorYRel:=Vector[row](numelems(SpeicherlisteYRel),SpeicherlisteYRel);
    PunkteRelGK:= plot(SpeichervektorX,SpeichervektorYRel,style = point, color=red, legend = ["GK"]);
   
    if numelems (SpeicherlisteX)>1 then
      InterpolationsfunktionRel:=Spline(SpeichervektorX,SpeichervektorYRel,n);
      GraphRelGK:= plot(InterpolationsfunktionRel, n=nUnten..nOben, color=red);
    end if;
  end if;
end if;
NichtexistenzGK:=nOben-nUnten+1-numelems(SpeicherlisteX);

SpeicherlisteX:=[]; # analoges Vorgehen für PZP
SpeicherlisteYAbs:=[];
SpeicherlisteYRel:=[];
for î from nUnten to nOben do
  if Erweiterung(Unten,Oben,f,g,[-1,VorgegebeneKnoten[î]],î) <> false then
    SpeicherlisteX:=[op(SpeicherlisteX),î]; #Stützstellen definieren                                   
    SpeicherlisteYAbs:=[op(SpeicherlisteYAbs),Erweiterung(Unten,Oben,f,g,[-1,VorgegebeneKnoten[î]],î)-TatsächlicherWert]; #Bestimmen des absoluten Fehlers von f für n=î
      if abs(TatsächlicherWert) > 10^(-10) then #Bestimmen des relativen Fehlers von f1 falls                                                          #dieser definiert ist
      SpeicherlisteYRel:=[op(SpeicherlisteYRel),abs(SpeicherlisteYAbs[-1]/TatsächlicherWert)];
    end if;
  end if;
end do;
if numelems(SpeicherlisteX)>0 then
  SpeichervektorX:=Vector[row](numelems(SpeicherlisteX),SpeicherlisteX);
  SpeichervektorYAbs:=Vector[row](numelems(SpeicherlisteYAbs),SpeicherlisteYAbs);
  PunkteAbsPZP:= plot(SpeichervektorX,SpeichervektorYAbs,style = point, color=orange);
    #  Generierung des Punktgraphen, der sich aus den absoluten Fehlern von f ergibt
  if numelems(SpeicherlisteX)>1 then
    InterpolationsfunktionAbs:=Spline(SpeichervektorX,SpeichervektorYAbs,n);
      #  Splines aus Stützpunkten, die sich aus den absoluten Fehlern von f ergeben
    GraphAbsPZP:= plot(InterpolationsfunktionAbs, n=nUnten..nOben, color=orange);
      #  Generierung des Graphen, der sich aus dem Spline aus den absoluten Fehlern von f ergibt
  end if;
end if;

if abs(TatsächlicherWert) > 10^(-10) then
  # falls der relative Fehler definiert ist analoges Vorgehen für die relativen Fehler
  if numelems(SpeicherlisteX)>0 then
    SpeichervektorYRel:=Vector[row](numelems(SpeicherlisteYRel),SpeicherlisteYRel);
    PunkteRelPZP:= plot(SpeichervektorX,SpeichervektorYRel,style = point, color=orange, legend = ["PZP"]);
    
    if numelems (SpeicherlisteX)>1 then
      InterpolationsfunktionRel:=Spline(SpeichervektorX,SpeichervektorYRel,n);
      GraphRelPZP:= plot(InterpolationsfunktionRel, n=nUnten..nOben, color=orange);
    end if;
  end if;
end if;
NichtexistenzPZP:=nOben-nUnten+1-numelems(SpeicherlisteX);
SpeicherlisteX:=[];# analoges Vorgehen für PZM
SpeicherlisteYAbs:=[];
SpeicherlisteYRel:=[];
for î from nUnten to nOben do
  if Erweiterung(Unten,Oben,f,g,[VorgegebeneKnoten[î],1],î) <> false then
    SpeicherlisteX:=[op(SpeicherlisteX),î]; #Stützstellen definieren                                   
    SpeicherlisteYAbs:=[op(SpeicherlisteYAbs),Erweiterung(Unten,Oben,f,g,[VorgegebeneKnoten[î],1],î)-TatsächlicherWert]; #Bestimmen des absoluten Fehlers von f für n=î
      if abs(TatsächlicherWert) > 10^(-10) then #Bestimmen des relativen Fehlers von f1 falls                                                          #dieser definiert ist
      SpeicherlisteYRel:=[op(SpeicherlisteYRel),abs(SpeicherlisteYAbs[-1]/TatsächlicherWert)];
    end if;
  end if;
end do;
if numelems(SpeicherlisteX)>0 then
  SpeichervektorX:=Vector[row](numelems(SpeicherlisteX),SpeicherlisteX);
  SpeichervektorYAbs:=Vector[row](numelems(SpeicherlisteYAbs),SpeicherlisteYAbs);
  PunkteAbsPZM:= plot(SpeichervektorX,SpeichervektorYAbs,style = point, color=blue, legend = ["PZM"]);
    #  Generierung des Punktgraphen, der sich aus den absoluten Fehlern von f ergibt
  if numelems(SpeicherlisteX)>1 then
    InterpolationsfunktionAbs:=Spline(SpeichervektorX,SpeichervektorYAbs,n);
      #  Splines aus Stützpunkten, die sich aus den absoluten Fehlern von f ergeben
    GraphAbsPZM:= plot(InterpolationsfunktionAbs, n=nUnten..nOben, color=blue);
      #  Generierung des Graphen, der sich aus dem Spline aus den absoluten Fehlern von f ergibt
  end if;
end if;

if abs(TatsächlicherWert) > 10^(-10) then
  # falls der relative Fehler definiert ist analoges Vorgehen für die relativen Fehler
  if numelems(SpeicherlisteX)>0 then
    SpeichervektorYRel:=Vector[row](numelems(SpeicherlisteYRel),SpeicherlisteYRel);
    PunkteRelPZM:= plot(SpeichervektorX,SpeichervektorYRel,style = point, color=blue, legend = ["PZM"]);
 
    if numelems (SpeicherlisteX)>1 then
      InterpolationsfunktionRel:=Spline(SpeichervektorX,SpeichervektorYRel,n);
      GraphRelPZM:= plot(InterpolationsfunktionRel, n=nUnten..nOben, color=blue);
    end if;
  end if;
end if;
NichtexistenzPZM:=nOben-nUnten+1-numelems(SpeicherlisteX);
SpeicherlisteX:=[]; #analoges Vorgehen für PY
SpeicherlisteYAbs:=[];
SpeicherlisteYRel:=[];
for î from nUnten to nOben do
  if Erweiterung(Unten,Oben,f,g,[-1,VorgegebeneKnoten[î],1],î-1) <> false then
    SpeicherlisteX:=[op(SpeicherlisteX),î]; #Stützstellen definieren                                   
    SpeicherlisteYAbs:=[op(SpeicherlisteYAbs),Erweiterung(Unten,Oben,f,g,[-1,VorgegebeneKnoten[î],1],î-1)-TatsächlicherWert]; #Bestimmen des absoluten Fehlers von f für n=î
      if abs(TatsächlicherWert) > 10^(-10) then #Bestimmen des relativen Fehlers von f1 falls                                                          #dieser definiert ist
      SpeicherlisteYRel:=[op(SpeicherlisteYRel),abs(SpeicherlisteYAbs[-1]/TatsächlicherWert)];
    end if;
  end if;
end do;
if numelems(SpeicherlisteX)>0 then
  SpeichervektorX:=Vector[row](numelems(SpeicherlisteX),SpeicherlisteX);
  SpeichervektorYAbs:=Vector[row](numelems(SpeicherlisteYAbs),SpeicherlisteYAbs);
  PunkteAbsPY:= plot(SpeichervektorX,SpeichervektorYAbs,style = point, color=purple, legend = ["PY"]);
    #  Generierung des Punktgraphen, der sich aus den absoluten Fehlern von f ergibt
  if numelems(SpeicherlisteX)>1 then
    InterpolationsfunktionAbs:=Spline(SpeichervektorX,SpeichervektorYAbs,n);
      #  Splines aus Stützpunkten, die sich aus den absoluten Fehlern von f ergeben
    GraphAbsPY:= plot(InterpolationsfunktionAbs, n=nUnten..nOben, color=purple);
      #  Generierung des Graphen, der sich aus dem Spline aus den absoluten Fehlern von f ergibt
  end if;
end if;

if abs(TatsächlicherWert) > 10^(-10) then
  # falls der relative Fehler definiert ist analoges Vorgehen für die relativen Fehler
  if numelems(SpeicherlisteX)>0 then
    SpeichervektorYRel:=Vector[row](numelems(SpeicherlisteYRel),SpeicherlisteYRel);
    PunkteRelPY:= plot(SpeichervektorX,SpeichervektorYRel,style = point, color=purple, legend = ["PY"]);

    if numelems (SpeicherlisteX)>1 then
      InterpolationsfunktionRel:=Spline(SpeichervektorX,SpeichervektorYRel,n);
      GraphRelPY:= plot(InterpolationsfunktionRel, n=nUnten..nOben, color=purple);
    end if;
  end if;
end if;
NichtexistenzPY:=nOben-nUnten+1-numelems(SpeicherlisteX);
print(display({GraphAbsGK,PunkteAbsGK,GraphAbsPZP,PunkteAbsPZP, GraphAbsPZM,PunkteAbsPZM, GraphAbsPY,PunkteAbsPY}, title= "Absoluter Fehler", titlefont=["ROMAN",18]));
if abs(TatsächlicherWert) > 10^(-10) then
  print(display({GraphRelGK,PunkteRelGK,GraphRelPZP,PunkteRelPZP, GraphRelPZM,PunkteRelPZM, GraphRelPY,PunkteRelPY}, title= "Relativer Fehler", titlefont=["ROMAN",18]));
end if;
Text("Häufigkeit der Nichtexistenz: GK ",NichtexistenzGK, ", PZP ",NichtexistenzPZP, ", PZM ", NichtexistenzPZM, ", PY ", NichtexistenzPY);
      
 end proc

An example of how it should not look like is this:

Plotting(-1,1,2*x^2+2,1,3,10)

On a side note, Maple's return is

Warning, `GraphRelGK` is implicitly declared local to procedure `Plotting`
Warning, `GraphRelPZP` is implicitly declared local to procedure `Plotting`
Warning, `GraphRelPZM` is implicitly declared local to procedure `Plotting`
Warning, `GraphRelPY` is implicitly declared local to procedure `Plotting`

even though I did declare them.

Any suggestions on that minor issue? And on how to construct a math table which allows for a symbol like - for nonexistence?

Thank you in advance!

Hi,

I would be really grateful if someone can help me in solving the below attached problem in maple.

Thanks in advance.

Hi,

I was trying to differentiate an equation with two random variables with respect to my decision variable. But it took 4-5 hrs still MAPLE could not evaluate nor show an error.

 Later I want to find q* from the FOC and wish to substitute it back into the equation. I will be grateful if someone please help me in doing both.

doubt_1.mw

Thanks in Advance

Hi, 
I want to solve a simple newsvendor problem for general Demand Probability Distribution with maple code.
let, 
Underage cost 
                              c[u]

Overage cost 
                              c[o]

Demand (
D;
) follows pdf 
f(.);
 and CDF 
F;

(.);

Decision varible is quantity 
q;

So cost function is 

TC := piecewise(q-D >= 0, c[o]*(q-D), q-D < 0, c[u]*(-q+D));
 piecewise(0 <= q - D, c[o] (q - D), q - D < 0, c[u] (-q + D))

After taking expectation of TC and derivating it w.r.t q we will get 
q;
* as critical fractile i.e.

q*= F^(-1) (
                             c[u]    
                          -----------
                          c[u] + c[o]
)      i.e. F inverse (
c[u]/(c[u]+c[o]);
)

Can someone please help me how to find q* in MAPLE using solve() or any other function?

Thanks in Advance!

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!

Please I have an issue with the attached plot code. Can you kindly help to correct it? 

restart:
interface(rtablesize=infinity):
B:=<<0,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.49,"0","0.05","0.1","0.15","0.2","0.25","0.3","0.35","0.4","0.45","0.49">|
	<14.73,14.4,14,13.4,12.67,11.67,10.4,8.67,6,3,0,"14.73","14.4","14","13.4","12.67","11.67","10.4","8.67","6","3","0">|
     <-0.007072,0.013309,0.033707,0.054125,0.074571,0.095056,0.115597,0.136218,0.156956,0.177867,0.199036,0.22059,0.242719,0.265702,0.289932,0.31592,0.344214,0.375124,0.408175,0.441761,0.473484,0.501857>|
	<1.34E+01,1.33E+01,1.33E+01,1.32E+01,1.31E+01,1.31E+01,1.30E+01,1.29E+01,1.28E+01,1.26E+01,1.25E+01,1.22E+01,1.19E+01,1.14E+01,1.08E+01,9.86E+00,8.58E+00,6.90E+00,4.90E+00,2.81E+00,1.00E+00,-2.86E-01>>:

B:
 
 plot([B[..,[1, 2]],B[1..1,[1, 2]], B[.., [3, 4]],B[1..1,[3, 4]], B[..,[5, 6]],B[1..1,[5, 6]],B[.., [7, 8]],B[1..1,[7, 8]],
 	  B[..,[9, 10]],B[1..1,[9, 10]], B[.., [11, 12]],B[1..1,[11, 12]],B[..,[13, 14]],B[1..1,[13, 14]],B[.., [15, 16]],B[1..1,[15, 16]],
 	  B[..,[17, 18]],B[1..1,[17, 18]], B[.., [19, 20]],B[1..1,[19, 20]],B[..,[21, 22]],B[1..1,[21, 22]]],
 	  legend = ["","Experimental","","Simulation"],
 	  style = ["line","line","line","line","line","line","line","line","line","line","line","line","line","line","line","line",
 	 		"line","line","line","line","line","line"],
 	  color=[blue,red], labels=[`V (V)`, `Jsc (mA/cm^2)`]);
 

Download Graph_Example.mw

Hi there!

I am trying to tell Maple to print something one time once a condition has been met for a variable that ranges from a to b, let's say from 1 to 9 in this case. Testbestanden is to be printed, if (F[i]+F[i+1]<20) is true for every i, and Durchgefallen is to be printed as soon as one counterexample is found. The following code illustrates what I need, it doesn't work, however. What is the correct syntax? Switching the lines "if (F[i]+F[i+1]<20)" and "for i from 1 to 9" leads to Maple printing "Testbestanden" for every i that fulfils the statement which is not what I need either.

for i from 1 to 10 do
  F[i]:=i
end do;

if (F[i]+F[i+1]<20)
for i from 1 to 9
then print(Testbestanden)

else print(Durchgefallen)
end if
end do

Thank you in advance! :)

Hello
I solved an equation using Maple . but in answer there is RootOf statement that I can't change to a simple statement without RootOf.

Download spr.mw

How can I simplify to an answer without RootOf?

Thank u

Hi
I want to say maple convert this radicals to a single radical expr.

convert to :

thank u .

In the following code, I want the base of the logarithm in the label to appear at the conventional position not beside the logarithm. How can I achieve this? Thank you all in anticipation of your educative response and do stay safe.

restart;

B:=<<39,76,151,301,601>|<7.71E-8,5.43E-9,3.55E-10,2.11E-11,1.32E-12>|
	<26,51,101,201,401>|<6.46E-3,1.17E-4,1.88E-6,2.96E-8,4.46E-10>|
	<26,51,101,201,401>|<2.74E-4,6.34E-6,1.16E-7,1.85E-9,2.92E-11>|
	<26,51,101,201,401>|<6.48E-4,4.39E-5,2.99E-6,1.88E-7,1.18E-8>>:

for i from 1 to 5 do
   B[i, 2] := log[10](B[i, 2]):		  			
   B[i, 4] := log[10](B[i, 4]): 		
   B[i, 6] := log[10](B[i, 6]): 		
   B[i, 8] := log[10](B[i, 8]):	
       

   
end do:  # computing the log of the max-error
B: # This is the table of values we'll plot.

T:=plot([B[..,[1, 2]],B[1..1,[1, 2]], B[.., [3, 4]],B[1..1,[3, 4]], 
	 B[..,[5, 6]],B[1..1,[5, 6]],B[.., [7, 8]],B[1..1,[7, 8]]], 
	legend = ["","BFFM","", "BHT","", "BHTRKNM", "", "BNM"],
	#title="Efficiency Curve for Example 1",
	style = ["pointline","point","pointline","point","pointline","point","pointline","point"], 
	symbolsize = 15,axes = framed, 
	symbol = [box,box, circle,circle,diamond,diamond,solidcircle, solidcircle],
	color=[red, red, blue,blue, black, black, green, green], 
	axis = [gridlines = [colour = green, majorlines = 1,linestyle = dot]], 
	labels = ["NFE", "log[10](Max Err)"]);

I have written the following attached code to use Euler explicit method to solve the following IVP

diff(y(x), x) = 2*(1+x)-y(x), y(2) = 5
With Exact solution  y(x) = 2*x+exp(-x)/exp(-2)

However, I found out that my exact results are not correct while the numerical results are okay. What have I done wrong in the code? Can someone modify the code?

Thank you and kind regards.

Hi all, we know Maple provided discrim method to find discriminant of a polynomial 

I want to write a similar method with independent variable is ,... my code is below

Some examples

Maple already similar method? If not have, we can improve performance it?

Thank you very much.

Hi,

I have a random variable that follows Uniform(1,4). Now I have a function which is of the following type:

g := a*alpha+b*t/alpha+exp(alpha)

where,

A := RandomVariable(Uniform(c, d));
                 RandomVariable(Uniform(c, d))
f := proc (alpha) options operator, arrow; PDF(A, alpha) end proc;
alpha -> PDF(A, alpha)
#Defining expectation fuction
E := proc (alpha) options operator, arrow; int(alpha*f(alpha), alpha = c .. d) end proc;
alpha -> int(alpha f(alpha), alpha = c .. d)
#g is a function of random variable &alpha;, where a and b are parameter

now I want to find the expectation of g and the first derivative of expectation of g,

E(g)

diff(E(g), t)

I understand the way I have defined E(alpha) is improper. But please understand my intent and help! here is the maple file also doubt_1.mw

Dear all,

Please I want only 8 points to show on this curve, how do I specify it?

plot(ln(1+sin(Pi*x)), x = 0 .. 1, legend = numerical, style = point, symbol = box, color = blue, symbolsize = 15, numpoints = 8);

Thank you all and kind regards.

Please do keep safe amidst this global pandemic.