Romberg-Integration<Text-field layout="Heading 1" style="Heading 1">Matrix</Text-field>restart:Funktion eingebenf:=(x)->sin(x);Grenzena:=0;
b:=Pi;Matrixgr\366ssen:=6;BerechnungR:=Matrix(n,n):## AusgangsWert R[1,1]
R[1, 1]:=(b - a) / 2 * (f(a) + f(b)):
## Trapezregel anwenden (die Werte in der ersten Spalte entsprechen der Trapezregel)
for k from 2 to n by 1 do
h:=(b - a) / (2^(k - 1));
R[k, 1]:= 1 / 2 * (R[k - 1, 1] + h * 2 * (sum(f(a + (2 * i - 1)* h), i=1..2^(k - 2))));
end do:
## Richardson Extrapolation (alle Werte ab der ersten Spalte wurden durch dieses Beschleunigungsverfahren erzeugt!)
for j from 2 to n by 1 do
for k from j to n by 1 do
R[k,j] := R[k, j-1] + (R[k,j-1] - R[k-1, j-1])/(4^(j-1)-1);
end do:
end do:Ausgabe (Resultat steht unten rechts in der Matrix)evalf(R);<Text-field layout="Heading 1" style="Heading 1">Alternativ mit Fehlerschranke</Text-field>restart:
Digits:=20:Funktion eingebenf:=(x)->sin(x);
Grenzena:=0;
b:=Pi;Fehlereps:=10^(-10);Matrixgr\366ssen:=10:
Berechnung & Ausgabe (beste Ann\344herung zuunterst)err:=100:
R:=Matrix(n,n):
R[1,1]:=(b-a)/2 *(f(a)+f(b)):
for k from 2 by 1 while ( err > eps ) do
h:= (b-a)/(2^(k-1)) :
R[k,1]:= 1/2*(R[k-1,1]+h*2*(sum(f(a+(2*i-1)*h), i=1..2^(k-2)))) :
for j from 2 by 1 while j <= k do
R[k,j]:= R[k,j-1] + (R[k,j-1] - R[k-1,j-1])/(4^(j-1)-1) :
end do:
err:= abs( evalf( R[k,k] ) - evalf( R[k-1,k-1] ) );
printf("%2.1d. Romberg-Approximation :", k);
printf("%20.15f\n", evalf(R[k, k]));
end do: