Potenzreihen sind nicht sonderlich gut geeignet für sowas weil man
schnell Überläufe aus dem Wertebereich bekommt auch wenn der Endwert im
Rechenbereich liegt.
Beispiel: exp(-10) ≈ 0, hier werden die Summanden x^k/k! in der
Reihenentwicklung recht groß und neben dem Überlauf-Problem handelt man
sich noch Ungenauigkeit durch Auslöschung ein.
http://de.wikipedia.org/wiki/Auslöschung_(numerische_Mathematik%29
Das Überlauf-Problem bekommt man in den Griff, indem man ein Polynam
anders darstellt und auswertet, zB zur Basis von Bernsteinpolynomen. Im
Endeffekt führt das zu Bézierkurven, die besser bekannt sein dürften.
Für exp(-t), t \in {0, 4/3}
Kann man die Funktion zB so berechnen:1 | #include <stdio.h>
|
2 | #include <math.h>
|
3 |
|
4 | double y[] = { 1, 0.56, 0.385, 0.2636 };
|
5 |
|
6 | double castelj (double t0)
|
7 | {
|
8 | double t = t0*3.0/4.0; // t in {0,1}
|
9 |
|
10 | double y00 = y[0] + t * (y[1] - y[0]);
|
11 | double y01 = y[1] + t * (y[2] - y[1]);
|
12 | double y02 = y[2] + t * (y[3] - y[2]);
|
13 |
|
14 | double y10 = y00 + t * (y01 - y00);
|
15 | double y11 = y01 + t * (y02 - y01);
|
16 |
|
17 | double y20 = y10 + t * (y11 - y10);
|
18 |
|
19 | printf ("y[%f] = %f, %f\n", t0, y20, exp(-t0)-y20);
|
20 | }
|
21 |
|
22 | int main()
|
23 | {
|
24 | castelj (0);
|
25 | castelj (0.25);
|
26 | castelj (0.5);
|
27 | castelj (0.75);
|
28 | castelj (1);
|
29 | castelj (1.25);
|
30 | castelj (4.0/3.0);
|
Die Ausgabe des Programms ist
1 | y[0.000000] = 1.000000, 0.000000
|
2 | y[0.250000] = 0.779056, -0.000255
|
3 | y[0.500000] = 0.605649, 0.000882
|
4 | y[0.750000] = 0.471418, 0.000948
|
5 | y[1.000000] = 0.368003, -0.000124
|
6 | y[1.250000] = 0.287042, -0.000537
|
7 | y[1.333333] = 0.263600, -0.000003
|
Am interessantesten hier ist die letzte Spalte, die den absoluten Fehler
der Näherung vom Grade 3 zu exp(-t) angibt.
Wichtig sind eigentlich nur die Auswertung (kein Hornerschema etc!)
und die Konstanten y[].
Die Umstellung auf Festpunkt-Arithmetik ist dann nur noch ein bisschen
Technik :-)
