 |
von = getlong();
printf("\nGeben Sie die obere Intervallgrenze ein:");
bis = getlong();
n = von;
while(n <= bis)
{
if (is_prim(n))
printf(" %ld\n", n);
n++;
}
}
int is_prim(long int i)
/* liefert TRUE, falls i eine Primzahl ist */
{
long int i;
long int j;
if (i < 2) return(FALSE);
j = 2;
while (j*j <= i)
{
if (i % j == 0)
return (FALSE);
else if (j == 2)
j = 3;
else
j = j + 2;
}
return (TRUE);
}
void primfaktoren(long int i)
/* Druckt Primfaktoren einer Zahl */
{
long int i;
long int j = 2;
while (j*j <= i)
{
while ((i%j) != 0)
{
if (j == 2)
j = 3;
else
j = j + 2;
}
i = i/j;
printf(" %ld", j);
}
if (i > 1)
printf(" %ld", i);
printf("\n");
}
long int getlong ()
/* Ziffernfolge lesen und in long konvertieren */
{
long int i = 0;
int c;
c = getchar();
while (isdigit(c))
{
i = 10*i + c - '0';
c = getchar();
}
return(i);
}
Auflösen arithmetischer Ausdrücke
U. a. ist die Aufgabe von Compilern das Aufbrechen von Formeln in eine für die
lineare Maschinensprache geeignete klammerlose Form (Postfix-Notation). Auch die direkte
Interpretation einer Formeleingabe ist eine häfig gestellte aufgabe, z.B. bei
Kalkulationsprogrammen oder Formelinterpretern. Wie ein Arithmetischer Ausdruck aufgebaut
ist, kann mit dem Schon bekannten Syntaxdiagrammen dargestellt werden - oder man verwendet
die Textvariante, die Backus-Naur-Darstellung. Beim vorgesehenen Formelinterpreter
sollen die vier Grundrechenarten sowie Klammerung realisiert werden. Natürlich gilt
auch hier "Punkt vor Strich". Eine Formel kann dann folgendermaßen definiert werden:
Ausdruck ::= Term { + Term | - Term }
Term ::= Faktor { * Faktor | / Faktor }
Faktor ::= Zahl | ( Ausdruck )
Zahl ::= [-] Ziffer {Ziffer} [. {Ziffer}]
Ziffer ::= 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
Dabei bedeutet:
- ::=: "wird definiert als"
- { ... }: Was in geschweiften Klammern steht, darf beliebig oft
wiederholt werden.
- [ ... ]: Was in eckigen Klammern steht, darf auch fehlen.
- |: Alternative Auswahl, z. B. darf Faktor entweder eine Zahl oder
eine Ausdruck in runden Klammern sein.
Die Ziffern sowie die Zeichen '(', ')', '+', '-', '*', '/' sind Terminalsymbole.
Zur Programmierung: Die folgenden Funktionen geben genau die Definition wieder
für Ausdruck, Term und Faktor müssen jeweils die Prozeduren definiert werden,
die sich dann gegenseitig rekursiv aufrufen.
Es muß immer nur das nächste Eingabezeichen betrachtet werden, um die Aufgabe zu
lösen daher reicht getchar() aus. Das gelesene Zeichen c ist global definiert.
Die folgenden Ansätze sind in einer Art "Pseudocode" skizziert, und ohne Deklarationen.
Ausdruck:
{
wert = term;
solange (c in ['+','-'])
{
falls (c = '+') wert = wert + term;
falls (c = '-') Wert = wert term;
}
Ausdruck = wert;
}
Term:
{
wert = faktor;
solange (c in ['*','/'])
{
falls (c = '*') wert = wert * faktor;
falls (c = '/') Wert = wert / faktor;
}
Term = wert;
}
Faktor:
{
Lies ein Zeichen c von der Eingabe;
wenn (c in ['0'..'9']) wert = Zahl;
wenn (c = '(')
{
wert = Ausdruck;
Lies schlieîende Klammer;
}
Faktor = Wert;
}
Zahl:
{
Lies eine Zahl mit Vorzeichen zeichenweise ein;
}
Damit das Bearbeiten eines Ausdrucks nicht nur von der Eingabe, sondern auch z. B.
aus einem String heraus leicht möglich ist, definieren wir noch eine Funktion
int lies() die dem Einlesen eines Zeichens dient.
#include <stdio.h>
#define FALSE 0
#define TRUE 1
int c; /* eingelesenes Zeichen */
float wert; /* Ergebniswert */
/* Funktionen-Vorwaertsdeklaration */
float ausdruck();
float term();
float faktor();
float zahl();
int lies();
main()
{
do
{
printf("Eingabe : "); /* Eingabeaufforderung */
wert = ausdruck(); /* Einlesen und Berechnen */
printf("Ergebnis: %0.5f\n\n",wert); /* 5 Nachkommastellen */
}
while (c != EOF); /* Ende mit Ctrl-D oder Ctrl-Z */
}
int lies()
/* Einlesen eines Zeichens, Leerzeichen übergehen */
{
while((c = getchar()) == ' ');
return(c);
}
float ausdruck()
/* Ausdruck ::= Term { +Term || -Term } */
{
float wert;
wert = term();
while (c == '+' || c == '-')
{
if (c == '+') wert = wert + term();
if (c == '-') wert = wert - term();
}
return(wert);
}
float term()
/* Term ::= Faktor { *Faktor || /Faktor } */
{
float wert;
wert = faktor();
while (c == '*' || c == '/')
{
if (c == '*') wert = wert * faktor();
if (c == '/') wert = wert / faktor();
}
return(wert);
}
float faktor()
/* Faktor :== Zahl || ( Ausdruck ) */
{
float wert;
c = lies(); /* erst hier wird gelesen */
if (('0' <= c && c <= '9') || (c == '-'))
wert = zahl();
if (c == '(') /* rekursiver Aufruf! */
{
wert = ausdruck();
c = lies(); /* ')' wird ueberlesen */
}
return(wert);
}
float zahl()
/* Zahl ::= [-] Ziffer {Ziffer} [.{Ziffer}] */
{
float wert, teiler;
int negativ = FALSE; /* Merker fr neg. Vorzeichen */
int bruch = FALSE; /* Merker fr Dezimalpunkt */
wert = 0.0; teiler = 1.0;
if (c == '-') /* Vorzeichen merken */
{
negativ = TRUE;
c = lies();
}
while (('0' <= c && c <= '9') || (c == '.'))
{ /* Zahl konvertieren */
if (c == '.') /* Dezimalpunkt */
bruch = TRUE;
else
{
wert = wert * 10.0 + (c - '0'); /* weitere Stelle dazu */
if (bruch) teiler = teiler*10.0; /* fr Nachkommastellen */
}
c = lies(); /* n„chstes Zeichen */
}
if (negativ) wert = -wert;
if (bruch) /* nur dann ist teiler != 0 */
return(wert/teiler); /* Dezimalpunkt beachten */
else
return(wert);
}
Nullstellenberechnung
Die Nullstelle einer Funktion f(x) soll näherungsweise bestimmt werden.
Das Aufsuchen der Nullstelle einer Funktion f(x) entspricht der Lösung
der Gleichung f(x) = O. Ist das Intervall [a, b] bekannt, in dem
sich die Nullstelle befindet, so kann dieses halbiert werden und aufgrund des
Funktionswertes an der Halbierungsstelle überprüft werden, ob sich die
Nullstelle im linken oder im rechten Teilintervall befindet. Dieses
Intervall wird nach dem gleichen Verfahren wieder halbiert. Dieses
als Bisektion bezeichnete Verfahren wird solange wiederholt, bis das
Intervall genügend klein ist.
Zur Programmierung dieses Verfahrens ist es günstig, eine Funktion zu
vereinbaren, die als Ergebnis den Näherungswert der Nullstelle einer
Funktion float f() liefert:
float nullstelle(float a, float b)
{
float X;
do
{
X = (a + b )/2.0;
if (f(a)*f(b) > 0)
a = X;
else
b = X;
}
while ((b - a) > EPS);
return(X);
Die Parameter a und b werden hier als lokale Hilfsvariable verwendet.
Mit EPS wird dabei die Genauigkeit bezeichnet, ein Wert, der
nahe beim kleinsten darstellbaren Wert für float liegt. Definition
beispielsweise über #define EPS 1.0-10.
Anstatt das Intervall in jedem Iterationsschritt zu halbieren, kann es
auch im Verhältnis der Funktionswerte an den Intervallgrenzen geteilt
werden. Dieses als Regula Falsi bezeichnete Verfahren führt im allgemeinen
rascher zum Ziel. In der obigen Funktion muß zu diesem Zweck die Anweisung
X = (a + b)/2.0;
durch
X = a - f(a)*(b-a)/(f(b)-f(a))
ersetzt werden. Um die wiederholte Berechnung der Funktionswerte f(a) und f(b)
zu vermeiden, wurden im nachfolgenden Programmbeispiel Hilfsvariable fa, fb
und fx verwendet.
float nullstelle(float a, float b)
{
float X, fa, fb, fx;
fa = f(a);
fb = f(b);
do
{
X = a - fa*(b - a)/(fb - fa);
fx = f(X);
if (fa*fx > 0)
a = X; fa = fx;
else
b = X; fb = fx;
}
while ((b - a) > EPS);
return(X);
Das Verfahren von NEWTON erlaubt die näherungsweise Bestimmung der Nullstelle
einer Funktion, deren Ableitung bekannt ist. Von einem Näherungswert
X0 ausgehend, wird eine verbesserte Näherung X1 bestimmt,
indem die Tangente an die Funktion an der Stelle f(X0) mit der
X-Achse geschitten wird:
Die verbesserte Näherung X1 berechnet sich zu
X1 = X0 - f(X0)/f'X0
Die Vorschrift wird solange wiederholt, bis sich keine nennenswerte Verbesserung
mehr ergibt. Man kann die Ableitung der Funktion annähern, indem man den
Differenzenquotienten verwendet (df = (f(X+dH) - f(X))/dH;). Bei Funktionen,
deren Ableitungsfunktion bekannt ist, kann man bei der Zuweisung an Fs
auch der Funktionswert der Ableitungsfunktion berechnet werden.
double Newton (double X0, double Eps, int ItMax)
/* X0: Startwert der Nullstellenbestimmung */
/* Eps: Genauigkeit der Nullstellenbestimmung */
/* ItMax: Maximale Anzahl der Iterationen */
{
double Xi, Fi, Fs;
double dH = 1.e-6;
int it = 0;
Xi = X0; /* Startwert */
Fi = f(X0); /* Funktionswert bei X0 */
Fs = (f(X0+dH) - f(X0))/dH; /* numerische Ableitung f'(X0) */
do
{
/* verschwindende Ableitung abfangen */
if (fabs(Fs) < Eps)
return Xi;
/* Begrenzung der Iterationen */
if (it > ItMax)
return Xi;
/* Genauigkeit ueberpruefen */
if (fabs(Fi) < Eps)
return Xi;
/* Berechnung des neuen Funktionswerts */
Fi = f(Xi);
Fs = (f(Xi+dH) - f(Xi))/dH; /* numerische Ableitung */
Xi = Xi - Fi/Fs;
it++;
}
while (1);
return 0.0;
}
Numerische Integration
Zur numerischen Integration einer Funktion kann diese in kleine Intervalle
unterteilt werden. Innerhalb der Intervalle vird der Funktionsgraph durch
Geraden angenähert.
Der Wert des Integrals wird durch die Summe der Tapezfl„chen angenähert:
Zur Programmierung dieses Näherungsverfahrens wird eine FUnktion verwendet,
bei der die Intervallgrenzen sowie die Anzahl der Teilintervalle vorgegeben
werden. Die zu integrierende Funktion wird definiert als Funktion float f().
float integral(float a, float b, int anz)
{
float sum, h;
int i;
h = (b - a)/anz;
sum = (f(a) + f(b))/2.0;
for(i = 1; i < anz; i++)
sum = sum + f(a + float(i)*h);
return(h * sum);
}
Genauere Resultate erzielt man, ween beim Integreren nicht Geraden, sondern
Parabelstücke verwendet werden. Die von SIMPSON gefundene Näherungsformel
lautet dann
Die folgende Programmvariante verwendet die Simpson-Formel:
float integral(float a, float b, int anz)
{
float sum, h;
int i, k;
h = (b - a)/anz;
k = 2;
sum = f(a) + f(b);
for(i = 1; i < anz; i++)
{
k = 6 - k; /* 4, 2, 4, 2, ... */
sum = sum + float(k)*f(a + float(i)*h);
}
return(h/3.0 * sum);
}
Gausselimination (lin. Gleichungen)
Eine der ältesten direkten Methoden zur Lösung von linearen
Gleichungssystemen ist die Gauss-Elimination. Gegeben ist ein Gleichungssystem
Ax = b
Wird die k-te Gleichung nach der Unbekannten xk aufgelöst, gilt
für akk ungleich 0:
Bei der Gaußschen Eliminationsmethode werden die Ausgangsgleichungen so
manipuliert, daß die Koeffizientenmatrix in ihrer Hauptdiagonalen den
Wert 1 und unterhalb der Hauptdiagonalen nur Nullen enthält.
Zwei Grundtypen von Matrixoperationen werden in der Gaußsehen
Eliminationsmethode benutzt: skalare Multiplikation und Addition. Jede
Gleichung kann mit einer Konstanten multipliziert werden,
ohne das Ergebnis zu verändern. Dies ist äquivalent mit der
Multiplikation einer Zeile der Koeffizientenmatrix und des zugehörigen
Elementes im Konstantenvektor mit demselben Wert. Jede
Gleichung kann auch durch die Summe zweier Gleichungen ersetzt werden.
Die folgenden Schritte werden die Benutzung der Gaußschen Ellmination anhand
der folgenden Gleichungen erläutert:
Koeffizienten- Lösungs-
Matrix A Vektor b
13 -8 -3 20
-8 10 -1 -5
-3 -1 11 0
Ziel der Elimination ist es, aus Ax = b eine äquivalentes System
Cx = d mit C als obere Dreiecksmatix zu erzeugen:
Dieses System kann dann durch Rückwärtseinseten gelöst werden:
Rechnen wir das Beispiel einmal durch.
Die erste Variable wird von allen außer der ersten Gleichung eliminiert.
Ziel dieser Manipulation ist es, den Wert 1 als oberstes Element der ersten Spalte
zu erzeugen. Die übrigen Positeonen der Spalte werden 0. Dazu wird die gesamte
erste Zeile durch das erste Element in der Zeile (das Pivot-Element) dividiert.
Dies erzeugt den Wert 1 in der ersten Positeon der Diagonalen. Die erste Gleichung
erhält dann folgendes Aussehen:
[ 1 -0.61 -0.23] [1.5]
Die erste Unbekannte wird von der zweiten Zeile eliminiert, indem die neue erste
Zeile mit dem ersten Element in der zweiten Zelle multipliziert und dann von der
zweiten Zeile subtrahiert wird. Als neue zweite Zeile ergibt sich dann:
Zeile 2: [ -8.00 10.00 -1.00 ] [- 5]
Minus Zeile 1 mal -8: [ -8.00 4.88 1.84 ] [ 12]
-----------------------------
Ergibt: [ 0.00 5.12 0.84 ] [ 7]
In ähnlicher Weise eliminiert man die erste Variable aus der dritten Gleichung.
Die drei Gleichungen haben nun folgendes Aussehen:
[ 1.00 -0.61 -0.23 ] [1.5]
[ 0 5.12 2.84 ] [ 7]
[ 0 2.83 10.31 ] [4.6]
Im nächsten Schritt wird der Wert 1 in der zweiten Positeon der
zweiten Zeile (dem neuen Pivot-Element) erzeugt. Die zweite
Zeile wird dazu durch das zweite Element dividiert. Die zweite
Variable wird von der dritten Gleichung eliminiert, indem eine
Null in der zweiten Positeon direkt unter dem Pivot-Element
erzeugt wird. Die drei Gleichungen sehen dann so aus:
[ 1.00 -0.61 -0.23 ] [1.5]
[ 0 1.00 -0.56 ] [1.4]
[ 0 0 8.7 ] [8.7]
Der letzte Schritt dieser Phase besteht darin, den Wert 1 in der dritten
Pivot-Positeon zu erhalten. Dies erreicht man durch Division der dritten
Gleichung durch den Pivot-Wert. Insgesamt erhält man dann:
[ 1.00 -0.61 -0.23 ] [1.5]
[ 0 1.00 -0.56 ] [1.4]
[ 0 0 1.00 ] [1.0]
Das Ergebnis entspricht den drei Gleichungen:
x - 0.6ly - 0.23z = 1.5
y - 0.56z = 1.4
z = 1
Die dritte Gleichung kann direkt gelöst werden, da sie nur eine
Unbekannte hat. Das Ergebnis lautet:
z = 1
Die zweite Gleichung kann folgendermaßen umgeformt werden:
y = 1.4 + 0.56z
Durch Einsetzung des Wertes von z ergibt sich für y der Wert 2.
Die Einsetzung von y und z in die erste Gleichung liefert für x den
Wert 3. Diese Phase der Berechnung Wird als "Rücksubstitution
bezeichnet. Das Programm stellt sich dann so dar:
#include <stdio.h>
#define N 3 /* Spalten/Zeilen der Matrix */
/* = Anzahl der Gleichungen */
int main(void)
{
int i, j, k; /* Schleifenzaehler */
double f; /* Hilfsvariable */
double a[N][N]; /* Koeffizientenmatrix */
double b[N]; /* rechte Seite */
double x[N]; /* Loesung */
/* Einlesen */
for(i = 1; i <= N; i++)
{
for(j = 1; j <= N; j++)
scanf("%lg", &a[i-1][j-1]);
scanf("%lg", &b[i-1]);
}
/* Elimination */
for(i = 1; i <= (N-1); i++)
{
for(k = i+1; k <= N; k++)
{
f = a[k-1][i-1] / a[i-1][i-1];
for(j = i+1; j <= N; j++)
a[k-1][j-1] -= f*a[i-1][j-1];
b[k-1] -= f*b[i-1];
}
}
/* Ruecksubstitution */
for(i = N; i >= 1; i--)
{
for(j = i+1; j <= N; j++)
b[i-1] -= a[i-1][j-1] * x[j-1];
x[i-1] = b[i-1] / a[i-1][i-1];
}
/* Ausgabe */
for(i = 1; i <= N; i++)
printf("x%d = %g\n", i, x[i-1]);
return 0;
}
Ein Problem ergibt sich, wenn die Pivot-Elemente relativ klein sind, weil dann
- bedingt durch die begrenzte Rechengenauigkeit - sich die Rundungsfehler relativ
stark auf das Endergebnis auswirken.
Die Genauigkeit der Gaußschen Eliminationsmethode kann dadurch
verbessert werden, daß man zwei Zeilen derart vertauscht,
daß das Element mit dem größten absoluten Wert das
Pivot-Element wird. Nehmen wir z. B. an, daß die vorhergehenden drei
Gleichungen ursprünglich in einer anderen Reihenfolge geschrieben
worden wären:
[ -3 -1 11 ] [ 0 ]
[ 13 -8 -3 ] [ 20 ]
[ -8 10 -1 ] [ -5 ]
Die oberste Gleichung könnte dann durch -3 dividiert werden,
um den Wert 1 in die Plvot-Positeon zu bringen. Das Ergebnis wird
jedoch genauer, wenn wir zuerstst die erste und zweite Zeile vertauschen.
Dies bringt den größeren Wert 13 in die Plvot-Positeon.
Nach Elimination der ersten Variablen aus der zweiten und dritten
Gleichung lautet das Ergebnis:
[ 1 -0.6 -0.2 ] [1.5]
[ 0 -2.8 10.3 ] [4.6]
[ 0 5.1 -2.8 ] [7.3]
Nun sollte man die zweite und dritte Gleichung vertauschen, um das
größere Element 5.1 in die zweite Pivot-Positeon zu bringen.
Es gibt einen weiteren Grund, warum Zeilen vertauscht werden
müssen. Erscheint nämlich ein Null-Element (oder ein Wert nahe Null)
in der Hauptdiagonalen, dann ist es nicht möglich, die Zeile durch dieses
Pivot-Element zu dividieren ohne große Rundungsfehler zu erhalten.
Der Austausch dieser Zeile mit einer darunterliegenden entfernt die Null aus
der Pivot-Positeon.
Das Programm kann mit relativ wenig Aufwand erweitert werden:
#include <stdio.h>
#define N 3 /* Spalten/Zeilen der Matrix */
/* = Anzahl der Gleichungen */
int main(void)
{
int i, j, k; /* Schleifenzaehler */
int merk; /* Merker fuer Zeilentausch */
double f; /* Hilfsvariable */
double a[N][N]; /* Koeffizientenmatrix */
double b[N]; /* rechte Seite */
double x[N]; /* Loesung */
/* Einlesen */
for(i = 1; i <= N; i++)
{
for(j = 1; j <= N; j++)
scanf("%lg", &a[i-1][j-1]);
scanf("%lg", &b[i-1]);
}
/* Elimination */
for(i = 1; i <= (N-1); i++)
{
merk = i;
for(k = i+1; k <= N; k++)
if (abs(a[k-1][i-1]) > abs(a[merk-1][i-1]))
merk = k;
if (merk != i) /* vertauschen */
{
for (j = i+1; j <= N; j++)
{
f = a[i-1][j-1];
a[i-1][j-1] = a[merk-1][j-1];
a[merk-1][j-1] = f;
}
f = b[i-1];
b[i-1] = b[merk-1];
b[merk-1] = f;
}
for(k = i+1; k <= N; k++)
{
f = a[k-1][i-1] / a[i-1][i-1];
for(j = i+1; j <= N; j++)
a[k-1][j-1] -= f*a[i-1][j-1];
b[k-1] -= f*b[i-1];
}
}
/* Ruecksubstitution */
for(i = N; i >= 1; i--)
{
for(j = i+1; j <= N; j++)
b[i-1] -= a[i-1][j-1] * x[j-1];
x[i-1] = b[i-1] / a[i-1][i-1];
}
/* Ausgabe */
for(i = 1; i <= N; i++)
printf("x%d = %g\n", i, x[i-1]);
return 0;
}
Um alle Fährnissen standzuhalten, sollte bei der Suche nach dem größten
Pivot-Element auch gleich noch geprüft werden, ob die Matrix nicht singulär
wird, d. h. ob es betragsmäßig größer als eine vorgegebene
untere Schranke ist.
|
|
|
DIENSTE
