The code with proper format: :)
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