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!
