Genau genommen: Taylorentwicklung für Sinus und Cosinusfunktionen
Um numerische Aufgaben zu lösen, ist es oft nötig, bis auf die letzte Stelle zu rechnen, um Folgefehler so gering wie möglich zu halten. Dies nützt allerdings nichts, wenn die implementierten Funktionen schon Fehler liefern, die 3 oder 4 10er Potenzen unterhalb der Rechnergenauigkeit liegen, wie es bei den Sinus-und Cosinusfunktionen von Pascal + der Fall ist.
Einen Anhaltspunkt für die Ungenauigkeit hat man nur an den Stellen 0, pi/2, pi... (z. B. sin(pi/2) = 1.0000000032). Um im gesamten Intervall die Funktionen zu testen, ohne einen genauen Wert zu haben, kann man den Einheitskreis verwenden. Flier bilden der Sinus und der Cosinus zusammen mit dem Radius ein rechtwinkliges Dreieck. Nach dem Satz von Pythagoras folgt:
r = sqrt(sin^2 + cos^2)
Da im Einheitskreis der Radius die Länge 1 hat, erhält man mit 1—r einen Fehler, der indirekt etwas über die Genauigkeit der Sinus- und der Cosinusfunktion aussagt. Dies ist im Programm SINTEST.PAS dargestellt. (Natürlich sollte die Wurzelfunktion auch auf ihre Genauigkeit überprüft werden). Nun zur Behebung dieses Mangels.
Aufbereiten des Arguments x: Bei einem Aufruf der cos(x) Funktion wird zu x pi/2 hinzugezählt und sin(x) aufgerufen, da die Cosinusfunktion nur eine um pi/2 verschobene Sinusfunktion ist. Bei einem Aufruf der sin(x) Funktion wird der x-Wert, der beliebig reell sein kann, auf einen Bereich von 0 bis pi/2 reduziert, und das Vorzeichen aufbewahrt. Dieses Intervall genügt vollkommen, um sämtliche Funktionswerte zu ermitteln.
Berechnung von sin(x) im Intervall [0..pi/2]: Die Funktion sns berechnet die Funktionswerte von Sinus durch Taylorreihenentwicklung für Sinus und Cosinus.
Diese lautet für sin x:
x - x^3/3! + x^5/5! - x^7/7! + - ...
und für cos x:
1 - x^2/2! + x^4/4! - x^6/6! + - ...
Da die einzelnen Glieder der Entwicklungen sehr schnell gegen Null streben, kann man den jeweils letzten berechneten Wert mit der Rechnergenauigkeit vergleichen, und hat somit ein relativ gutes Abbruchkriterium. Es wurde neben der Genauigkeit auch Wert auf eine erträgliche Geschwindigkeit gelegt, deshalb mag die Art der Programmierung etwas seltsam erscheinen. Da Taylorreihen am Intervallrand sehr viele Glieder (und auch Rechenzeit) benötigen, um die Funktion anzunähern, wurde im Bereich 0 < = x < pi/4 nach der sin Reihe entwickelt und im Bereich pi/4 < = x < pi/2 nach der cos-Reihe im Intervall [0..pi/4], Bei einer Rechnerge'nauigkeit von 1e-12 heißt das maximal 8 Glieder. Um weitere Rechenzeit einzusparen, wurden die Koeffizienten der Reihen (1/i!) ausgerechnet und in Felder abgespeichert. Um Fehler zu vermeiden, werden die einzelnen Glieder zuerst in einem weiteren Feld abgespeichert und erst am Ende von der kleinsten Zahl aufwärts zusammengezählt. Wer sich nicht vorstellen kann, welche Fehler auftreten können, teste einmal die beiden kleinen Programme: Test 1 und Test 2
Das richtige Ergebnis lautet
1.000000500el0.
Alles klar?
program test1;
var i:integer;
x:real;
begin
x: = 0;
for i: = 1 to 10000 do
x: = x + 0.05;
x:= x + 1e10;
writeln(x);
end.
program test2;
var i:integer;
x:real;
begin
x: = 1e10;
for i: = 1 to 10000 do
x: = x + 0.05;
writeln(x);
end.
test1/test2: Zwei kleine Programme zum Genauigkeitstest. Vergleichen Sie selbst.
Verwendung der Routinen: Will man nun diese Routinen verwenden, muß man nur am Anfang des Prozedurvereinbarungsteils die Datei mit (*$I sincos.pas *) einfügen. Der Pascalcompiler ersetzt dann selbsttätig die alten Funktionen durch die neuen.
program sintestO;
(**********************************************************)
(* Programm zum testen der Genauigkeit der Sinusfunkticn. *)
(**********************************************************)
var df,pi,b,x1,x2:real;
(*$I sincos.pas*)
begin
pi:=arctan(1)*4;
b:=0;
while b <= 2*pi do
begin
x1:=cos(b);
x2:=sin(b);
df:=x2-x1;
writeln(b,' r: ' ,sqrt(x1*x1+x2*x2),' Fehler:' ,1-sqrt(x1*x1+x2*x2));
b:=b+1e-2;
end;
writeln(chr(7));
end.
Listing 1: Testen der Genauigkeit
Listing 2: Taylorentwicklung in Pascal
(**********************************************************)
(* SIN und COS Routinen fuer Pascal unter Verwendung von *)
(* Taylor Reihen. *)
(* *)
(* 07.07.1987 *)
(* Peter Sadoni *)
(**********************************************************)
function sns(x:real) :real;
forward;
function sin(x:real):real; (* Einsprung fuer sin x *)
(* Reduzierung von x auf den Bereich 0..pi/2 *)
var vz:integer;
begin (* x 'liegt zwischen: *)
vz:=1; (* *)
x:=x-trunc(x/6.283185307180)*6.283185307180; (* -2*pi uni 2*pi *)
if x<0 then (* *)
x: =x+6.283185307180; (* 0 und 2*pi *)
if x>3.141592653590 then (* 0 uni pi mit *)
begin (* vz = Vorzeichen von sin*)
vz:=-1;
x: =x-3.141592653590; (* *)
end;
if x>1.570796326795 then (* *)
x:=3.141592653590 -x; (* 0 und pi/2 *)
sin:=vz*sns(x); (* Aufruf von eigentl.*)
end; (* Sinusfunkticn *)
function cos(x:real) :real; (*Einsprung fuer cos x*)
(* Unwandlung vcn Cosinus in Sinus *)
begin
x:=x-1.570796326795; (* um Pi/2 verschieben*)
cos:=-sin(x); (* Aufruf von sin x *)
end;
function sns;
(* eigentliche Sinusfunktion zur Berechnung von sin x mit 0 <= x <- pi/2 *)
const eps=1e-12; (* Rechnergenauigkeit *)
type vektor = array[1..10] of real;
var x2,hlp1,hlp2:real;
n,i:integer;
koeff,hlp3:vektor;
procedure sinbelegen(var sinkoeff:vektor);
(* Belegen vcn sinkoeff mit den Koeffizienten der Taylor-Reihe fuer Sinus *)
begin
sinkoeff[1]:=1;
sinkoeff[2]:=-0.166666666667; (*-1/6 *)
sinkoeff[3]:=8.333333333333e-3; (* 1/120 *)
sinkoeff[4]:=-1.984126984127e-4; (*-1/5040 *)
sinkoeff[5]:=2.755731922399e-6; (* 1/362880 *)
sinkoeff[6]:=-2.505210838544e-8; (*-1/39916800 *)
sinkoeff[7]:=1.605904383682e-ll; (* 1/6227020.8e3 *)
sinkoeff[8]:=>-7.647163731820e-13; (*-1/1.307674368e12 *)
end;
procedure cosbelegen(var coskoeff:vektor);
(* Belegen von coskoeff mit den Koeffizienten der Taylor-Reihe fuer Cosinus*)
begin
coskoeff[1]:=-0.5; (*-1/2 *)
coskoeff[2]:=4.166666666667e-2; (* 1/24 *)
coskoeff[3]:=-1.388888888889e-3; (*-1/720 *)
coskoeff[4]:=2.480158730159e-5; (* 1/40320 *)
coskoeff[5]:=-2.755731922399e-7; (*-1/3628800 *)
coskoeff[6]:=2.087675698787e-9; (* 1/479001600 *)
coskoeff[7]:=-1.147074559773e-ll; (*-1/8.71782912e10 *)
coskoeff[8]:=4.779477332387e-14; (* 1/2.0922789888e13 *)
end;
begin
n:=1;
if x>=0.785398163397 then (* wann x aus pi/4..pi/2 *)
begin (* dann: *)
cosbelegen(koeff); (* Cosinus Reihenentwicklg *)
x:=1.570796326795-x; (* mit x=pi/2-x *)
x2:=x*x; (* x'2 einmal, berechnen *)
hlp1:=x2; (* Variable fuer x^n *)
hlp2:=1; (* Variable fuer Wert des *)
(* Polynoms *)
end
else (* wenn x aus 0..pi/4 *)
begin (* dann: *)
sinbelegen (koeff); (* Sinus Reihenentwicklung*)
x2:=x*x; (* x^2 einmal berechnen *)
hlp1:=x; (* Variable fuer x^n *)
hlp2:=0; (* Variable fuer Wert des *)
(* Polynoms *)
end;
hlp3[n]:=koeff[n]*hlp1; (* 1. Glied der Bitwicklg *)
while abs(hlp3[n]) >= eps do (* Abbruchbedingung *)
begin
n:=n+1; (* naechstes Glied *)
hlp1:=hlp1*x2; (* ermitteln *)
hlp3[n]:=koeff[n]*hlp1; (* *)
end;
for i:=n-1 downto 1 do (* Addition der einzelnen *)
hlp2:=hlp2+hlp3[i]; (* Glieder *)
sns:=hlp2; (* Rueckgabe des *)
end; (* Näherungswertes *)
Peter Sadoni