dev cpp - Startseite TU Ilmenau · Hinweise zur Nutzung der Programmieroberﬂ¨ache Dev-C++, der...

TU Ilmenau, Institut fur MathematikFG Numerische Mathematik und InformationsverarbeitungPD Dr. rer. nat. habil. W.Neundorf

Datei: DEV CPP.TEX

Hinweise zur Nutzung der Programmieroberflache Dev-C++,der Formelsammlung zur Numerischen Mathematik in ANSI C undC++ (

”Aachener Bibliothek“) sowie zur Erstellung von Program-

men und Projekten in C++ und C

1 Dev-C++

Version Dev-C++ 4.9.6.0

Die Oberflache mit den GNU-Compilern g++ und gcc lasst ein bequemes Arbeiten mit derEntwicklungsumgebung (Windows IDE) zu. Alle Einstellungen konnen von dort aus erfolgen.Die Arbeit soll an einfachen Beispielen erlautert werden. Nach diesem Muster konnen andereProgramme gleichfalls erstellt und abgearbeitet werden.

1.1 Aufgabe 1

Demoprogramm a1c.cpp im Verzeichnis C:\d\Neundorf\nwptexte\tech phy\05\cpp gpp gcc

// a1c.cpp

// Rechnen mit komplexen Zahlen - Newton-Verfahren fuer Polynome

#include <iostream.h> // deklariert Objekte wie cin, cout,...

#include <stdio.h> // kann entfallen

#include <stdlib.h> // kann entfallen

#include <math.h> // kann entfallen

#include <conio.h> // deklariert Objekte wie getch, printf, scanf, ...

struct komplex

{ double re,im; };

// Addition: z = z1+z2

komplex add(komplex z1,komplex z2)

{ komplex z;

z.re = z1.re+z2.re;

z.im = z1.im+z2.im;

return z; }

// Subtraktion: z = z1-z2

komplex sub(komplex z1,komplex z2)

{ komplex z;

z.re = z1.re-z2.re;

z.im = z1.im-z2.im;

return z; }

// Betrag |z|

double kabs(komplex z)

{ return sqrt(z.re*z.re+z.im*z.im); }

// Multiplikation: z = z1*z2

komplex mult(komplex z1,komplex z2)

{ komplex z;

z.re = z1.re*z2.re-z1.im*z2.im;

z.im = z1.re*z2.im+z1.im*z2.re;

return z; }

// Division: z = z1/z2

komplex div(komplex z1,komplex z2)

{ komplex z,h; double nenner;

h.re = z2.re; h.im = -z2.im;

nenner = z2.re*z2.re+z2.im*z2.im;

z.re = (z1.re*h.re-z1.im*h.im)/nenner;

1

z.im = (z1.re*h.im+z1.im*h.re)/nenner;

return z; }

int main()

{ komplex c1,c2,h,p,ps;

komplex a[11];

char vz;

int n,i;

cout<<"Rechnen mit komplexen Zahlen - Newton-Verfahren\n";

cout<<"pn(x)=a[n]x^n+a[n-1]x^(n-1)+...+a[1]x+a[0] = 0\n";

cout<<"Polynomgrad n (n<=10): "; cin>>n;

cout<<"\nEingabe Polynomkoeffizienten:\n";

for(i=0;i<=n;i++)

{ cout<<"\na["<<i<<"].re: ";

cin>>a[i].re;

cout<<"a["<<i<<"].im: ";

cin>>a[i].im;

}

cout<<"\nStartwert Newton: x[0].re= ";

cin>>c1.re;

cout<<"Startwert Newton: x[0].im= ";

cin>>c1.im;

// Newton-Verfahren, Nullstellenberechnung fuer komplexe Polynome

// unter Nutzung des zweizeiligen Hornerschema

do

{ c2 = c1;

p = a[n]; ps = p;

for(i=n-1;i>0;i--)

{ h = mult(p,c2); p = add(h,a[i]); // p =p*x0+a[i];

h = mult(ps,c2); ps = add(h,p); // ps=ps*x0+p;

}

h = mult(p,c2); p = add(h,a[0]); // p=p*x0+a[0];

c1 = sub(c2,div(p,ps));

}

while(kabs(sub(c1,c2))>1e-10);

vz = (c1.im<0?’-’:’+’);

cout.precision(16);

cout<<"\nNullstelle c1: "<<c1.re<<vz<<fabs(c1.im)<<"*I\n";

getch();

return 0;

}

Die dortigen Include-Anweisungen betreffen Standard-Header-Dateien.Bei der Rechnung zahlreicher Beispiele und Ausgabe vieler Daten sind im Quelltext soge-nannte Haltepunkte einzufugen. Das konnen solche Anweisungen sein wie getch(), dazunoch ev. das Loschen des BS mittels clrscr(). Dabei ist im Demoprogramm noch dieInclude-Anweisung fur die Header-Datei conio.h einzutragen, also

#include <conio.h> /*Direkte MSDOS Console I/O*/

Zur Ausfuhrung sind die folgenden Schritte zu machen.

1. Menupunkt Datei

⇒ Projekt oder Datei offnen... ⇒ Datei offnen

Auswahlfenster: Suchen in:

Dateiname:

Dateityp:

Schalter Offnen

Programm a1c.cpp erscheint in der Anzeigeleiste.

2. Menupunkt Ausfuhren

⇒ Kompilieren

2

Fuehrt g++.exe... aus

g++.exe "C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c.cpp"

-o "C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c.exe"

-s -I"C:\Dev-Cpp\include"

-I"C:\Dev-Cpp\include\g++-3"

-I"C:\Dev-Cpp\include"

-L"C:\Dev-Cpp\lib"

Ausfuehrung beendet

Kompilierung erfolgreich

Im Verzeichnis erscheint die ausfuhrbare Datei a1c.exe.


⇒ Ausfuhren

Es offnet sich das Eingabe- und Ergebnisfenster

2 C:\d\Neundorf\nwptexte\tech phy\05\cpp gpp gcc\a1c.exe

Rechnen mit komplexen Zahlen - Newton-Verfahrenpn(x)=a[n]x n+a[n-1]x (n-1)+...+a[1]x+a[0] = 0Polynomgrad n (n<=10): 2

Eingabe Polynomkoeffizienten:...


⇒ Kompilieren u. Ausfuhren

Dasselbe Programm stellen wir als a1c .c bereit und ubersetzen es. Dazu wird er Compilergcc.exe aufgerufen. Als erstes kommt eine Meldung zur Header-Datei iostream.h, die inC nicht vorgesehen ist, weil es die Objekte cin, cout,... nicht gibt.

4 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c iostream.h: No such file or directory.

Wir kommentieren die entsprechende Zeile heraus und machen den nachsten Versuch, jedochbei zahlreichen weiteren Fehlermeldungen.

14 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c parse error before ‘add’

14 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c parse error before ‘z1’

C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c [Warning] In function ‘add’:

15 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c ‘komplex’ undeclared (first use in this function)

[Build Error] (Each undeclared identifier is reported

[Build Error] only once for each function it appears in.)

15 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c parse error before ‘z’

16 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c ‘z’ undeclared (first use in this function)

16 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c ‘z1’ undeclared (first use in this function)


[Build Error] At top level:

21 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c parse error before ‘sub’


C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c [Warning] In function ‘sub’:








C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c [Warning] In function ‘kabs’:



32 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c parse error before ‘mult’


C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c [Warning] In function ‘mult’:

3







39 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c parse error before ‘div’


40 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c conflicting types for ‘div’

350 C:\Dev-Cpp\include\stdlib.h previous declaration of ‘div’

C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c [Warning] In function ‘div’:



41 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c ‘h’ undeclared (first use in this function)




C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c [Warning] In function ‘main’:


48 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c parse error before ‘c1’

53 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c ‘cout’ undeclared (first use in this function)

55 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c ‘cin’ undeclared (first use in this function)

60 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c ‘a’ undeclared (first use in this function)

66 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c ‘c1’ undeclared (first use in this function)

73 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c ‘c2’ undeclared (first use in this function)

74 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c ‘p’ undeclared (first use in this function)

74 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c ‘ps’ undeclared (first use in this function)

76 C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c_.c ‘h’ undeclared (first use in this function)

Man erkennt, dass eine Reihe von Veranderungen bzw. Korrekturen notwendig ist, um einlauffahiges C-Programm zu erhalten.

(1) In C++ reicht als Typangabe fur eine Strukturvariable der Typbezeichner (hier komplex),in C ist der Typ durch struct komplex zu notieren. Dadurch kommen die meisten Fehler-meldungen.(2) Die Ein-/Ausgabeanweisungen cin, cout sind nicht deklariert und in C durch die Kom-mandos scanf, printf zu ersetzen. scanf erwartet als Argumente Adressen, so dass manden Adressoperator & direkt vor den Variablennamen stellt.(3) div ist in der Bibliothek stdlib.h als Standardfunktion vordefiniert und steht im Kon-flikt mit der gleichnamigen Nutzerfunktion. Diese sollte man anders bezeichnen (hier divi).(4) Die meisten Header-Dateien werden im Ubersetzungsprozess durch Verlinken mit denBibliotheken Dev-Cpp\include und Dev-Cpp\lib schon eingebunden. Damit ist auch dieFunktion getch() verfugbar.

Nach den Korrekturen erhalt man das C-Programm a1c.c.

// a1c.c

// Rechnen mit komplexen Zahlen - Newton-Verfahren fuer Polynome

//#include <stdio.h>

//#include <stdlib.h>

#include <math.h> // muss bleiben

//#include <conio.h>

struct komplex

{ double re,im; };

// Addition: z = z1+z2

struct komplex add(struct komplex z1,struct komplex z2)

{ struct komplex z;

z.re = z1.re+z2.re;

z.im = z1.im+z2.im;

return z; }

// Subtraktion: z = z1-z2

struct komplex sub(struct komplex z1,struct komplex z2)

{ struct komplex z;

z.re = z1.re-z2.re;

z.im = z1.im-z2.im;

return z; }

4

// Betrag |z|

double kabs(struct komplex z)

{ return sqrt(z.re*z.re+z.im*z.im); }

// Multiplikation: z = z1*z2

struct komplex mult(struct komplex z1,struct komplex z2)

{ struct komplex z;

z.re = z1.re*z2.re - z1.im*z2.im;

z.im = z1.re*z2.im + z1.im*z2.re;

return z; }

// Division: z = z1/z2

struct komplex divi(struct komplex z1,struct komplex z2)

{ struct komplex z,h;

double nenner;

h.re=z2.re; h.im=-z2.im;

nenner=z2.re*z2.re+z2.im*z2.im;

z.re = (z1.re*h.re - z1.im*h.im)/nenner;

z.im = (z1.re*h.im + z1.im*h.re)/nenner;

return z; }

int main()

{ struct komplex c1,c2,h,p,ps;

struct komplex a[11];

char vz;

int n,i;

printf("Rechnen mit komplexen Zahlen - Newton-Verfahren\n");

printf("pn(x)=a[n]x^n+a[n-1]x^(n-1)+...+a[1]x+a[0] = 0\n");

printf("Polynomgrad n (n<=10): ");

scanf("%i",&n);

printf("\nEingabe Polynomkoeffizienten:\n");

for (i=0;i<=n;i++)

{ printf("\na[%i].re: ",i);

scanf("%lf",&a[i].re);

printf("a[%i].im: ",i);

scanf("%lf",&a[i].im);

}

printf("\nNewton-Verfahren");

printf("\nStartwert: x[0].re= ");

scanf("%lf",&c1.re);

printf(" x[0].im= ");

scanf("%lf",&c1.im);

// Newton-Verfahren, Nullstellenberechnung fuer komplexe Polynome unter Nutzung des zweizeiligen Hornerschema

do

{ c2 = c1;

p = a[n]; ps = p;

for(i=n-1;i>0;i--)

{ h = mult(p,c2); p = add(h,a[i]); // p =p*x0+a[i];

h = mult(ps,c2); ps = add(h,p); // ps=ps*x0+p;

}

h = mult(p,c2); p = add(h,a[0]); // p=p*x0+a[0];

c1 = sub(c2,divi(p,ps));

}

while(kabs(sub(c1,c2))>1e-10);

vz = (c1.im<0?’-’:’+’);

printf("\nNullstelle: %17.15lf%c%17.15lf*I\n",c1.re,vz,fabs(c1.im));

getch();

return 0;

}

Man ubersetzt es (⇒ a1c.exe) und fuhrt es aus.

Fuehrt gcc.exe... aus

gcc.exe "C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c.c"

-o "C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\a1c.exe"

-s -I"C:\Dev-Cpp\include"

-L"C:\Dev-Cpp\lib"

Ausfuehrung beendet


5

1.2 Aufgabe 2


// a4c.cpp

// Numerische Integration - zusammengesetzte Trapezregel

#include <iostream.h> // deklariert Objekte wie cin, cout, ...

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <conio.h> // deklariert Objekte wie getch, printf, scanf, ...

// Integranden

double f1(double x)

{ return sin(x*x);}

double f2(double x)

{ return sin(x)/x; }

double f3(double x)

{ return sqrt(8.0*x)-x*x/8; }

double f4(double x)

{ return x-x*x*x; }

double f5(double x)

{ return 2.0*exp(-x*x)/sqrt(3.141592653589793); }

// Zusammengesetzte Trapezregel

double Trapez(double f(double),double a,double b,int n)

{

int k;

double I,x,h;

x = a;

h = (b-a)/n;

I = 0.5*(f(a)+f(b));

for(k=1;k<n;k++)

{

x = x+h;

I = I+f(x);

}

return h*I;

}

int main()

{

double a,b;

int n;

a = 0; b = 1;

n = 10;

cout<<"Integralberechnungen mit zusammengesetzter Trapezregel\n";

cout<<"\nIntegrand f5(x)=2 exp(-x*x)/sqrt(Pi)\n";

cout<<"n="<<n<<"\n";

cout.precision(16);

for(b=1;b<4.1;b=b+1.0)

{

cout<<"Integralwert Tn(f5,0,"<<b<<") = "<<Trapez(f5,a,b,n)<<"\n";

}

n = 20;

cout<<"\nn="<<n<<"\n";

for(b=1;b<4.1;b=b+1)

{

cout<<"Integralwert Tn(f5,0,"<<b<<") = "<<Trapez(f5,a,b,n)<<"\n";

}

cout<<"\nExakter Wert = Int(f5,0,unendlich) = 1\n";

getch();

return 0;

}

Alle Parameter fur die Integrationsformel sind im Programm definiert, so dass keine Eingabevorgesehen ist.Nach Ubersetzen und Ausfuhren von a4c.exe offnet sich das Ergebnisfenster

6

2 C:\d\Neundorf\nwptexte\tech phy\05\cpp gpp gcc\a4c.exe

Integralberechnungen mit zusammengesetzter Trapezregel Tn(f,a,b)

Integrand f5(x)=2 exp(-x∗x)/sqrt(Pi)n=10Integralwert Tn(f5,0,1) = 0.842008717Integralwert Tn(f5,0,2) = 0.995048543Integralwert Tn(f5,0,3) = 0.999971913Integralwert Tn(f5,0,4) = 0.999999973

n=20Integralwert Tn(f5,0,1) = 0.842527817Integralwert Tn(f5,0,2) = 0.995253490Integralwert Tn(f5,0,3) = 0.999976360Integralwert Tn(f5,0,4) = 0.999999981

Exakter Wert = Int(f5,0,unendlich) = 1

Wiederrum sind einige Veranderungen bzw. Korrekturen notwendig, um ein analoges lauffahi-ges C-Programm zu erhalten.

// a4c.c

// Numerische Integration - zusammengesetzte Trapezregel



#include <math.h>


// Integranden

double f1(double x)

{ return sin(x*x);}

...

double f5(double x)

{ return 2.0*exp(-x*x)/sqrt(3.141592653589793); }

// Zusammengesetzte Trapezregel

double Trapez(double f(double),double a,double b,int n)

{ int k;

double I,x,h;

...

}

int main()

{ double a,b;

int n;

a = 0; b = 1;

n=10;

printf("Integralberechnungen mit zusammengesetzter Trapezregel Tn(f,a,b)\n");

printf("\nIntegrand f5(x)=2 exp(-x*x)/sqrt(Pi)\n");

printf("n=%d\n",n);

for(b=1;b<4.1;b=b+1.0)

{

printf("Integralwert Tn(f5,0,%1.0lf)=%11.9lf\n",b,Trapez(f5,a,b,n));

}

n = 20;

printf("\nn=%d\n",n);

for(b=1;b<4.1;b=b+1)

{

printf("Integralwert Tn(f5,0,%1.0lf)=%11.9lf\n",b,Trapez(f5,a,b,n));

}

printf("\nExakter Wert = Int(f5,0,unendlich) = 1\n");

getch();

return 0;

}

7

1.3 Aufgabe 3


// a5c.cpp

// Koerpergewicht

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <conio.h> /* deklariert Objekte wie getch, printf, scanf */

int main(void) /* Gibt ’int’ zurueck */

{

short groesse;

float gewicht, ideal;

char auswahl, ret;

char ende; /* Kontrollvariable */

do /* Schleifenanfang */

{

// clrscr(); /* Diese Funktion loescht den Bildschirm! */

// /* Headerdatei: <conio.h> */

/* Achtung: clrscr() funktioniert nicht mit Visual C++, Dev-C++,... */

printf("\n\n\t Bitte geben Sie Ihre Koerpergroesse (in cm) ein: ");

scanf("%d%c", &groesse, &ret);

printf("\n\t Bitte geben Sie Ihr Gewicht (in kg) ein: ");

scanf("%f%c", &gewicht, &ret);

printf("\n\t Bitte geben Sie Ihr Geschlecht an:");

printf("\n\t m = Mann");

printf("\n\t f = Frau\n\t "); // !=m ist Frau

auswahl = getch();

printf("%c\n", auswahl);

/* Hier erfolgt nun der Aufruf unserer neuen Funktion: */

ideal = errechne_ideal(groesse, auswahl);

printf("\n\t Ihr Idealgewicht waere: %5.1f kg", ideal);

if(ideal<gewicht)

printf("\n\n\t Sie haben %5.1f kg Uebergewicht!", gewicht - ideal);

else

printf("\n\n\t Sie haben %5.1f kg Untergewicht!", ideal - gewicht);

printf("\n\n\t Moechten Sie das Programm beenden (j/n)?");

ende = getch();

}while(ende!=’j’); /* Abbruchbedingung */

/* Schleifenende */

printf("\n\n\t Taste druecken...");

getch();

return 0; /* Ein Programm das 0 zurueck gibt */

/* wurde ordnungsgemaess beendet */

}

// Definition des Prototyps

float errechne_ideal(int groesse, char geschlecht)

{

float ideal;

ideal = groesse-100;

if(geschlecht==’m’)

ideal = ideal/100*90;

else

ideal = ideal/100*85;

return ideal; /* gibt den Inhalt von ’ideal’ zurueck */

}

Alles scheint in Ordnung zu sein, aber bei dem Versuch das Programm zu kompilieren weigertsich g++ den Quelltext zu ubersetzen.

[Warning] In function ‘int main()’:

implicit declaration of function ‘int errechne_ideal(...)’

8

Das liegt daran, dass die Funktion errechne_ideal() erst nach der Hauptfunktion definiertwurde. So mancher C-Compiler lasst sich uberlisten, indem man die Funktion einfach vorder Hauptfunktion definiert. Die korrekte Losung besteht jedoch in der Verwendung einesPrototypes, an dem der Compiler erkennen kann das irgendwo weiter unten im Quelltextnoch eine entsprechende Funktion definiert ist.

// Deklaration des Prototyps

float errechne_ideal(int groesse, char geschlecht);

Diesen fugt man als Zeile zwischen die #include-Anweisungen und die main-Funktion ein.Dies ist der ubliche Ort an dem Prototypen platziert werden sollten.

Da die Eingabe von Daten mit der <Enter>-Taste beendet wird, haben wir das <Enter>-Zeichen extra einer Variablen zugeordnet.

scanf("%d%c", &groesse, &ret);

Damit sichert man, dass der Tastaturpuffer geleert wird. Man kann dies auch mit demAufruf der Methoden flush() bzw. fflush(stdin) mit dem Standardbezeichner stdin ausden Header-Files stdio.h oder conio.h tun. Oft passiert dies jedoch automatisch mit denAbschluss der Eingabe einfach durch die <Enter>-Taste und nach der Anweisung cout<<.

Es sind keine Veranderungen bzw. Korrekturen notwendig, um daraus ein lauffahiges C-Programm zu erhalten. Man bezeichnet einfach die Datei um in a5c.c.

1.4 Aufgabe 4

Demoprogramm matsum.cpp im VerzeichnisC:\d\Neundorf\nwptexte\tech phy\05\cpp gpp gcc

// matsum.cpp

// Summe c = a+b zweier Matrizen der Dimension (m,n)

#include <iostream.h> // deklariert Objekte wie cin, cout,...

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <conio.h> // deklariert Objekte wie getch, printf, scanf,...

void summe(int m,int n,long double x[][3],long double y[][3],long double z[][3]); // long double = extended

int main()

{

long double a[2][3] = {0.987654321e-12,-2,1,3,4,7};

long double b[2][3] = {2.0e-5,4,8,-5,-5,-5};

// man teste und beobachte die Ausgabe

// long double b[2][3] = {2.0e-4,4,8,-5,-5,-5};

// long double b[2][3] = {2.0e-3,4,8,-5,-5,-5};

// ...

long double c[2][3] = {1.0,1,1,1,1,1};

cout << "\nSumme c = a+b von Matrizen" ;

cout << "\n============================" << endl;

summe(2,3,a,b,c);

cout.precision(16); // double !

cout << "c[0][0] = " << c[0][0] << endl;

cout << "c[0][1] = " << c[0][1] << endl;

cout << "c[0][2] = " << c[0][2] << endl;

cout << "c[1][0] = " << c[1][0] << endl;

cout << "c[1][1] = " << c[1][1] << endl;

cout << "c[1][2] = " << c[1][2] << endl;

getch();

return 0;

}

9

void summe(int m, // Zeilenzahl

int n, // Spaltenzahl

long double x[][3], // erste Matrix

long double y[][3], // zweite Matrix

long double z[][3]) // Summenmatrix

{

for(int i=0;i<m;++i)

for(int k=0;k<n;++k)

z[i][k] = x[i][k]+y[i][k];

}

Ausgewahlte Ergebnisse

Summe c = a+b von Matrizen

============================

c[0][0] = 2.000000098765432e-05

c[0][1] = 2

c[0][2] = 9

c[1][0] = -2

c[1][1] = -1

c[1][2] = 2

---------------------------------------------------------------------------------

cout.precision(16);

long double a[2][3] = {0.987654321e-12,-2,1,3,4,7};

long double b[2][3] = {2.0e-6,4,8,-5,-5,-5}; => c[0][0] = 2.000000987654321e-06

long double b[2][3] = {-2.0e-5,4,8,-5,-5,-5}; => c[0][0] = -1.999999901234568e-05

long double b[2][3] = {2.0e-5,4,8,-5,-5,-5}; => c[0][0] = 2.000000098765432e-05

long double b[2][3] = {2.0e-4,4,8,-5,-5,-5}; => c[0][0] = 0.0002000000009876543

long double b[2][3] = {2.0e-3,4,8,-5,-5,-5}; => c[0][0] = 0.002000000000987654

long double b[2][3] = {2.0e-2,4,8,-5,-5,-5}; => c[0][0] = 0.02000000000098765

long double b[2][3] = {2.0e-1,4,8,-5,-5,-5}; => c[0][0] = 0.2000000000009877

long double b[2][3] = {2.0e0,4,8,-5,-5,-5}; => c[0][0] = 2.000000000000988

long double b[2][3] = {2.0e1,4,8,-5,-5,-5}; => c[0][0] = 20.00000000000099

Man bemerkt den Wechsel im Ausgabeformat bei 16 gultigen dezimalen Mantissenstellen.Sobald die auszugebene Dezimalzahl betragsmaßig kleiner als 10−4 ist, erfolgt ein Ubergangzum Gleitpunktformat mit skalierter Mantisse.

Um ein analoges lauffahiges C-Programm zu erhalten, sind einige Veranderungen bzw. Kor-rekturen notwendig.

1. Der Datentyp long double steht in C nicht zur Verfugung. Es wird zwar, wenn mandamit arbeitet, kein syntaktischer Fehler angezeigt, aber bei Uberprufung der Initialisierungder Felder erkennt man ihre falsche Belegung.2. In C kann die Schleifenvariable nicht direkt im Initialisierungsteil deklariert werden, diesmuss im umliegenden Bereich erfolgen.

// matsum.c

// Summe c = a+b zweier Matrizen der Dimension (m,n)



#include <math.h>


10

void summe(int m,int n, double x[][3], double y[][3],

double z[][3]); // nicht in C: long double = extended

int main()

{

double a[2][3] = {0.987654321e-12,-2,1,3,4,7};

double b[2][3] = {2.0e-5,4,8,-5,-5,-5};

double c[2][3] = {1.0,1,1,1,1,1};

// auch moeglich

// double c[2][3] = {{1.0,1,1},

// {1,1,1}};

printf("\nSumme c = a+b von Matrizen");

printf("\n============================\n");

summe(2,3,a,b,c);

printf("c[0][0] = %17.15e\n",c[0][0]);

printf("c[0][1] = %.0lf\n",c[0][1]);

printf("c[0][2] = %.0lf\n",c[0][2]);

printf("c[1][0] = %.0lf\n",c[1][0]);

printf("c[1][1] = %.0lf\n",c[1][1]);

printf("c[1][2] = %.0lf\n",c[1][2]);

getch();

return 0;

}

void summe(int m, // Zeilenzahl

int n, // Spaltenzahl

double x[][3], // erste Matrix

double y[][3], // zweite Matrix

double z[][3]) // Summenmatrix

{ int i,k;

for(i=0;i<m;++i) // Deklaration von i,k nicht im Initialisierungsteil

for(k=0;k<n;++k)

z[i][k] = x[i][k]+y[i][k];

}

Ergebnisse

Summe c = a+b von Matrizen

============================

c[0][0] = 2.000000098765432e-005

c[0][1] = 2

c[0][2] = 9

c[1][0] = -2

c[1][1] = -1

c[1][2] = 2

11

1.5 Aufgabe 5

Demoprogramm tnewton1.c im Verzeichnis C:\d\Neundorf\nwptexte\tech phy\05\dev cpp

Die Definitionen und Prototypen, also die Prozedurkopfe usw., findet man in den Header-Files basis.h, u_proto.h, tfunc1.h als Vorabdeklarationen und werden damit bekanntgemacht.Ihre eigentliche Definition erfolgt in den Funktionen fnewtonc.c, tfunc1.c, basis.c,die nach der Hauptfunktion main stehen.

// tnewton1.c

/*-------------------------------------------------------------------*/

/* Test program for newton */

/*-------------------------------------------------------------------*/

//#include <math.h> ... in Header-File basis.h

#include <conio.h> // wegen getch()

#include "basis.h"

#include "u_proto.h"

#include "tfunc1.h"

int main (void)

{

int i, rc, iter;

REAL x, f;

REAL (*fct) (REAL), (*fctd) (REAL);

char *text;

WriteHead ("Newton Method for real, nonlinear functions");

for(i=1;i<=6;i++)

{

switch (i)

{

case 1: x = ONE;

fct = f1;

fctd = f1s;

text = "f = 5*x-exp(x)";

break;

case 2: x = TWO;

fct = f2;

fctd = f2s;

text = "f = (((((x-6)*x+15)*x-20)*x+15)*x-6)*x+1";

break;

case 3: x = ONE;

fct = f3;

fctd = f3s;

text = "f = sin(x)";

break;

case 4: x = ONE;

fct = f4;

fctd = f4s;

text = "f = 1+sin(x)";

break;

case 5: x = TWO;

fct = f5;

fctd = f5s;

text = "f = exp(x)-(1.0+x+x*x*0.5)";

break;

case 6: x = TWO;

fct = f6;

fctd = f6s;

text = "f = (x-1.0)*(x-1.0)*(sin(PI*x)-log(fabs(2.0*x/(x+1.0)))";

break;

default: return (1);

}

printf("%s\n", text);

printf("Start value x0 = ");

printf(FORMAT_LF, x);

printf("\n");

rc = newton(fct,fctd,&x,&f,&iter);

12

printf("Return code = % d\n", rc);

printf("Root = "); printf (FORMAT_2016LF, x);

printf("\nFunction value = "); printf (FORMAT_LE, f);

printf("\nIterations = % d\n\n", iter);

}

WriteEnd();

getch();

return (0);

}

// fnewton.c

/* Loesung nichtlinearer Gleichungen

Numerische Verfahren zur Loesung nichtlinearer Gleichungen

Newton-Verfahren fueur einfache und mehrfache Nullstellen */

#include "basis.h"


#define ITERMAX 300 /* Maximale Iterationszahl */

#define ABSERR ZERO /* Zugelassener Absolutfehler */

#define RELERR (REAL)((REAL)128.0 * MACH_EPS) /* Zugelassener Relativfehler */

#define FCTERR (REAL)((REAL)4.0 * MACH_EPS) /* Max. Fehler im Funktionwert*/

int newton /* Eindimensionales Newton Verfahren .........*/

( REALFCT fct, /* Funktion ........................*/

REALFCT fderv, /* 1. Ableitung ....................*/

REAL * x, /* Startwert / Loesung .............*/

REAL * fval, /* Funktionswert an Loesung........ */

int * iter /* Iterationszahl ..................*/

)

/* Die Funktion newton realisiert das Newton-Iterationsverfahren *

* zur Loesung der Gleichung fct(x) = 0. *

* Die Funktion fct und deren 1. Ableitung muessen als Parameter *

* uebergeben werden. */

...

// tfunc1.c

/* Definition of test functions for newton, pegasus and roots */

#include "basis.h"


#include "tfunc1.h"

REAL f1 (REAL x) /* f(x) = 5*x -exp(x) */

{ return ((REAL)5.0*x-EXP(x)); }

REAL f1s (REAL x)

{ return ((REAL)5.0-EXP(x)); }

REAL f1ss (REAL x)

{ return (-EXP(x)); }

...

// basis.c

/* grundlegende Funktionen: Definitionsdatei */

#include "basis.h" /* wegen NULL, freopen, stdout, fprintf, stderr, stdin, */

/* SQRT, EXP, sqrt, MACH_EPS, POSMAX, epsquad, maxroot, */

/* pi, ATAN, sqr, umleiten, readln, intervall, horner, */

/* norm_max, skalprod, copy_vector, REAL, ONE, TWO, FOUR, */

/* ZERO, HALF, FABS, boolean, FOUR, basis, mach_eps, */

/* epsroot, exp_1, posmin, sqrtlong, comdiv, comabs, */

/* quadsolv, SetVec, CopyVec, ReadVec, WriteVec, SetMat, */

/* CopyMat, ReadMat, WriteMat, WriteHead, WriteEnd, */

/* LogError, fgetc, stdin, SWAP */

...

Die Programme sind aus der Bibliothek FORMELN\CNUM8 (Formelsammlung zur Nu-merischen Mathematik in ANSI C und C++, Aachener Bibliothek) entnommen und zusam-mengestellt.Das Programm tnewton1.c wir mit dem Compiler gcc ubersetzt (→ tnewton1.exe) undman fuhrt es aus.

13

Das Ergebnisfenster ist

2 C:\d\Neundorf\nwptexte\tech phy\05\dev cpp\tnewton1.exe——————————————————————————–Newton Method for real, nonlinear functions——————————————————————————–

f = 5∗x-exp(x)Start value x0 = 1.000000Return code = 0Root = 0.2591711018190737Function value = -1.758576e-016Iterations = 5

f = (((((x-6)∗x+15)∗x-20)∗x+15)∗x-6)∗x+1Start value x0 = 2.000000Return code = 2Root = 1.0059508814409825Function value = 4.451994e-014Iterations = 300

f = sin(x)Start value x0 = 1.000000Return code = 0Root = 0.0000000000000000Function value = 0.000000e+000Iterations = 5

f = 1+sin(x)Start value x0 = 1.000000Return code = 0Root = -1.5707963497819439Function value = 2.642201e-016Iterations = 26

f = exp(x)-(1.0+x+x∗x∗0.5)Start value x0 = 2.000000Return code = 0Root = 0.0000148767867855Function value = 4.037569e-016Iterations = 162

f = (x-1.0)∗(x-1.0)∗(sin(PI∗x)-log(fabs(2.0∗x/(x+1.0)))Start value x0 = 2.000000Return code = 0Root = 2.0981196117437797Function value = -2.541919e-016Iterations = 5

——————————————————————————–

14

2 Aufbau der Numerik-Bibliothek

In Abschnitt 1, Aufgabe 5, ist eine Programmbibliothek verwendet worden.Die Aachener Bibliothek ist hierarchisch aufgebaut, wohl strukturiert und nach Teildiszipli-nen der Numerischen Mathematik angelegt.

1. Die Programme befinden sich auf dem Novell-Netz RUNGE (Volume share) des IfMathQ:\Neundorf\stud m93\FORMELN\ANSICNUM, ...\CNUM8

Die neuere Version CNUM8 om 26.04.1996 enthalt mehr und verbesserte numerischeAlgorithmen.Lokal ist das Bibliotheksverzeichnis z. B. unter

c:\d\Neundorf\FORMELN\ANSICNUM, ...\CNUM8

2. Literaturhinweis- Engeln-Mullges, G., Reutter, F.: Numerik-Algorithmen mit ANSI-C-Programmen.

Wissenschaftsverlag Mannheim 1993. Standort: 50/51=MAT 93 A 7745.- Engeln-Mullges, G., Reutter, F.: Formelsammlung zur Numerischen Mathematik

mit FORTRAN 77-Programmen. Bibliogr. Institut Mannheim 1988 (auch fur TP, C).

3. Informationen zum Aufbau der Bibliothek, ihrer Verzeichnisstruktur und zu ihrer Be-nutzung (z. B. Anwendung und Einstellungen fur unterschiedliche Compiler, Benut-zungsarten) findet man in der Datei LIESMICH.TXT im jeweiligen Verzeichnis.

4. Analog gibt es Verzeichnisse von Versionen der Numerik-Tools fur andere Program-miersprachen, z. B. fur Turbo Pascal im Verzeichnis \TPNUM mit ahnlicher Struktur.

5. Inhalte der einzelnen Unterverzeichnisse

Verzeichnis Inhalt

\CNUM8\02 Quellen der Verfahren fur nichtlineare Gleichungen\CNUM8\02\TST Quellen der zugehorigen Testprogramme

\CNUM8\03 Quellen der Verfahren fur Polynomnullstellen\CNUM8\03\TST Quellen der zugehorigen Testprogramme\CNUM8\03\EIN Polynomkoeffizienten

\CNUM8\04 Quellen der Verfahren fur LGS, EWP\CNUM8\04\TST Quellen der zugehorigen Testprogramme\CNUM8\04\EIN Matrizen, Vektoren, EW, EV

\CNUM8\05.................. Quellen weiterer Verfahren\CNUM8\17 Quellen der Verfahren fur AWP von DGL und DGL-Syst.\CNUM8\18 Quellen der Verfahren fur RWP von DGL\CNUM8\18\TST\CNUM8\18\EIN

\CNUM8\BASIS Quellen wichtiger Basisprogramme\CNUM8\INC enthalt samtliche Header-Dateien\CNUM8\LIB Library\CNUM8\REST i.Allg. nicht benotigt\CNUM8\MAKEFILE.MK\CNUM8\LIESMICH.TXT

15

3 Projekte in Borland-C

3.1 Borland-C++3.1 (DOS)

Der Borland-C-Compiler lasst ein bequemes Arbeiten mit der Entwicklungsumgebung (DOSIDE) zu. Alle Einstellungen konnen von dort aus erfolgen.

1. Beispiel: Nullstellenbestimmung mit dem Newton-Verfahren

Die mogliche Arbeit mit dem Compiler C++3.1 sol an einem ersten einfachen Beispielerlautert werden. Nach diesem Muster konnen andere Programme gleichfalls erstelltund abgearbeitet werden.

Aufgabe: das Demoprogramm TNEWTON.C im Verzeichnis \CNUM8\02\TST istzu ubersetzen und auszufuhren. Dazu sind die folgenden Schritte empfehlenswert.

(a) Man lege im Verzeichnis \CNUM8 ein Arbeitsverzeichnis fur EXE- und OBJ-Dateien an, im weiteren mit \EXE bezeichnet.

(b) Man starte den Borland-Compiler, und wechsle in das Verzeichnis \CNUM8\02\TST.

(c) Man lade das Demoprogramm TNEWTON.C. Die dortigen Include-Anweisungenfur die Header-Dateien enthalten schon wichtige Hinweise dafur, welche Files zumProjekt gehoren.

#include <basis.h> /*Grundlegender Deklarationsteil*/#include <u_proto.h> /*Vordeklaration aller Bibliotheksfunktionen*/#include <tfunc1.h> /*Testfunktionen*/

(d) Offnen eines Turbo C - Projektfiles TNEWTON.PRJ mittelsProject, Open project.

(e) Man fuge mittels Project, Add item die folgenden benotigten Files in das Pro-jekt ein.

\CNUM8\02\TST\TNEWTON.C\CNUM8\02\FNEWTON.C\CNUM8\02\TST\TFUNC1.C\CNUM8\BASIS\BASIS.C

Noch einmal der Hinweis: Erkenne ev. die benotigten Files aus den Include-Anweisungen der Anwendungsbeispiele.

(f) Man stelle mittels Options, Directories die zusatzlichen (!) Verzeichnisse ein.

Include Directories : ;..\CNUM8\INCLibrary : keine Anderung oder ev. ;..\CNUM8\LIBOutput : ..\CNUM8\EXESource : keine Anderung oder ev. ..\CNUM8\02\TST

Das standardmaßige Include-Verzeichnis ist meist ..\BC31\INCLUDE.Die Standardbibliotheken sind ..\BC31\LIB, ..\BC31\CLASSLIB\LIB,..\BC31\OWL\LIB, ..\BC31\LIB\STARTUP.

(g) Man ubersetzt mittels Compiler, Build all das Demoprogramm TNEWTON.C,d.h. das Projekt TNEWTON.PRJ. Der Ablauf erfolgt in den Schritten Compiling

und Linking. Die ausfuhrbare Datei TNEWTON.EXE steht im Verzeichnis\CNUM8\EXE zur Verfugung.

16

(h) Bei erfolgreicher Ubersetzung ist eine Kontrolle der in TNEWTON.C einbezo-genen Header-Dateien moglich mittels Project, Include files.

(i) Start des ubersetzten Projekts gleich in der Entwicklungsumgebung mittels Run,

Run.Die letzten Ergebnisse sind mittels Window, Output zu besichtigen. Bei derRechnung zahlreicher Beispiele und Ausgabe vieler Daten uber mehrere Bildschir-me sind im Quelltext sogenannte Haltepunkte einzufugen. Das konnen solche An-weisungen sein wie getch(), dazu noch ev. das Loschen des BS mittels clrscr().Dabei ist moglicherweise im Demoprogramm noch die Include-Anweisung fur dieHeader-Datei conio.h einzutragen, also

#include <conio.h> /*Direkte MSDOS Console I/O*/

(j) Es empfiehlt sich, vor Verlassen der Entwicklungsumgebung das Projekt zu spei-chern. Mittels Options, Save, Project wird die Projektdatei TNEWTON.PRJin das Verzeichnis \CNUM8\02\TST gespeichert. Auch die gesetzten Directorieswerden gespeichert. Damit ist eine erneute Behandlung dieses Projekts schnellmoglich.

(k) Naturlich kann man die Datei TNEWTON.EXE aus der DOS- oder Windows-Oberflache starten.

2. Beispiel: Losung einer Anfangswertaufgabe fur Systemevon gewohnlichen Differentialgleichungen mittels Runge-Kutta-und anderer Verfahren mit Schrittweitensteuerung

Aufgabe: das Demoprogramm M AWP.C im Verzeichnis \CNUM8\17\TST ist zuubersetzen und auszufuhren. Dazu sind ahnlich wie im obigen Beispiel die folgendenSchritte auszufuhren. Wir beschranken uns auf die Angabe einiger Unterschiede.

(a) Das Demoprogramm ist M AWP.C.

(b) Offnen eines Projektfiles M AWP.PRJ mittels Project, Open project.

(c) Man fuge mit Project, Add item die folgenden benotigten Files in das Projektein.

\CNUM8\17\TST\M AWP.C\CNUM8\17\AWP.C\CNUM8\17\TST\T DGLS.C\CNUM8\BASIS\VMBLOCK.C\CNUM8\BASIS\BASIS.C

(d) Es gibt noch ein Verzeichnis \CNUM8\17\EIN mit Eingabedateien der Parameterder verschiedenen Beispiele.

(e) Normalerweise erfolgt eine Eingabe im Dialog. Gemaß der im Quelltext T DGLS.Cvorliegenden Beispiele sind die Beispielnummer und danach die Verfahrenspara-meter nrbeisp, εa, εr, nrmeth, x0, y1(x0), y2(x0), h, xend, fmax einzugeben.

Die Eingabeaufforderung kann wahlweise angezeigt oder unterdruckt werden. Umdiese anzuzeigen, ist im Demoprogramm die Anweisung

#define INTERAKTIV

noch zu erganzen.Haltepunkte wie getch() sind bei der Ausgabe großer Datenmengen als auchnach Fehlermeldungen sinvoll.

17

Ubersetzen mit dem Borland-C-Compiler C++3.1 (Windows)

Unter Windows sollte man ein eigenstandiges Projekt erstellen.Ein

”Zusammenspiel“ mit DOS ist nicht ganz problemlos.

Dazu die Bemerkungen:

1. Es konnen Projekte, die unter DOS IDE erstellt wurden, auch in die Windows IDE laden.Dann erscheint ein Dialogfenster mit der Frage, ob die Optionen des DOS IDE umzukonfi-gurieren sind.

2. EXE-Dateien von Projekten unter Windows IDE nur unter Windows aufrufen.3. Verwendet man ein Projekt aus Windows IDE in der Entwicklungsumgebung DOS, ist

aufgrund von Voreinstellungen durch Windows die Ubersetzung Compiler, Build all

moglich, nicht aber die Ausfuhrung Run, Run.

3.2 Borland-C++4.52, C++5.02

Die Compiler lassen ein bequemes Arbeiten mit der Entwicklungsumgebung (IDE) zu. Manbemerke jedoch eine Reihe von Unterschieden zur Oberflache von C++3.1.Alle Einstellungen konnen von der IDE aus erfolgen, auch die Wahl von Windows-Oberflacheoder DOS-Standard.

1. Erstellen eines neuen Projekts.Die einzelnen Schritte seien grob umrissen.

(a) File, New, Project, New Target (C++5.02)Projekt, Neues Projekt, Neues Ziel (C++4.52)Projektverzeichnis und -name, Zielname eintragen. Die Typkennung ist ide.Zur Illustration sei der gewahlte Projektname ..\CNUM8\17\TST\proj0001.ide

Weiterhin ist der Zieltyp anzugeben. Unter den zahlreichen Varianten sind hierzwei gebrauchliche genannt:

- Anwendung[.exe] mit Umgebung DOS (Standard)- EasyWin[.exe] mit Umgebung Windows 3.x(16),

Ausgabe in einem Windows-Fenster

Es wird nun das Fenster Projekt geoffnet mit den Standardvorgaben· - proj0001[.exe] (LinkTarget)· |− 2 proj0001[.cpp] (CppCompile)|· |− 2 proj0001[.def] (SourceInclude)|· −2 proj0001[.rc] (CompilerResources)

Dort kann nun das Projekt konfiguriert werden, indem insbesondere uberflussigeKnoten (Programme) geloscht und neue Knoten hinzugefugt werden. Zwecks Zu-sammenstellung der Projektfiles (hierarchische Struktur) gibt es die- Arbeit mit der rechten Maustaste- oder Tastenbetatigung (Einf, Entf, ↑, ↓, ... ).Ein zusammengestelltes Beispielprojekt ist im Punkt 3 angegeben.

Achtung!

Beim Zieltyp Anwendung[.exe] in der DOS-Umgebung (Plattform=DOS(Standard))sind die Einstellungen in den Standardbibliotheken zu kontrollieren. Insbesonde-re ist bei Nutzung von Graphik in C-Programmen explizit die Bibliothek BGIhinzuzufugen.

18

Vorgehensweise:

- Im Projektfenster mittels rechter Maustaste auf proj0001[.exe] klicken.- TargetExpert auswahlen.- Im TargetExpert-Fenster die Standardbibliothek BGI erganzen.

(b) Optionen, Projekt

Im Fenster Projektoptionen sind Verzeichnisse zu aktualisieren.

- QuelltextverzeichnisseInclude : ;..\CNUM8\INC

Bibliothek : keine Anderung oder ev. ;..\CNUM8\LIBQuelltext : ..\CNUM8\17\TST

- AusgabeverzeichnisseOBJ-Dateien : ..\CNUM8\EXE

EXE-Datei : ..\CNUM8\EXE

(c) Optionen, Speichern

Im Fenster Optionen speichern der Konfiguration des Projekts.

(d) Projekt, Projekt neu compilieren

Hier erfolgt die Ubersetzung des Gesamtprojekts.Im Fenster Compilerstatus wird der Ablauf dieses Prozesses angezeigt.Im Fenster Meldung erscheinen Informationen zur Ubersetzung sowie zu ev.Warnungen und Fehlern.

(e) Debug, Ausfuhren (oder Strg F9 bzw.”Blitz“-Ikone)

Starten des Projekts. Der Eingabedialog sowie die Ergebnisdarstellung erfolgenin einem Windows-Fenster. Vor einem wiederholten Start des Projekts ist diesesi. Allg. neu zu offnen (außer bei Fehlermeldugen).

(f) Nach Beendigung der Arbeit am Projekt konnen einige Dateien geloscht werden,so z. B. im Verzeichnis \CNUM8\17\TST die Dateien- proj0001.dsw (Turbo C Context File)- proj0001.csm (Konfigurationsfile)- proj0001.˜de (Bak-Datei von proj0001.ide)- proj0001.obr (Browser State File)

Achtung!

Auch ein Programm, das nur Standard-Header-Dateien enthalt, sollte als Projekt de-finiert werden (nicht Bedingung), um damit den Zugriff und die Einbindung dieser injedem Fall zu sichern.

Beispiel:

C++ Programm matbau6a.cpp enthalt die Header-Files uber die Include-Anweisungen

#include <stdio.h>#include <iostream.h>#include <malloc.h>#include <stdlib.h>#include <math.h>

Das Projekt hat die kurze Darstellung

· - proj0002[.exe]||· −2 matbau6a[.cpp]

19

2. Arbeit mit vorhandenem Projekt.

Als Beispiel diene das Programm M AWP.C im Verzeichnis ..\CNUM8\17\TST.

(a) Datei, Offnen

von ..\CNUM8\17\TST\M AWP.C

(b) Projekt, Projekt offnen

von ..\CNUM8\17\TST\PROJ0001.IDEIm Projektfenster ist folgender Baum mit den Files des Projekts.

· - proj0001[.exe]· |− 2 m awp[.c]· |− 2 t dgls[.c]· |− 2 ..\awp[.c]· |− 2 ..\..\basis\basis[.c]|· −2 ..\..\basis\vmblock[.c]

(c) Debug, Ausfuhren (oder Strg F9 bzw.”Blitz“-Ikone)

Starten des Projekts

3. Zusammenfassung mehrerer unabhangiger Programme in einem Projekt.

Als Beispiel dienen die Programme matbau6a.cpp ... matbau6e.cpp.Dazu kann man folgenden Baum im Projektfenster erzeugen.

· - proj0003[.exe]· |− 2 matbau6a[.cpp]· |− 2 matbau6b[.cpp]· |− 2 matbau6c[.cpp]· |− 2 matbau6d[.cpp]|· −2 matbau6e[.cpp]

Da aber immer nur eines dieser Programme aktiv sein kann und ausgefuhrt werdenkann, mussen die anderen mittels- Anglicken der rechten Maustaste,- Lokale Optionen bearbeiten ... (Edit local options)Optionen, Attribute

- Vom Elternknoten ausschließen (Build Attributes, Exclude from parents)vom Elternknoten (hier matbau6a.cpp) ausgeschlossen werden. Anstelle der Markie-rung “·“ erhalten die ausgegrenzten Knoten die neue Markierung “ “.

Die Darstellung im Projektfenster ist dann

· - proj0003[.exe]· |− 2 matbau6a[.cpp]|− 2 matbau6b[.cpp]|− 2 matbau6c[.cpp]|− 2 matbau6d[.cpp]|−2 matbau6e[.cpp]

Naturlich kann auch jedes andere Programm so aktiviert werden.

20

4 Projekte in Dev-C++

Wir demonstrieren die Erstellung von Projekten und den Umgang damit anhand von einigenBeispielen. Dabei wollen wir schrittweise vorgehen unter Verwendung der Informationen ausden vorherigen Abschnitten.

4.1 Aufgabe 1

In den Demoprogrammen wird der Gauß-Algorithmus zur Losung des LGS Ax = b in ver-schiedenen Versionen implementiert. Die Programme befinden sich im Verzeichnis

C:\d\Neundorf\nwptexte\tech phy\05\dev cpp\projekt1.

Variante 1

Demoprogramme gauss1.c und gauss1.cpp (außer den Dateinamen keine Unterschiede)

/* gauss1.c bzw. gauss1.cpp */

/* Resttableau-Gauss-Algorithmus mit Spaltenpivotisierung und Zeilenvertauschung */

/* Tausch mit ersten NNE in der Spalte */

/* Ax=b, A=>oberes U, unteres Dreieck gleich Null */

/* Ux=c Loesung mit Rueckwaertseinsetzen */

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <conio.h>

#define Nmax 10

// Vowaertsdeklarationen, Prototypen

int Gauss();

void Einlesen();

void Rausschreiben();

int SucheNichtNull(int);

void NulleSpalte(int);

void BerechneX();

void Vertausche(int, int);

void Zeile_abziehen(int, int);

int N;

double A[Nmax][Nmax],b[Nmax],x[Nmax]; // statische Felder

int main()

{ int flag;

Einlesen();

flag = Gauss();

if(flag!=0)

{

printf("Gleichungssystem nicht loesbar!\n");

getch();

return(-1);

}

Rausschreiben();

return(0);

}

/*---------------------------------------------------------------------------*/

int Gauss()

{ int i,j,flag;

for(i=0;i<N;i++)

{

flag = SucheNichtNull(i);

if(flag!=0) return(-1);

NulleSpalte(i);

}

BerechneX();

return(0);

}

21

void Einlesen() // Dateiarbeit

{ FILE *datin;

int i,j;

datin = fopen("gauss1-input.dat","r");

fscanf(datin,"%i",&N);

for(i=0;i<N;i++)

for(j=0;j<N;j++) fscanf(datin,"%lf",&A[i][j]);

for(i=0;i<N;i++) fscanf(datin,"%lf",&b[i]);

fclose(datin);

}

/*---------------------------------------------------------------------------*/

void Rausschreiben() // Dateiarbeit

{ FILE *datout;

int i,j;

datout = fopen("gauss1-output.dat","w");

for(i=0;i<N;i++) fprintf(datout,"%lf\n",x[i]);

fclose(datout);

}

/*---------------------------------------------------------------------------*/

int SucheNichtNull(int i)

{ int k;

for(k=i;k<N;k++)

if(A[k][i]!=0.0)

{

Vertausche(i,k);

return(0);

}

return(-1); // kein NNE in der i-ten Spalte, A singulaer

}

/*---------------------------------------------------------------------------*/

void NulleSpalte(int i)

{ int k;

for(k=i+1;k<N;k++) Zeile_abziehen(i,k); // i-tes Resttableau modifizieren

}

/*---------------------------------------------------------------------------*/

void BerechneX(void) // Rueckwaertseinsetzen fuer Loesung

{ int k,l;

for(k=N-1;k>=0;k--)

{

x[k] = b[k];

for(l=k+1;l<N;l++) x[k] = x[k]-A[k][l]*x[l];

x[k] = x[k]/A[k][k];

}

}

/*---------------------------------------------------------------------------*/

void Vertausche(int i, int j)

{ double temp;

int k;

for(k=i;k<N;k++)

{

temp = A[i][k];

A[i][k] = A[j][k];

A[j][k] = temp;

}

temp = b[i];

b[i] = b[j];

b[j] = temp;

}

/*---------------------------------------------------------------------------*/

void Zeile_abziehen(int i, int j)

{ double quotient;

int k;

quotient = A[j][i]/A[i][i];

for(k=i;k<N;k++) // Resttableau modifizieren

A[j][k] = A[j][k]-quotient*A[i][k]; // dabei i-te Spalte zu Null machen

b[j] = b[j]-quotient*b[i];

}

22

Eingabedatei: gauss1-input.dat

2

0 1

2 0

1

4

Ausgabedatei: gauss1-output.dat

2.000000

1.000000

Variante 2

Demoprogramme gauss21.cpp, gauss22.cpp und gauss21.c

Anders als bei Variante 1 wollen wir den Gauß-Algorithmus in Form des Resttableau-Algorithmus als kompakte Funktion und die Ein-/Ausgabe als selbstandige Funktionen de-finieren und dabei die Moglichkeit der Eingabe- und Ruckgabeparameter nutzen. Deshalbsind dazu einige Erlauterungen zum Parameterkonzept sinnvoll.

Es gibt in C++/C mehrere Ubergabeverfahren fur Parameter.

(a) Werteparameter, call by valueDies ist der ubliche Weg, wo dem Parameter beim Aufruf der Funktion eine Variable uber-geben wird. Der Parameter wird dann mit dem Wert dieser Variablen initialisiert und stellteine Kopie der Variablen dar (erfordert zusatzlichen Speicherplatz), mit der die Funktionarbeiten kann. Die ubergebene Variable wird nicht von der aufrufenden Funktion verandert.

(b) Referenzparameter, call by referenceDabei wird als Argument eine Adresse ubergeben, und uber diese Adresse kann die auf-rufende Funktion auf ein Objekt zugreifen. Dieses sogenannte var-Objekt ist sowohl in derFunktion als auch in ihrer Umgebung verfugbar. Vorgenommene Veranderungen in der Funk-tion werden damit auch an die Umgebung zuruckgegeben (Ruckgabeparameter).

Fur Felder ist das Zeigerkonzept und damit die Adressubergabe fest in der Sprache veran-kert. Fur einfache Variablen muss ein Zeiger-Parameter (speichert und verwaltet Adressen)verwendet werden, der auf das var-Objekt zugreift und es auch verandert.

Drei Moglichkeiten der Verwendung von Zeigern als Parameter und damit der Ubergabe vonAdressen in Parametern gibt es.

- Direkte Ubergabe der Adresse einer Variablen (&variable),- analog dazu die Ubergabe eines Zeigers,- Ubergabe einer Referenz (nur C++).

Betrachten wir die parameterabhangige Eingabefunktion

void eingabe(int n, double a[Nmax][Nmax], double b[Nmax])

{ int i,j;

printf("Matrixdimension n = ");

scanf("%d",&n);

printf("Matrixeingabe element- und zeilenweise\n");

for(i=0;i<n;i++)

for(j=0;j<n;j++) scanf("%lf",&a[i][j]);

printf("Rechte Seite elementweise\n");

for(i=0;i<n;i++) scanf("%lf",&b[i]);

}

23

Darin ist die Variable n ein Werteparameter sowie die Felder a und b automatisch Refe-renzparameter. Die Felder stehen nach ihrer Eingabe innerhalb der Funktion auch in derUmgebung (Rahmenprogramm) zur Verfugung. Die eingelesene Dimension n wird jedochnicht zuruckgegeben. Das soll nun verandert werden.

Ubergabeverfahren von Adresse oder Zeiger fur C++/C

Definition der Eingabefunktion

void eingabe(int *n, double a[Nmax][Nmax], double b[Nmax])

{ int i,j;


scanf("%d",n); // genauso scanf("%d",&*n);


for(i=0;i<*n;i++)

for(j=0;j<*n;j++)

scanf("%lf",&a[i][j]);


for(i=0;i<*n;i++)

scanf("%lf",&b[i]);

}

Aufruf der Eingabefunktion

...

int n,abb,t,i,j;

double a[Nmax][Nmax],b[Nmax],x[Nmax];

int ind[Nmax];

...

eingabe(&n,a,b);

...

Ubergabeverfahren von Referenz nur fur C++

Definition der Eingabefunktion

void eingabe(int &n, double a[Nmax][Nmax], double b[Nmax])

{ int i,j;


scanf("%d",&n);


for(i=0;i<n;i++)

for(j=0;j<n;j++)



for(i=0;i<n;i++)

scanf("%lf",&b[i]);

}

Aufruf der Eingabefunktion

...

int n,abb,t,i,j;


int ind[Nmax];

...

eingabe(n,a,b);

...

Bemerkenswert ist die Syntax fur die Funktion mit dem Referenzparameter.Die Funktion sieht der Werteparameter-Version sehr ahnlich. Nur in der Parameterdeklara-tion zeigt das Adresssymbol &, dass es sich um einen Referenzparameter handelt. Auch derAufruf ist identisch.

24

Demoprogramm gauss21.c mit Ubergabeverfahren von Zeiger

(genauso gauss21.cpp mit Ubergabeverfahren von Zeiger,analog gauss22.cpp mit Ubergabeverfahren von Referenz)

/* gauss21.c */


/* Tausch mit betragsgroessten NNE in der Spalte */

/* Ax=b, A => L\U, P, PA=LU */


#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <conio.h>

#define Nmax 10

// Ohne Prototypen, alle Definitionen stehen vor der main-Funktion


{ int i,j;


scanf("%d",n); // genauso richtig ist scanf("%d",&*n);

// Kontrolle

// printf("*n=%d\n",*n); // n

// printf("n=%d\n",n); // Adresse auf n

// printf("&*n=%d\n",&*n); // gleiches Ergebnis

// getch();


for(i=0;i<*n;i++)

for(j=0;j<*n;j++)



for(i=0;i<*n;i++)

scanf("%lf",&b[i]);

}

void ausgabe(int n, double x[Nmax], int ind[Nmax])

{ int i;

printf("Loesungsvektor\n");

for(i=0;i<n;i++)

printf("%lf ",x[i]);

printf("\nPermutationsvektor\n");

for(i=0;i<n;i++)

printf("%d ",1+ind[i]);

printf("\n");

}

void gauss(int n, int *t, double a[Nmax][Nmax], double b[Nmax], double x[Nmax], int ind[Nmax], int *sing)

{

double max,s,eps;

int i,j,m,index,h;

/* Initialisierung */

*sing = 0;

for(m=0;m<n;m++) /* Permutationsvektor */

ind[m] = m;

eps = 1e-12; /* Toleranz fuer Test auf Singularitaet */

/* Hinrechnung */

for(m=0;m<n;m++)

{

/* Pivotsuche */

max = 0.0;

for(i=m;i<n;i++)

{

s = a[i][m];

if(fabs(s)>fabs(max))

{

25

max = s;

index = i;

}

}

/* Test auf Singularitaet */

if(fabs(max)<eps)

{

*sing = 1;

*t = m;

printf("\nMatrix singulaer, Abbruch im Schritt = %d\n",m+1);

goto _LM1;

}

/* Zeilentausch */

if(index>m)

{

h = ind[m];

ind[m] = ind[index];

ind[index] = h;

for(i=0;i<n;i++)

{

s = a[m][i];

a[m][i] = a[index][i];

a[index][i] = s;

}

s = b[m];

b[m] = b[index];

b[index] = s;

}

/* Bestimmung der Spaltenelemente und Zeilenelemente */

for(i=m+1;i<n;i++)

{

s = a[i][m]/max;

a[i][m] = s;

for(j=m+1;j<n;j++)

a[i][j] = a[i][j]-s*a[m][j];

b[i] = b[i]-s*b[m];

}

}

/* Rueckrechnung */

for(i=n-1;i>=0;i--)

{

s = -b[i];

for(j=n-1;j>i;j--)

s += a[i][j]*x[j];

x[i] = -s/a[i][i];

}

_LM1: ;

}

int main()

{ int n,abb,t,i,j;


int ind[Nmax];

printf("Gauss-Algorithmus: Resttableau-Algorithmus\n");

printf("mit Spaltenpivotisierung und Zeilenvertauschung\n\n");

eingabe(&n,a,b);

gauss(n,&t,a,b,x,ind,&abb);

if(abb==1)

{


getch();

return(-1);

}

ausgabe(n,x,ind);

getch();

return(0);

}

26


2 C:\d\Neundorf\nwptexte\tech phy\05\dev cpp\projekt1\gauss21.exe

Gauss-Algorithmus: Resttableau-Algorithmusmit Spaltenpivotisierung und Zeilenvertauschung

Matrixdimension n = 2Matrixeingabe element- und zeilenweise1246Rechte Seite elementweise11

Loesungsvektor-2.000000 1.500000Permutationsvektor2 1

———————————————————————————————————-

Gauss-Algorithmus: Resttableau-Algorithmusmit Spaltenpivotisierung und Zeilenvertauschung

Matrixdimension n = 2Matrixeingabe element- und zeilenweise1212Rechte Seite elementweise11

Matrix singulaer, Abbruch im Schritt = 2Gleichungssystem nicht loesbar!

Im weiteren sollen zwei Programmiervarianten im Vordergrund stehen, namlich die Zeiger-Ubergabe in C, die auch in C++ machbar ist, sowie die Referenz-Ubergabe in C++.

• Als erstes verwenden wir fur die 3 Funktionen eingabe, ausgabe, gauss Prototypen.

C-Version mit Zeiger

/* gauss2a.c */



/* Ax=b, A => L\U, P, PA=LU */


#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <conio.h>

#define Nmax 10

// mit Prototypen, Vorwaertsdeklarationen

void eingabe(int*, double(*)[Nmax], double*); // 1.Par. = Zeiger

void gauss(int, int*, double(*)[Nmax], double*, double*, int*, int*); // 2.,7.Par. = Zeiger

27

void ausgabe(int, double*, int*);

int main()

{ int n,abb,t,i,j;


int ind[Nmax];



eingabe(&n,a,b);


if(abb==1)

{


getch();

return(-1);

}

ausgabe(n,x,ind);

getch();

return(0);

}

// Definitionen stehen nach der main-Funktion


{ int i,j;


scanf("%d",n);


for(i=0;i<*n;i++)

for(j=0;j<*n;j++)



for(i=0;i<*n;i++)

scanf("%lf",&b[i]);

}


{ int i;

printf("\nLoesungsvektor\n");

for(i=0;i<n;i++)



for(i=0;i<n;i++)


printf("\n");

}

void gauss(int n, int *t, double a[Nmax][Nmax], double b[Nmax],

double x[Nmax], int ind[Nmax], int *sing)

{

double max,s,eps;

int i,j,m,index,h;


*sing = 0;

for(m=0;m<n;m++) ind[m] = m; /* Permutationsvektor */


/* Hinrechnung */

for(m=0;m<n;m++)

{ /* Pivotsuche */

max = 0.0;

for(i=m;i<n;i++)

{ s = a[i][m];


{ max = s;

index = i;

}

}


if(fabs(max)<eps)

{ *sing = 1;

*t = m;

28

printf("\nMatrix singulaer, Abbruch im Schritt %d\n",m+1);

goto _LM1;

}

/* Zeilentausch */

if(index>m)

{ h = ind[m];


ind[index] = h;

for(i=0;i<n;i++)

{ s = a[m][i];


a[index][i] = s;

}

s = b[m];

b[m] = b[index];

b[index] = s;

}


for(i=m+1;i<n;i++)

{ s = a[i][m]/max;

a[i][m] = s;

for(j=m+1;j<n;j++)

a[i][j] = a[i][j]-s*a[m][j];

b[i] = b[i]-s*b[m];

}

}

/* Rueckrechnung */

for(i=n-1;i>=0;i--)

{ s = -b[i];

for(j=n-1;j>i;j--)

s += a[i][j]*x[j];

x[i] = -s/a[i][i];

}

_LM1: ;

}

C++-Version mit Referenz

/* gauss2a.cpp */

...


void eingabe(int &, double(*)[Nmax], double*); // 1.Par. = Referenz

void gauss(int, int &, double(*)[Nmax], double*, double*, int*, int &); // 2.,7.Par. = Referenz


int main()

{ int n,abb,t,i,j;


int ind[Nmax];



eingabe(n,a,b);

gauss(n,t,a,b,x,ind,abb);

if(abb==1)

{


getch();

return(-1);

}

ausgabe(n,x,ind);

getch();

return(0);

}

// Definitionen stehen nach der main-Funktion


...


...

void gauss(int n, int &t, double a[Nmax][Nmax], double b[Nmax], double x[Nmax], int ind[Nmax], int &sing)

...

29

• Als zweites erzeugen wir zu den Prototypen fur die 3 Funktionen eingabe, ausgabe, gauss

eine Header-Datei und nehmen diese ins Rahmenprogramm auf.

C-Version: Header-Datei gauss2b.h.

/* gauss2b.h */


#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <conio.h>

#define Nmax 10


void eingabe(int*, double(*)[Nmax], double*); // 1.Par. = Zeiger

void gauss(int, int*, double(*)[Nmax], double*, double*, int*, int*); // 2.,7.Par. = Zeiger


Diese liegt im selben Verzeichnis wie das einschließende Programm gauss2b.c und wirddort verwendet.

/* gauss2b.c */


#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <conio.h>

#include "gauss2b.h" // Header-Datei im aktuellen Verzeichnis

// Definitionen und Prototypen eingabe, gauss, ausgabe

int main()

...

// Definitionen der Prototypen stehen nach der main-Funktion

...

C++-Version: Header-Datei gauss2bp.h.

/* gauss2bp.h */


#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <conio.h>

#define Nmax 10


void eingabe(int &, double(*)[Nmax], double*); // 1.Par. = Referenz

void gauss(int, int &, double(*)[Nmax], double*, double*, int*, int &); // 2.,7.Par. = Referenz


Diese liegt im selben Verzeichnis wie das einschließende Programm gauss2b.cpp und wirddort verwendet.

/* gauss2b.cpp */


#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <conio.h>

#include "gauss2bp.h" // Header-Datei im aktuellen Verzeichnis


int main()

...

// Definitionen der Prototypen stehen nach der main-Funktion

...

30

• Als dritten und letzten Schritt nehmen wir die Header-Datei zu den 3 Prototypen, bildenaus deren Definitionen einen eigenstandigen Quelltext, um schließlich mit dem Rahmenpro-gramm ein Projekt daraus zu konstruieren.

C-Version:

Header-Datei gauss2b.h.Quelltext gauss2cc.c.

/* gauss2cc.c */




// Definitionen der Prototypen als eigenstaendiger Quelltext


...


...


...

Rahmenprogramm (Quelltext mit main-Funktion) gauss2c.c.

/* gauss2c.c */


#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <conio.h>



int main()

...

Projekt gauss2c mitRahmenprogramm gauss2c.c,Quelltext gauss2cc.c,Header-Datei gauss2b.h

Wir erstellen dieses Projekt.Die einzelnen Schritte seien grob umrissen.

1. Zunachst ist es sinnvoll, alle Dateien zu schließen.

Menupunkt Datei

⇒ Schließen oder Alles Schließen

2. Dann erstellt man das neue Projekt.

Menupunkt Datei

⇒ Neu ⇒ Projekt

Auswahlfenster Neues Projekt

Basic ⇒ Empty Project

—– Projekt Optionen: —————————————————Name: • C-Projekt

gauss2c

Ok

31

3. Damit offnet sich ein neues Auswahlfenster Create new project

Speichern in: C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\projekt1

Dateiname: gauss2c SpeichernDateityp: dev-c++ project(*.dev)

Es erscheint unter Projekt im Projektfenster folgender Baum des noch leeren Projektsgauss2c.

ℵ⋃

≡ gauss2c

Gleichzeitig wird im Verzeichnis die Projektdatei gauss2c.dev mit folgendem Inhaltangelegt.

[Project]

filename=C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\projekt1\gauss2c.dev

name=hilf

Folders=

CommandLine=

Order=

[Views]

ProjectView=1

4. Nun vervollstandigt man das Projekt um die Komponenten gauss2c.c und gauss2cc.c.

Menupunkt Projekt

⇒ Zum Projekt hinzufugen

Damit offnet sich das Auswahlfenster Unit offnen

Suchen in: projekt1 (aktuelles Verzeichnis)

Dateiname: gauss2c OffnenDateityp: All known Files

Genauso mit gauss2cc.

Im Projektfenster vergroßert sich der Baum des Projekts gauss2c. Es erscheint

2- –ℵ⋃

≡ gauss2c||− gauss2c.c|− gauss2cc.c

Hinzugefugte Komponenten konnen durch Anklicken als Quelltexte angezeigt werden.Gleichzeitig vergroßert sich im Verzeichnis die Projektdatei gauss2c.dev.

[Project]

FileName=gauss2c.dev

Name=gauss2c

UnitCount=2

Type=1

Ver=1

ObjFiles=

Includes=

Libs=

PrivateResource=

ResourceIncludes=

MakeIncludes=

Resources=

Compiler=

Linker=

IsCpp=0

Icon=

ExeOutput=

ObjectOutput=

OverrideOutput=0

32

OverrideOutputName=

Folders=

CommandLine=

Focused=1

Order=0,1

[Unit1]

FileName=gauss2c.c

Open=1

Folder=

Top=0

CursorCol=1

CursorRow=1

TopLine=1

LeftChar=1

[Unit2]

FileName=gauss2cc.c

Open=1

Folder=

Top=1

CursorCol=1

CursorRow=1

TopLine=1

LeftChar=1

[Views]

ProjectView=1

Das Projekt kann also konfiguriert werden, indem insbesondere uberflussige Knoten(Programme) geloscht und neue Knoten hinzugefugt werden (hierarchische Struktur).

5. Es ist manchmal sinnvoll bzw. notwendig, Verzeichnisse von eingebundenen Dateien(C-Includes, C++-Includes) zu erganzen.

Menupunkt Werkzeuge

⇒ Compiler Optionen

Damit offnet sich das Auswahlfenster Compiler Optionen

Unter Verzeichnisse und C-Includes kann im Fenster das Verzeichnis

C:\d\neundorf\nwptexte\tech_phy\05\dev_cpp\projekt1

eingetragen werden. Dann Hinzufugen und Ok .

6. Man klickt nun auf das Projekt gauss2c und kommt zum

• Menupunkt Ausfuhren

⇒ Kompilieren

Es folgen Informationen zum Ubersetzen und im Kompilier Log steht

Building Makefile: "C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\projekt1\Makefile.win"

Fuehrt make... aus

make.exe -f "C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\projekt1\Makefile.win" all

gcc.exe -c gauss2c.c -o gauss2c.o

-I"C:/Dev-Cpp/include"

-I"c:/d/Neundorf/nwptexte/tech_phy/05/dev_cpp/projekt1" -s

gcc.exe -c gauss2cc.c -o gauss2cc.o



gcc.exe gauss2c.o gauss2cc.o -o "gauss2c.exe"

-L"C:/Dev-Cpp/lib"



Ausfuehrung beendet


33

Im Verzeichnis stehen dann die ausfuhrbare Datei gauss2c.exe sowie die Objektda-teien gauss2c.o und gauss2cc.o.Sind die Objektdateien zu diesem Projekt schon vorhanden (z. B. durch wiederhol-ten Aufruf), so mussen diese nicht erzeugt werden und die Information verkurzt sichentsprechend. Legt man im Verzeichnis jedoch ein weiteres Projekt an, das zufalligdieselben Namen von Objektdateien verwendet, sollte man zuvor die “alten“ Objekteloschen.

Dazu kommt zum Projekt noch das sogenannte Makefile Makefile.win.

# Project: gauss2c

# Makefile created by Dev-C++ 4.9.6.0

CC = gcc.exe

WINDRES = windres.exe

RES =

OBJ = gauss2c.o gauss2cc.o $(RES)

LIBS = -L"C:/Dev-Cpp/lib"

INCS = -I"C:/Dev-Cpp/include" -I"C:/d/Neundorf/nwptexte/tech_phy/05/dev_cpp/projekt1"

BIN = gauss2c.exe

CFLAGS = $(INCS) -s

.PHONY: all all-before all-after clean clean-custom

all: all-before gauss2c.exe all-after

clean: clean-custom

rm -f $(OBJ) $(BIN)

$(BIN): $(OBJ)

$(CC) $(OBJ) -o "gauss2c.exe" $(LIBS) $(CFLAGS)

gauss2c.o: gauss2c.c

$(CC) -c gauss2c.c -o gauss2c.o $(CFLAGS)

gauss2cc.o: gauss2cc.c

$(CC) -c gauss2cc.c -o gauss2cc.o $(CFLAGS)

Das Makefile zeigt an, wie das Projekt aufgebaut ist, wie es abgearbeitet wir, wel-che Abhangigkeiten zwischen einzelnen Programmmodulen bestehen und welche Pro-grammbibliotheken eingebunden sind. Unter Unix bzw. Linux kann es einfach mit demKommando make aufgerufen werden. Unter Windows gibt es dies nicht.


⇒ Ausfuhren

Es folgt das Ausfuhren des Programms mit den Fenster fur Ein- und Ausgaben.





Fuehrt make... aus


make.exe: Nothing to be done for ‘all’.

Ausfuehrung beendet


Diese Kurzform ergibt sich, wenn das Makefile schon vorhanden ist.

34

C++-Version:

Header-Datei gauss2bp.h.Quelltext gauss2cc.cpp.

/* gauss2cc.cpp */




// Definitionen der Prototypen als eigenstaendiger Quelltext


...


...

void gauss(int n, int &t, double a[Nmax][Nmax], double b[Nmax], double x[Nmax], int ind[Nmax], int &sing)

...

Rahmenprogramm (Quelltext mit main-Funktion) gauss2c.cpp.

/* gauss2c.cpp */


#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <conio.h>



int main()

...

Projekt gauss2cpp mitRahmenprogramm gauss2c.cpp,Quelltext gauss2cc.cpp,Header-Datei gauss2bp.h

Der Umgang mit diesem Projekt, das hier dasselbe Verzeichnis benutzen soll, ist sehr ahn-lich zum vorhergehenden Projekt gauss2c. Es kann aber immer nur ein Projekt bearbeitetwerden.

Ein Unterschied besteht darin, dass man es naturlich mit C++ zu tun hat und der Compilerg++ verwendet wird.Dann beachte man die Modifikationen in den Funktionen eingabe, gauss bezuglich derVerwendung von Referenzparametern.

Da die Komponenten des Projekts gauss2cpp dieselben Bezeichner haben wie die vomProjekt gauss2c, ist es ratsam vor der Kompilierung bestehende “alte“ Objektdateiengauss2c.o und gauss2cc.o sicherheitshalber zu loschen.

Im vorliegenden Fall kann jedoch jedes der beiden Projekte jeweils mit den Objektdateiendes anderen arbeiten.

Im Projektfenster steht der Baum des Projekts gauss2cpp.

2- –ℵ⋃

≡ gauss2cpp||− gauss2c.cpp|− gauss2cc.cpp

Klickt man den Schalter Klassen an, so erscheinen im Projektfenster die Klassen (hierFunktionen) main, ausgabe, eingabe, gauss.

35

Variante 3

Demoprogramme gauss31.cpp, gauss32.cpp und gauss31.c

Ahnlich zur Variante 2 wollen wir den Gauß-Algorithmus in Form des verketteten Algorith-mus als kompakte Funktion und die Ein-/Ausgabe als selbstandige Funktionen definieren unddabei die schon beschriebenen Moglichkeiten der Eingabe- und Ruckgabeparameter nutzen.

Demoprogramm gauss31.c mit Ubergabeverfahren von Zeiger

(genauso gauss31.cpp mit Ubergabeverfahren von Zeiger,analog gauss32.cpp mit Ubergabeverfahren von Referenz)

/* gauss31.c */

/* Verketteter Gauss-Algorithmus mit Spaltenpivotisierung und Zeilenvertauschung */


/* Ax=b, A => (-L)\U, P, PA=LU */


#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <conio.h>

#define Nmax 10

// Ohne Prototypen, alle Definitionen stehen vor der main-Funktion


{ int i,j;


scanf("%d",n); // genauso richtig ist scanf("%d",&*n);


for(i=0;i<*n;i++)

for(j=0;j<*n;j++)



for(i=0;i<*n;i++)

scanf("%lf",&b[i]);

}


{ int i;

printf("\nLoesungsvektor\n");

for(i=0;i<n;i++)



for(i=0;i<n;i++)


printf("\n");

}


{

double max,s,eps;

int i,j,m,index,h;


*sing = 0;

for(m=0;m<n;m++) ind[m] = m; /* Permutationsvektor */


/* Hinrechnung */

for(m=0;m<n;m++)

{

/* Pivotsuche */

max = 0.0;

for(i=m;i<n;i++)

{ s = a[i][m];

for(j=0;j<m;j++)

s += a[i][j]*a[j][m];

a[i][m] = s;


{ max = s;

36

index = i;

}

}


if(fabs(max)<eps)

{ *sing = 1;

*t = m;

printf("\nMatrix singulaer, Abbruch im Schritt %d\n",m+1);

goto _LM1;

}

/* Zeilentausch */

if(index>m)

{ h = ind[m];


ind[index] = h;

for(i=0;i<n;i++)

{ s = a[m][i];


a[index][i] = s;

}

s = b[m];

b[m] = b[index];

b[index] = s;

}

/* Bestimmung der Spaltenelemente */

for(i=m+1;i<n;i++)

a[i][m] = -a[i][m]/max;

/* Bestimmung der Zeilenelemente */

for(i=m+1;i<n;i++)

{ s = a[m][i];

for(j=0;j<m;j++)

s += a[m][j]*a[j][i];

a[m][i] = s;

}

s = b[m];

for(j=0;j<m;j++)

s += a[m][j]*b[j];

b[m] = s;

}

/* Rueckrechnung */

for(i=n-1;i>=0;i--)

{ s = -b[i];

for(j=n-1;j>i;j--)

s += a[i][j]*x[j];

x[i] = -s/a[i][i];

}

_LM1: ;

}

int main()

{ int n,abb,t,i,j;


int ind[Nmax];

printf("Verketteter Gauss-Algorithmus\n");


eingabe(&n,a,b);


if(abb == 1)

{ printf("Gleichungssystem nicht loesbar!\n");

getch();

return(-1);

}

ausgabe(n,x,ind);

getch();

return(0);

}

Weitere Betrachtungen kann man entsprechend der Vorgehensweise in Variante 2 machen,so dass auf Einzelheiten verzichtet wird.Im Prinzip findet nur die Ersetzung des Resttableau-Gauß-Algorithmus durch den verkette-ten Gauß-Algorithmus statt.

37

4.2 Aufgabe 2

Aus der Aachener Bibliothek wollen wir den Komplex LGS mit Gauß-Algorithmus auswahlen.Das Demoprogramm TGAUSS.C befindet sich im Verzeichnis \CNUM8\04\TST und um-fasst inhaltlich mehr als den Gauß-Algorithmus.Die dortigen Include-Anweisungen fur die Header-Dateien

#include <basis.h> /*Grundlegender Deklarationsteil, Funktionskoepfe*/

#include <u_proto.h> /*Vordeklaration aller Bibliotheksfunktionen*/

#include <vmblock.h> /*Dynamische Vektoren und Matrizen*/

enthalten schon wichtige Hinweise dafur, welche Files zum Projekt gehoren.Dazu kommen also noch die Dateien

\CNUM8\04\FGAUSS.C /*Direkte Verfahren zur Loesung LGS*/

/*u.a. Gauss-Algorithmus mit Varianten*/

\CNUM8\BASIS\BASIS.C /*grundlegende Funktionen: Definitionsdatei*/

\CNUM8\BASIS\VMBLOCK.C /*Verwaltung von dynamischen Vektoren und Matrizen*/

\CNUM8\04\EIN /*Beispielmatrizen, hier nicht genutzt*/

Man erkennt eventuell die benotigten Files aus den Include-Anweisungen der Anwendungs-beispiele.Der Einfachheit halber kopiert man alle diese Files in das Arbeitsverzeichnis

C:\d\Neundorf\nwptexte\tech phy\05\dev cpp\projekt2

Dabei soll die Kleinschreibung der Bezeichner verwendet werden.

Nun stellt man das C-Projekt gemaß der Vorgehensweise in Abschnitt 4.1 zusammen.

Projekt tgauss mitRahmenprogramm tgauss.c,Quelltexte fgauss.c, basis.c, vmblock.c,Header-Dateien basis.h, u proto.h, vmblock.h

Im Projektfenster steht der Baum des Projekts tgauss.

2- –ℵ⋃

≡ tgauss||− basis.c||− fgauss.c||− tgauss.c|− vmblock.c

Klickt man den Schalter Klassen an, so erscheint im Projektfenster eine umfangreiche Listevon Klassen und Funktionen zum Projekt.

Man klickt nun auf das Projekt tgauss und kommt zum





Fuehrt make... aus


38

gcc.exe -c basis.c -o basis.o

-I"C:/Dev-Cpp/include" -I"C:/d/Neundorf/nwptexte/tech_phy/05/dev_cpp/projekt2" -s

gcc.exe -c vmblock.c -o vmblock.o


gcc.exe -c fgauss.c -o fgauss.o


gcc.exe -c tgauss.c -o tgauss.o


gcc.exe basis.o vmblock.o fgauss.o tgauss.o -o "tgauss.exe"

-L"C:/Dev-Cpp/lib" -I"C:/Dev-Cpp/include"


Ausfuehrung beendet


Im Verzeichnis stehen dann die Projektdatei tgauss.dev, die ausfuhrbare Datei tgauss.exe

sowie die Objektdateien tgauss.o, fgauss.o, basis.o, vmblock.o.


# Project: tgauss


CC = gcc.exe


RES =

OBJ = basis.o vmblock.o fgauss.o tgauss.o $(RES)


INCS = -I"C:/Dev-Cpp/include" -I"C:/d/Neundorf/nwptexte/tech_phy/05/dev_cpp/projekt2"

BIN = tgauss.exe

CFLAGS = $(INCS) -s


all: all-before tgauss.exe all-after

clean: clean-custom

rm -f $(OBJ) $(BIN)

$(BIN): $(OBJ)

$(CC) $(OBJ) -o "tgauss.exe" $(LIBS) $(CFLAGS)

basis.o: basis.c

$(CC) -c basis.c -o basis.o $(CFLAGS)

vmblock.o: vmblock.c

$(CC) -c vmblock.c -o vmblock.o $(CFLAGS)

fgauss.o: fgauss.c

$(CC) -c fgauss.c -o fgauss.o $(CFLAGS)

tgauss.o: tgauss.c

$(CC) -c tgauss.c -o tgauss.o $(CFLAGS)

Folgende Berechnungen konnen im Programm durchgefuhrt werden:

– inverse Matrix,– transponierte Inverse,– Gauß-Algorithmus mit rechter Seite als Einheitsvektoren,

dazu Determinante und Kontrolle der Spur der Matrix AA−1,– Gauß-Algorithmus fur LGS mit Spaltenpivotisierung und Zeilenvertauschung,

dazu LU -Faktorisierung von PA (zeilenvertauschte Matrix) und Permutationsvektor,– Gauß-Algorithmus fur LGS mit Nachiteration zur Verbesserung der Naherungslosung.

39

Beispiel Ax = b

A =

(

1 24 6

)

, b =

(

11

)

, x∗ = A−1b =

(

−21.5

)

det(A) = −2, A =

(

1 04 1

)(

1 20 −2

)

, trace(A) = spur(A) =n∑

i=1

aii,

A−1 =

(

−3 12 −0.5

)

, A−T =

(

−3 21 −0.5

)

p = (2, 1), P =

(

0 11 0

)

, PA =

(

4 61 2

)

=

(

1 00.25 1

)(

4 60 0.5

)

= LU

Eingabe- und Ergebnisfenster

Haltepunkte wie getch(); sind bei der Ausgabe großer Datenmengen als auch nach Fehler-meldungen sinvoll.

2 C:\d\Neundorf\nwptexte\tech phy\05\dev cpp\projekt2\tgauss.exe——————————————————————————————–Gauss Method——————————————————————————————–

Dimension of the input quadratic matrix n = 2Input matrix (elementwise and rowwise)1246

Dimension of the input matrix = 2

Input matrix:

1.000000 2.0000004.000000 6.000000

Transposed Inverse:

-3.000000 2.0000001.000000 -0.500000

——————————————————————————————–Gauss Method for multiple right sides——————————————————————————————–

Right sides = identity matrixInverse:

-3.000000 1.0000002.000000 -0.500000

Determinante = -2.000000e+000Trace (Matrix * Inverse) = 2.0000000000000000e+000

40

——————————————————————————————–Solution without iterative improvement, one right side vector——————————————————————————————–

Input vector (right side elementwise):11

A:1.000000 2.0000004.000000 6.000000

b:1.000000 1.000000

LU:4.000000 6.0000000.250000 0.500000

perm:1 0

Solution:x: -2.000000 1.500000

——————————————————————————————–Solution with iterative improvement——————————————————————————————–

Input vector (right side elementwise):11

A:1.000000 2.0000004.000000 6.000000

b:1.000000 1.000000

LU:4.000000 6.0000000.250000 0.500000

perm:1 0

Solution

x: -2.000000 1.500000

——————————————————————————————–Solution with post iteration for right sides = identity vectors——————————————————————————————–

b: 1.000000 0.000000x: -3.000000 2.000000

b: 0.000000 1.000000x: 1.000000 -0.500000

——————————————————————————————–

41

4.3 Aufgabe 3

In Abschnitt 1, Aufgabe 5, haben wir das Newtonsche Naherungsverfahren fur Nullstellenvon skalaren Funktionen als Demoprogramm tnewton1.c vorgestellt.Die Definitionen und Prototypen, also die Prozedurkopfe usw., findet man in den Header-Files basis.h, u_proto.h, tfunc1.h als Vorabdeklarationen und werden damit bekanntgemacht.Ihre eigentliche Definition erfolgt in den Funktionen fnewtonc.c, tfunc1.c, basis.c.Dazu brauchen wir die Dateien

\CNUM8\02\TST\TNEWTON.C /*Test program for newton*/

\CNUM8\02\TST\TFUNC1.C /*Definition of test functions for newton, pegasus,..*/

/*1. and 2. derivity*/

\CNUM8\02\FNEWTON.C /*Numerische Verfahren zur Loesung NLG: NV,...*/

\CNUM8\BASIS\BASIS.C /*grundlegende Funktionen: Definitionsdatei*/

Daraus machen wir nun das C-Projekt tnewton im Verzeichnis

C:\d\Neundorf\nwptexte\tech phy\05\dev cpp\projekt3,

in das alle Dateien (Bezeichner klein geschrieben) kopiert werden.

Projekt tnewton mitRahmenprogramm tnewton.c,Quelltexte tfunc1.c, fnewton.c, basis.c,Header-Dateien basis.h, u proto.h, tfunc1.h

Im Projektfenster steht der Baum des Projekts tnewton.

2- –ℵ⋃

≡ tnewton||− basis.c||− fnewton.c||− tfunc1.c|− tnewton.c

Man klickt nun auf das Projekt tnewton und kommt zum



Es folgen Informationen zum Ubersetzen, es entstehen im Verzeichnis die entsprechendenDateien *.o, *.exe und im Kompilier Log stehtBuilding Makefile: "C:\d\Neundorf\nwptexte\tech_phy\05\dev_cpp\projekt3\Makefile.win"

Fuehrt make... aus


gcc.exe -c basis.c -o basis.o


gcc.exe -c fnewton.c -o fnewton.o


gcc.exe -c tfunc1.c -o tfunc1.o


gcc.exe -c tnewton.c -o tnewton.o


gcc.exe basis.o fnewton.o tfunc1.o tnewton.o -o "tnewton.exe"


-I"C:/d/Neundorf/nwptexte/tech_phy/05/dev_cpp/projekt3" -s

Ausfuehrung beendet


42


2 C:\d\Neundorf\nwptexte\tech phy\05\dev cpp\tnewton1.exe——————————————————————————–Newton Method for real, nonlinear functions——————————————————————————–

f = 5∗x-exp(x)Start value x0 = 2.000000Return code = 1Root = 2.5426413577735265Function value = -3.448110e-014Iterations = 7

f = (((((x-6)∗x+15)∗x-20)∗x+15)∗x-6)∗x+1Start value x0 = 1.000000Return code = 0Root = 1.0000000000000000Function value = 0.000000e+000Iterations = 0

f = sin(x)Start value x0 = 2.000000Return code = 0Root = 3.1415926535897931Function value = 1.224606e-016Iterations = 6

f = 1+sin(x)Start value x0 = 2.000000Return code = 0Root = 4.7123890100448769Function value = 4.398608e-016Iterations = 26

f = exp(x)-(1.0+x+x∗x∗0.5)Start value x0 = 1.000000Return code = 0Root = 0.0000143482588343Function value = 3.555099e-016Iterations = 170

f = (x-1.0)∗(x-1.0)∗(sin(PI∗x)-log(fabs(2.0∗x/(x+1.0)))Start value x0 = 1.000000Return code = 0Root = 1.0000000000000000Function value = 0.000000e+000Iterations = 0

——————————————————————————–

43

4.4 Aufgabe 4

Nun soll ein etwas umfangreicheres Projekt zur Losung LGS mit den Gauß-Algorithmus zu-sammengestellt werden.

Der zentrale Datentyp ist ein Zeiger auf ein zweidimensionales Feld (Matrix) mit Kompo-nenten vom Typ double. Vektoren werden in dieser Konvention als Feld mit einer Spalteinterpretiert. Der Datentyp laMatrix zeigt auf eine Struktur, die als Informationen dieZeilen- und Spaltenmdimension sowie die Matrixelemente zeilenweise enthalt.Die Nummerierung der Zeilen und Spalten beginnt in C/C++ jeweils bei Null.

Der Resttableau-Gauß-Algorithmus mit Spaltenpivotisierung und Zeilenvertauschung ausAbschnitt 4.1 (vergl. Demoprogramm gauss21.c mit Ubergabeverfahren von Zeiger) mussauf diesen Datentyp angepasst werden.Die Koeffizientenmatrix A ist Eingabematrix und wird uberschrieben durch die LU -Zerlegungvon PA, welche dann zuruckgegeben wird. Deshalb werden Kopien von A und b angelegt.Der Losungsvektor x und der Permutationsvektor von p=perm fur die Zeilenvertauschungensind einspaltige reelle Ergebnismatrizen. Die Komponenten von p sind wertemaßig naturlicheZahlen, was in der formatierten Ausgabe berucksichtigt wird.Falls der Gauß-Algorithmus nicht durchfuhrbar ist, so weist der Ergebnisparameter sing6= 0darauf hin, wobei mit der Große t noch die Stufe (Schritt im Resttableau-Algorithmus) desAbbruchs angezeigt werden kann. Gleichzeitig ist der Stand der Berechnungen bis zum even-tuellen Abbruch abfragbar.Die Toleranz ε=eps fur den Test auf Singularitat ist in der Funktion definiert und gegebe-nenfalls zu verandern.

void gauss4c(int n, int *t, laMatrix *A, laMatrix b, laMatrix *x, laMatrix *perm, int *sing)

{

double max,s,eps,h;

int i,j,k,index;


*sing = 0;

*t = 0;

for(k=0;k<n;k++)

{

(*perm)[k][0] = k; /* Permutationsvektor */

(*x)[k][0] = 0;

}


/* Hinrechnung */

for(k=0;k<n;k++)

{

/* Pivotsuche */

max = 0;

for(i=k;i<n;i++)

{

s = fabs((*A)[i][k]);

if(s>max)

{

max = s;

index = i;

}

}


if(max<eps)

{

*sing = 1;

*t = k+1;

printf("\nMatrix singulaer, Abbruch im Schritt = %d\n",k+1);

goto _LM1;

}

44

/* Zeilentausch */

if(index>k)

{

h = (*perm)[k][0]; /* laVertauscheZe_(*perm,k,index); */

(*perm)[k][0] = (*perm)[index][0];

(*perm)[index][0] = h;

/*

for(i=0;i<n;i++)

{

s = (*A)[k][i];

(*A)[k][i] = (*A)[index][i];

(*A)[index][i] = s;

}

*/

laVertauscheZe_(*A,k,index);

s = b[k][0]; /* laVertauscheZe_(b,k,index); */

b[k][0] = b[index][0];

b[index][0] = s;

}


max = (*A)[k][k];

for(i=k+1;i<n;i++)

{

s = (*A)[i][k]/max;

(*A)[i][k] = s;

for(j=k+1;j<n;j++)

(*A)[i][j] = (*A)[i][j]-s*(*A)[k][j];

b[i][0] = b[i][0]-s*b[k][0];

}

}

/* Rueckrechnung */

for(i=n-1;i>=0;i--)

{

s = -b[i][0];

for(j=n-1;j>i;j--)

s += (*A)[i][j]*(*x)[j][0];

(*x)[i][0] = -s/(*A)[i][i];

}

_LM1: for(k=0;k<n;k++)

(*perm)[k][0] = 1+(*perm)[k][0];

}

In der Funktion gauss4c erkennt man den Funktionsaufruf laVertauscheZe_(*A,k,index)

zum Tausch der Pivotzeile “nach oben“.Diese, noch weitere Module zur linearen Algebra und die Zeigerdefinition sollen in kleinenProgrammbibliotheken zusammengefasst werden.

Alle Vorwartsdeklarationen und eine Sammlung von Prototypen befinden sich im Header-Filela_.h und sind im Wesentlichen selbstkommentierend.

// la_.h

/* Header-File fuer eine Bibliothek fuer Lineare Algebra */

#include <stdlib.h>

#include <stdio.h>

typedef double **laMatrix;

void laGlSysEinlesen_(laMatrix *A, laMatrix *b);

/* Liest die Matrix A und den Vektor b der rechten Seite

aus der Datei "Matrix.in" bzw. Eingabedatei ... ein.

Vektor = Spaltenvektor = n*1-Matrix */

void laGlSysAusgeben_(laMatrix A, laMatrix b);

/* Gibt die Matrix A und den Vektor b auf dem Bildschirm aus. */

45

void laMatAusgeben_(laMatrix A);

/* Gibt die Matrix A auf dem Bildschirm aus. */

void laVektAusgeben_(laMatrix x);

/* Gibt den reellen Vektor x (n*1-Matrix) auf dem Bildschirm aus. */

void laVektAusgebenI_(laMatrix x);

/* Gibt den reellen Vektor x mit ganzzahligen Werten der Komponenten auf dem Bildschirm aus. */

laMatrix laMatrixNeu_(int n, int m);

/* Erzeugt eine n*m-Matrix und speichert die Dimensionen in ihr ab. */

void laVernichte_(laMatrix A);

/* Gibt den Speicher der Matrix A wieder frei. */

void laDim_(laMatrix A, int *n, int *m);

/* Gibt die Dimensionen der Matrix A in n (Zeilen) und m (Spalten) zurueck. */

laMatrix laNull_(int n, int m);

/* Erzeugt eine n*m-Nullmatrix */

laMatrix laEins_(int n, int m);

/* Erzeugt eine n*m-Einheitsmatrix mit 1 auf der Diagonale */

void laMult_(laMatrix A, laMatrix B, laMatrix *C);

/* Multipliziert die Matrix A mit der Matrix B algebraisch und speichert das Ergebnis

in einer von der Funktion erzeugten Matrix ab. Auf diese Matrix zeigt nach dem Aufruf C.

Fehlermeldung und Programmabbruch, falls die inneren Dimensionen von A und B

nicht vertraeglich sind. */

void laAdd_(laMatrix A, laMatrix B, laMatrix *C);

/* Addiert die Matrizen A und B und speichert das Ergebnis in einer

von der Funktion erzeugten Matrix ab. Auf diese Matrix zeigt nach dem Aufruf C.

Fehlermeldung und Programmabbruch, falls die Dimensionen von A und B nicht vertraeglich sind. */

void laSub_(laMatrix A, laMatrix B, laMatrix *C);

/* Subtrahiert die Matrix B von der Matrix A und speichert das Ergebnis in einer

von der Funktion erzeugten Matrix ab. Auf diese Matrix zeigt nach dem Aufruf C.

Fehlermeldung und Programmabbruch, falls die Dimensionen von A und B nicht vertraeglich sind. */

void laSkMult_(double z, laMatrix A);

/* Multipliziert die reelle Zahl z mit der Matrix A. */

laMatrix laTransponiere_(laMatrix A);

/* Transponiert die Matrix A. Da die Matrix i. Allg. nicht quadratisch ist, muss

auch der Speicher umorganisiert werden. A zeigt also nach dem Aufruf auf einen

anderen Speicherbreich, der urspruengliche Speicherbereich wird freigegeben. */

void laVertauscheZe_(laMatrix A, int r, int s);

/* Vertauscht in der Matrix A die Zeile r mit Zeile s.

Fehlermeldung und Programmabbruch, falls r, s unzulaessig sind.*/

void laVertauscheSp_(laMatrix A, int r, int s);

/* Vertauscht in der Matrix A die Spalte r mit Spalte s.

Fehlermeldung und Programmabbruch, falls r, s unzulaessig sind.*/

void laKopiere_(laMatrix A, laMatrix B);

/* Kopiert die Matrix A in eine vor der Funktion erzeugte Matrix.

Auf diese Matrix zeigt nach dem Aufruf B. Dimensionskontrolle mit Fehlermeldung. */

void laKopiere1_(laMatrix A, laMatrix *B);

/* Kopiert die Matrix A in eine von der Funktion erzeugte Matrix.

Auf diese Matrix zeigt nach dem Aufruf B. */

void laElemMult_(laMatrix A, laMatrix B, laMatrix *C);

/* Multiplikation der Matrizen A und B elementweise und Speichern des Ergebnisses

in eine von der Funktion erzeugte Matrix. Auf diese Matrix zeigt nach dem Aufruf C.

Fehlermeldung und Programmabbruch, falls die Dimensionen von A und B

nicht vertraeglich sind. */

void laKreuzProd_(laMatrix A, laMatrix B, laMatrix *C);

/* Berechnet das Kreuzprodukt der 3*1-Matrizen (Vektoren) A und B und speichert das Ergebnis

in einer von der Funktion erzeugte 3*1-Matrix. Auf diese Matrix zeigt nach dem Aufruf C.

Bei unzulaessigen Dimensionen von A, B Abbruch mit einer Fehlermeldung. */

46

double laSkalProd_(laMatrix A, laMatrix B);

/* Berechnet das Skalarprodukt der n*1-Matrizen (Vektoren) A und B und gibt den Wert zurueck.

Bei unzulaessigen Dimensionen Abbruch mit Fehlermeldung. */

double laNorm_(char c, laMatrix A);

/* Berechnet die Norm von n*m-Matrizen bzw. n*1-Matrizen (Vektoren)

und gibt den Wert zurueck.

Matrix Vektor

Norm ’u’ : Zeilensummennorm, ’u’: Maximumnorm

Norm ’1’ : Spaltensummennorm, ’1’: Betragssummennorm

Norm ’F’ : Frobeniusnorm, ’2’: euklidische Norm

Bei unzulaessiger Normwahl Fehlermeldung und Programmabbruch. */

Die Definitionen sind in verschiedenen Quelltexten untergebracht.Zunachst ist es das sehr wichtige Modul la_.c mit vier Grundfunktionen, wobei die Funk-tion laTransponiere_ hier keine Anwendung findet.

// la_.c

/* Eine Bibliothek fuer Lineare Algebra */

/* Definitionen von laMatrixNeu_, laVernichte_, laDim_, laTransponiere_ */

#include <stdlib.h>

#include <stdio.h>

#include "la_.h"

laMatrix laMatrixNeu_(int n, int m)

{

unsigned i;

laMatrix A;

if((n<1)||(m<1))

{

printf("Fehler: Dimension(en) unzulaessig\n");

getch();

exit(EXIT_FAILURE);

}

if((A = malloc(sizeof(double*)*n+2*sizeof(int)))==NULL)

{

fprintf(stderr,"Zu wenig Speicher\n");

getch();

exit(EXIT_FAILURE);

}

*((int*)A) = n;

*(((int*)A)+1) = m;

A = (double**)(((int*)A)+2);

for(i=0;i<n;i++)

if((*(A+i) = malloc(sizeof(double)*m))==NULL)

{


getch();

exit(EXIT_FAILURE);

}

return(A);

}

/*----------------------------------------------------------------------*/

void laVernichte_(laMatrix A)

{

unsigned i,n;

n = *(((int*)A)-2);

for(i=0;i<n;i++)

free(*(A+i));

free(((int*)A)-2);

}

/*----------------------------------------------------------------------*/

void laDim_(laMatrix A, int *n, int *m)

{

*n = *(((int*)A)-2);

*m = *(((int*)A)-1);

}

47

/*----------------------------------------------------------------------*/

laMatrix laTransponiere_(laMatrix A)

{

int n,m,i,j;

laMatrix B;

laDim_(A,&n,&m);

B = laMatrixNeu_(m,n);

for(i=0;i<=n;++i)

for(j=0;j<=m;++j)

B[j][i] = A[i][j];

laVernichte_(A);

return(B);

}

Die Funktionen fur die Ein- und Ausgabe sind problemangepasst und im Modul laInOut_.c

zusammengefasst. Sie werden alle genutzt.Dort aufgerufen werden auch die Funktionen laMatrixNeu_, laDim_.

// laInOut_.c

/* Eine Bibliothek fuer Lineare Algebra, Ein/Ausgaben */

/* Definitionen von laGlSysEinlesen_, laGlSysAusgeben_, laMatAusgeben_, laVektAusgeben_, laVektAusgebenI_ */

#include <stdlib.h>

#include <stdio.h>


#include "la_.h"

/*----------------------------------------------------------------------------*/

void laGlSysEinlesen_(laMatrix *A, laMatrix *b)

{

FILE *matIn;

int i,j,n,m;

char s[80]="Matrix.in"; // Initialisierung

printf("Dateiname fuer Daten A,b zum LGS Ax=b : ");

scanf("%s",&s); // auch scanf("%s",s);

matIn = fopen(s,"r"); // matIn = fopen("Matrix.in","r");

if(matIn==NULL)

{

printf("Fehler: Datei nicht gefunden\n");

getch();

exit(EXIT_FAILURE);

}

fscanf(matIn,"%u %u",&n,&m);

*A = laMatrixNeu_(n,m);

*b = laMatrixNeu_(n,1);

for(i=0;i<n;i++)

for(j=0;j<m;j++)

fscanf(matIn,"%lf",&(*A)[i][j]);

for(i=0;i<n;i++)

fscanf(matIn,"%lf",&(*b)[i][0]);

fclose(matIn);

}

/*----------------------------------------------------------------------------*/

void laGlSysAusgeben_(laMatrix A, laMatrix b)

{

int i,j,n,m;

laDim_(A,&n,&m);

for(i=0;i<n;i++)

{

for(j=0;j<m;j++)

printf(" %10.6lf ",A[i][j]);

printf(" | %10.6lf \n",b[i][0]);

}

}

48

void laMatAusgeben_(laMatrix A)

{

int i,j,n,m;

laDim_(A,&n,&m);

for(i=0;i<n;i++)

{

for(j=0;j<m;j++)

printf(" %10.6lf ",A[i][j]);

printf("\n");

}

}

/*----------------------------------------------------------------------------*/

void laVektAusgeben_(laMatrix x)

{

int i,n,m;

laDim_(x,&n,&m);

for(i=0;i<n;i++)

printf(" %10.6lf ",x[i][0]);

printf("\n");

}

/*----------------------------------------------------------------------------*/

void laVektAusgebenI_(laMatrix x)

{

int i,n,m;

laDim_(x,&n,&m);

for(i=0;i<n;i++)

printf(" %2.0lf ",x[i][0]);

printf("\n");

}

Die Beispiele werden als ASCII-Dateien in der Datenreihenfolge n,m, A(1..n, 1..m), b(1..n), ...bereitgestellt. Der Dateiname fur die Daten A, b zum LGS wird eingelesen. Als Testvariantenzeigen wir LGS mit regularen und singularen Koeffizientenmatrizen(reg.: matrix23, matrix3, matrix40, sing.: matrix57, matrix58).

matrix23 matrix3 matrix40 matrix57 matrix58

-------- ------------------- ---------------- ------------ --------

3 3 2 2 5 5 3 3 3 3

1 -2 2 1.441969 1.040807 2 -3 1 -2 4 1 -2 2 2 2 2

-1 1 -1 1.040807 0.751250 -4 6 -2 5 -6 -1 2 -1 1 3 0

-2 -2 1 0.401162 6 -9 3 -4 10 -2 4 1 2 4 1

1 0.289557 2 -4 3 2 -3 1 5

-1 Loesung: -2 5 -3 2 -1 0 4

-3 1 7 3 7

Loesung: -1 -10 Loesung: Loesung:

1 EW: 17 1 keine

1 2.193218999999544 5 1 EW:

1 4.5595081932092e-13 1 1 0,1,5

EW: Loesung: EW:

1,1,1 14,14,13,0,2 0,2+-i*2.646

EW: 0.914

5.5713+-i*8.844

-0.0285+-i*0.489

Von den anderen selbstandigen Definitionen benotigen wir eigentlich nur die FunktionenlaKopiere_, laVertauscheZe_ und wegen einiger Kontrollrechnungen auch laKopiere1_,laAdd_, laNorm_. Aber es sollen jedoch alle zum Projekt hinzugefugt werden.

49

// laKopiere_.c, B = A

#include "la_.h"

void laKopiere_(laMatrix A, laMatrix B)

{

int n,m,r,s,i,j;

laDim_(A,&n,&m);

laDim_(B,&r,&s);

if((n!=r)||(m!=s))

{

printf("Fehler: Matrizen sind nicht von gleicher Dimension\n");

getch();

exit(EXIT_FAILURE);

}

for(i=0;i<n;i++)

for(j=0;j<m;j++)

B[i][j] = A[i][j];

}

// laVertauscheZe_.c, Zeilenvertauschung in Matrix, Zeilen r<->s

#include "la_.h"

void laVertauscheZe_(laMatrix A, int r, int s)

{

int n,m,j;

double h;

laDim_(A,&n,&m);

if((r<0)||(r>n-1)||(s<0)||(s>n-1))

{

printf("Fehler: Zeilennummer(n) unzulaessig\n");

getch();

exit(EXIT_FAILURE);

}

for(j=0;j<m;j++)

{

h = A[s][j];

A[s][j] = A[r][j];

A[r][j] = h;

}

}

---------------------------------------------------------------------

// laKopiere1_.c, B = A

#include "la_.h"

void laKopiere1_(laMatrix A, laMatrix *B)

{

int n,m,i,j;

laDim_(A,&n,&m);

*B = laMatrixNeu_(n,m);

for(i=0;i<n;i++)

for(j=0;j<m;j++)

(*B)[i][j] = A[i][j];

}

// laAdd_.c, C = A+B

#include "la_.h"

void laAdd_(laMatrix A, laMatrix B, laMatrix *C)

{

int n,m,k,l,i,j;

laDim_(A,&n,&m);

laDim_(B,&k,&l);

if((n!=k)||(m!=l))

{

printf("Fehler: Matrixdimension(en) nicht gleich\n");

getch();

exit(EXIT_FAILURE);

}

*C = laMatrixNeu_(n,m);

for(i=0;i<n;i++)

50

for(j=0;j<m;j++)

(*C)[i][j] = A[i][j]+B[i][j];

}

// laNorm_.c, Berechnet die Norm von n*m-Matrizen bzw. n*1-Matrizen (Vektoren)

#include <math.h>

#include "la_.h"

double laNorm_(char cc, laMatrix A)

{

int n,m,i,j;

double no,h;

no = 0;

if((cc!=’u’)&&(cc!=’1’)&&(cc!=’2’)&&(cc!=’F’))

{

printf("Fehler: Falsche Norm\n");

getch();

exit(EXIT_FAILURE);

}

laDim_(A,&n,&m);

// Fallunterscheidung m=1 und m>1

if(m==1)

switch (cc)

{

case ’u’: for(i=0;i<n;i++)

if (fabs(A[i][0])>no)

no = fabs(A[i][0]);

break;

case ’1’: for(i=0;i<n;i++)

no += fabs(A[i][0]);

break;

case ’2’: for(i=0;i<n;i++)

no += A[i][0]*A[i][0];

no = sqrt(no);

break;

}

else

switch (cc)

{

case ’u’: for(i=0;i<n;i++)

{

h = 0;

for(j=0;j<m;j++)

h += fabs(A[i][j]);

if (h>no) no = h;

}

break;

case ’1’: for(j=0;j<m;j++)

{

h = 0;

for(i=0;i<n;i++)

h += fabs(A[i][j]);

if (h>no) no = h;

}

break;

case ’F’: for(i=0;i<n;i++)

for(j=0;j<m;j++)

no += A[i][j]*A[i][j];

no = sqrt(no);

break;

}

return(no);

}

---------------------------------------------------------------------

// laEins_.c, Erzeugung der Einheitsmatrix

#include "la_.h"

laMatrix laEins_(int n, int m)

{

int i,min;

laMatrix M;

51

M = laNull_(n,m);

if(n>m) min = m; else min = n;

for (i=0;i<min;i++)

M[i][i] = 1;

return M;

}

// LaElemMult_.c, C = A*B elementweise

#include "la_.h"

void laElemMult_(laMatrix A, laMatrix B, laMatrix *C)

{

int n,m,k,l,i,j;

laDim_(A,&n,&m);

laDim_(B,&k,&l);

if((n!=k)||(m!=l))

{

printf("Fehler: Matrizen sind nicht von gleicher Dimension\n");

getch();

exit(EXIT_FAILURE);

}


for(i=0;i<n;i++)

for(j=0;j<m;j++)

(*C)[i][j] = A[i][j]*B[i][j];

}

// laKreuzProd_.c, Kreuzprodukt zweier 3*1-Matrizen (3-dimensionale Vektoren)

#include "la_.h"

void laKreuzProd_(laMatrix A, laMatrix B, laMatrix *C)

{

int na,ma,nb,mb;

if((A==*C)||(B==*C))

{

fprintf(stderr,"Fehler: Gleiche Zeiger\n");

getch();

exit(EXIT_FAILURE);

}

laDim_(A,&na,&ma);

laDim_(B,&nb,&mb);

if((na!=3)||(ma!=1)||(nb!=3)||(mb!=1))

{

fprintf(stderr,"In A,B Zeilenanzahl<>3, Spaltenanzahl<>1\n");

exit(EXIT_FAILURE);

}

*C = laMatrixNeu_(3,1);

(*C)[0][0] = A[1][0]*B[2][0]-A[2][0]*B[1][0];

(*C)[1][0] =-A[0][0]*B[2][0]+A[2][0]*B[0][0];

(*C)[2][0] = A[0][0]*B[1][0]-A[1][0]*B[0][0];

}

// laMult_.c, C = AB algebraisch

#include "la_.h"

void laMult_(laMatrix A, laMatrix B, laMatrix *C)

{

int i,j,k,n,m,r,s;

double h;

laDim_(A, &n, &m);

laDim_(B, &r, &s);

if(m!=r)

{

printf ("Fehler: innere Dimensionen unzulaessig\n");

getch();

exit(EXIT_FAILURE);

}

*C = laMatrixNeu_(n,s);

for(i=0;i<n; i++)

for(j=0;j<s;j++)

{

52

h = 0;

for(k=0;k<m;k++)

h += A[i][k]*B[k][j];

(*C)[i][j] = h;

}

}

// laNull_.c, Erzeugung der Nullmatrix

# include "la_.h"

laMatrix laNull_(int n, int m)

{

int i,j;

laMatrix M;

if((n<1)||(m<1))

{


getch();

exit(EXIT_FAILURE);

}

M = laMatrixNeu_(n,m);

for (i=0;i<n;i++)

for (j=0;j<m;j++)

M[i][j] = 0;

return M;

}

// laSkalProd_.c, Skalarprodukt der n*1-Matrizen (Vektoren) A und B

#include "la_.h"

double laSkalProd_(laMatrix A, laMatrix B)

{

int n,m,k,l,i;

double erg;

laDim_(A,&n,&m);

laDim_(B,&k,&l);

if((n==k)&&(m==1)&&(l==1))

{

erg=0;

for(i=0;i<n;i++)

erg += A[i][0]*B[i][0];

return (erg);

}

printf("Fehler: Falsche Dimension(en)\n");

getch();

exit(EXIT_FAILURE);

}

// laSkMult_.c, Multiplikation der Zahl z mit der Matrix A

#include "la_.h"

void laSkMult_(double z, laMatrix A)

{

int n,m,i,j;

laDim_(A,&n,&m);

for(i=0;i<n;i++)

for(j=0;j<m;j++)

A[i][j] = A[i][j]*z;

}

// laSub_.c, C = A-B

#include "la_.h"

void laSub_(laMatrix A, laMatrix B, laMatrix *C)

{

int n,m,k,l,i,j;

laDim_(A,&n,&m);

laDim_(B,&k,&l);

if((n!=k)||(m!=l))

{

53

printf("Fehler: Matrixdimension(en) nicht gleich\n");

getch(); exit(EXIT_FAILURE);

}


for(i=0;i<n;i++)

for(j=0;j<m;j++)

(*C)[i][j] = A[i][j]-B[i][j];

}

// laVertauscheSp_.c, Spaltenvertauschung in Matrix, Spalten r<->s

#include "la_.h"

void laVertauscheSp_(laMatrix A, int r, int s)

{

int n,m,i;

double h;

laDim_(A,&n,&m);

if((r<0)||(r>m-1)||(s<0)||(s>m-1))

{

printf("Fehler: Spaltennummer(n) unzulaessig\n");

getch(); exit(EXIT_FAILURE);

}

for(i=0;i<n;i++)

{

h = A[i][r];

A[i][r] = A[i][s];

A[i][s] = h;

}

}

Nun fehlt nur noch der Quelltext (Rahmenprogramm) meinprog_.c mit der main-Funktion.In ihm befinden sich die notwendigen include-Anweisungen, die Funktion gauss4c, dieEingaben, Berechnungen, Ausgaben sowie gewisse Kontrollrechnungen zwecks Uberprufungder Richtigkeit einiger Funktionen.

// meinprog_.c

#include<stdlib.h>

#include<stdio.h>

#include<math.h>

//#include<conio.h>

#include "la_.h"

void gauss4c(int n, int *t, laMatrix *A, laMatrix b, laMatrix *x, laMatrix *perm, int *sing)

...

int main()

{

int n,m,t,sing;

double erg;

laMatrix A,b,x,LU,c,perm,C1,C2,C3;

char kont,prot;

printf("Resttableau-Gauss-Algorithmus mit Spaltenpivotisierung und Zeilenvertauschung\n");

printf("Tausch mit betragsgroessten NNE in der Spalte\n");

printf("Ax=b, A => L\\U, P, PA=LU \n");

printf("Ux=c Loesung mit Rueckwaertseinsetzen\n");

printf("-----------------------------------------------------------------------------\n\n");

laGlSysEinlesen_(&A,&b);

printf("Kontrollausgabe A,b (j/n)? ");

kont = getch();

printf("\n");

if(kont==’j’)

{

laDim_(A,&n,&m);

printf("\nDimension A(n,m): n=%i, m=%i\n",n,m);

printf("\nA | b\n");

laGlSysAusgeben_(A,b);

printf("-----------------------------------------------------------------------------\n\n");

getch();

}

54

laDim_(A,&n,&m);

if(n!=m)

{

printf("\nDimension n<>m, keine Berechnung");

getch();

return (-1);

}

LU = laMatrixNeu_(n,m);

laKopiere_(A,LU);

c = laMatrixNeu_(n,1);

laKopiere_(b,c);

x = laMatrixNeu_(n,1);

perm = laMatrixNeu_(n,1);

gauss4c(n,&t,&LU,c,&x,&perm,&sing);

if(sing==0)

{

printf("Loesungsvektor x\n");

laVektAusgeben_(x);

printf("-----------------------------------------------------------------------------\n");

getch();

}

printf("\nErgebnisprotokoll (j/n)? ");

prot = getch();

if(prot==’j’)

{

printf("\n---------------------------------------------------------------------------\n");

if (sing==1)

printf("Achtung! Vorzeitiger Abbruch im Schritt t = %i\n",t);

printf("Dimension A(n,m): n=%i, m=%i\n",n,m);

printf("\nA | b\n");

laGlSysAusgeben_(A,b);

if(sing==0)

{

printf("\nLoesungsvektor x\n");

laVektAusgeben_(x);

}

printf("\nPermutationsvektor perm\n");

laVektAusgebenI_(perm);

printf("\nZerlegung PA=L\\U\n");

laMatAusgeben_(LU);

printf("-----------------------------------------------------------------------------\n");

getch();

}

printf("\n\nKontrollrechnungen (j/n)? ");

kont = getch();

if(kont==’j’)

{

printf("\nC1 = LU mittels laKopiere_\n");

C1 = laMatrixNeu_(n,m);

laKopiere_(LU,C1);

laMatAusgeben_(C1);

printf("\nC2 = LU mittels laKopiere1_\n");

laKopiere1_(LU,&C2);

laMatAusgeben_(C2);

printf("\nC3 = A+C2\n");

laAdd_(A,C2,&C3);

laMatAusgeben_(C3);

printf("\nNormen von b\n");

laVektAusgeben_(b);

printf("||b||_u=%10.6lf, ||b||_1=%10.6lf, ||b||_2=%10.6lf\n",laNorm_(’u’,b),laNorm_(’1’,b),laNorm_(’2’,b));

printf("\nNormen von A\n");

laMatAusgeben_(A);

printf("||A||_u=%10.6lf, ||A||_1=%10.6lf, ||A||_F=%10.6lf\n",laNorm_(’u’,A),laNorm_(’1’,A),laNorm_(’F’,A));

getch();

laVernichte_(C1); laVernichte_(C2); laVernichte_(C3);

}

// Speicherplatz freigeben

laVernichte_(A); laVernichte_(LU);

laVernichte_(b); laVernichte_(c); laVernichte_(x); laVernichte_(perm);

}

55

Zusammenstellung des C-Projekts gauss4c im Arbeitsverzeichnis

C:\d\Neundorf\nwptexte\tech phy\05\dev cpp\projekt4

Im Projektfenster steht der Baum des Projekts. Die wirklich verwendeten Module sind fetthervorgehoben.

2- –ℵ⋃

≡ gauss4c||− la .c||− laAdd .c||− laEins_.c||− laElemMult_.c||− laInOut .c||− laKopiere .c||− laKopiere1 .c||− laKreuzProd_.c||− laMult_.c||− laNorm .c||− laNull_.c||− laSkalProd_.c||− laSkMult_.c||− laSub_.c||− laVertauscheSp_.c||− laVertauscheZe .c|− meinprog .c

Klickt man den Schalter Klassen an, so erscheint im Projektfenster eine umfangreiche Listevon Klassen und Funktionen zum Projekt.

Man klickt nun auf das Projekt gauss4c und kommt zum





Fuehrt make... aus


gcc.exe -c la_.c -o la_.o


gcc.exe -c laAdd_.c -o laAdd_.o


gcc.exe -c laEins_.c -o laEins_.o


gcc.exe -c laElemMult_.c -o laElemMult_.o


gcc.exe -c laInOut_.c -o laInOut_.o


gcc.exe -c laKopiere_.c -o laKopiere_.o


gcc.exe -c laKopiere1_.c -o laKopiere1_.o


gcc.exe -c laKreuzProd_.c -o laKreuzProd_.o


gcc.exe -c laMult_.c -o laMult_.o


gcc.exe -c laNorm_.c -o laNorm_.o

56


gcc.exe -c laNull_.c -o laNull_.o


gcc.exe -c laSkalProd_.c -o laSkalProd_.o


gcc.exe -c laSkMult_.c -o laSkMult_.o


gcc.exe -c laSub_.c -o laSub_.o


gcc.exe -c laVertauscheSp_.c -o laVertauscheSp_.o


gcc.exe -c laVertauscheZe_.c -o laVertauscheZe_.o


gcc.exe -c meinprog_.c -o meinprog_.o


gcc.exe la_.o laAdd_.o laEins_.o laElemMult_.o laInOut_.o laKopiere_.o laKopiere1_.o laKreuzProd_.o

laMult_.o laNorm_.o laNull_.o laSkalProd_.o laSkMult_.o laSub_.o laVertauscheSp_.o laVertauscheZe_.o

meinprog_.o -o "gauss4c.exe"


-I"C:/d/Neundorf/nwptexte/tech_phy/05/dev_cpp/projekt4" -s

Ausfuehrung beendet


Im Verzeichnis stehen dann die Projektdatei gauss4c.dev, die ausfuhrbare Dateigauss4c.exe sowie alle Objektdateien *.o der Projektquellen.


# Project: gauss4c


CC = gcc.exe


RES =

OBJ = la_.o laAdd_.o laEins_.o laElemMult_.o laInOut_.o laKopiere_.o laKopiere1_.o laKreuzProd_.o

laMult_.o laNorm_.o laNull_.o laSkalProd_.o laSkMult_.o laSub_.o laVertauscheSp_.o laVertauscheZe_.o

meinprog_.o $(RES)


INCS = -I"C:/Dev-Cpp/include"

BIN = gauss4c.exe

CFLAGS = $(INCS) -s


all: all-before gauss4c.exe all-after

clean: clean-custom

rm -f $(OBJ) $(BIN)

$(BIN): $(OBJ)

$(CC) $(OBJ) -o "gauss4c.exe" $(LIBS) $(CFLAGS)

la_.o: la_.c

$(CC) -c la_.c -o la_.o $(CFLAGS)

laAdd_.o: laAdd_.c

$(CC) -c laAdd_.c -o laAdd_.o $(CFLAGS)

laEins_.o: laEins_.c

$(CC) -c laEins_.c -o laEins_.o $(CFLAGS)

laElemMult_.o: laElemMult_.c

$(CC) -c laElemMult_.c -o laElemMult_.o $(CFLAGS)

laInOut_.o: laInOut_.c

$(CC) -c laInOut_.c -o laInOut_.o $(CFLAGS)

laKopiere_.o: laKopiere_.c

$(CC) -c laKopiere_.c -o laKopiere_.o $(CFLAGS)

laKopiere1_.o: laKopiere1_.c

$(CC) -c laKopiere1_.c -o laKopiere1_.o $(CFLAGS)

laKreuzProd_.o: laKreuzProd_.c

$(CC) -c laKreuzProd_.c -o laKreuzProd_.o $(CFLAGS)

laMult_.o: laMult_.c

$(CC) -c laMult_.c -o laMult_.o $(CFLAGS)

laNorm_.o: laNorm_.c

$(CC) -c laNorm_.c -o laNorm_.o $(CFLAGS)

laNull_.o: laNull_.c

57

$(CC) -c laNull_.c -o laNull_.o $(CFLAGS)

laSkalProd_.o: laSkalProd_.c

$(CC) -c laSkalProd_.c -o laSkalProd_.o $(CFLAGS)

laSkMult_.o: laSkMult_.c

$(CC) -c laSkMult_.c -o laSkMult_.o $(CFLAGS)

laSub_.o: laSub_.c

$(CC) -c laSub_.c -o laSub_.o $(CFLAGS)

laVertauscheSp_.o: laVertauscheSp_.c

$(CC) -c laVertauscheSp_.c -o laVertauscheSp_.o $(CFLAGS)

laVertauscheZe_.o: laVertauscheZe_.c

$(CC) -c laVertauscheZe_.c -o laVertauscheZe_.o $(CFLAGS)

meinprog_.o: meinprog_.c

$(CC) -c meinprog_.c -o meinprog_.o $(CFLAGS)

Obwohl die Bibliotheken Dev-Cpp/lib und Dev-Cpp/include mit ihren zahlreichen Header-Files ins Projekt eingebunden sind, ist es notwendig, in den entsprechenden Modulen nochdie Include-Anweisungen zu belassen, also

#include<stdlib.h> // fuer EXIT_FAILURE

#include<stdio.h> // fuer stdin, stderr, ...

#include<math.h> // fuer mathematische Standardfunktionen

Beispiel 1 Ax = b, A regular

A =

1 −2 2−1 1 −1−2 −2 1

, b =

1−1−3

, eindeutige Losung x∗ = A−1b =

111

,

p = (3, 1, 2), P =

0 0 11 0 00 1 0

(Permutationsmatrix),

PA =

−2 −2 1

1 −2 2

−1 1 −1

=

1 0 0

−1

21 0

1

2−2

31

−2 −2 1

0 −3 5

2

0 0 1

6

= LU.


Haltepunkte wie getch() sind bei der Ausgabe großer Datenmengen als auch nach Fehler-meldungen sinvoll.

2 C:\d\Neundorf\nwptexte\tech phy\05\dev cpp\projekt4\gauss4c.exe

Resttableau-Gauss-Algorithmus mit Spaltenpivotisierung und ZeilenvertauschungTausch mit betragsgroessten NNE in der SpalteAx=b, A => L\U, P, PA=LUUx=c Loesung mit Rueckwaertseinsetzen————————————————————————————————————

Dateiname fuer Daten A,b zum LGS Ax=b : Matrix23Kontrollausgabe A,b (j/n)?Dimension A(n,m): n=3, m=3A | b

1.000000 -2.000000 2.000000 | 1.000000-1.000000 1.000000 -1.000000 | -1.000000-2.000000 -2.000000 1.000000 | -3.000000

————————————————————————————————————

Loesungsvektor x1.000000 1.000000 1.000000

58

Fortsetzung


Ergebnisprotokoll (j/n)?————————————————————————————————————Dimension A(n,m): n=3, m=3

A | b1.000000 -2.000000 2.000000 | 1.000000-1.000000 1.000000 -1.000000 | -1.000000-2.000000 -2.000000 1.000000 | -3.000000

Loesungsvektor x1.000000 1.000000 1.000000

Permutationsvektor perm3 1 2

Zerlegung PA=L\U-2.000000 -2.000000 1.000000-0.500000 -3.000000 2.5000000.500000 -0.666667 0.166667

————————————————————————————————————

Kontrollrechnungen (j/n)?C1 = LU mittels laKopiere

-2.000000 -2.000000 1.000000-0.500000 -3.000000 2.5000000.500000 -0.666667 0.166667

C2 = LU mittels laKopiere1-2.000000 -2.000000 1.000000-0.500000 -3.000000 2.5000000.500000 -0.666667 0.166667

C3 = A+C2-1.000000 -4.000000 3.000000-1.500000 -2.000000 1.500000-1.500000 -2.666667 1.166667

Normen von b1.000000 -1.000000 -3.000000

||b|| u= 3.000000, ||b|| 1= 5.000000, ||b|| 2= 3.316625

Normen von A1.000000 -2.000000 2.000000-1.000000 1.000000 -1.000000-2.000000 -2.000000 1.000000

||A|| u= 5.000000, ||A|| 1= 5.000000, ||A|| F= 4.582576

Fur die Ergebnisdarstellung bei großeren LGS bzw. mit mehr Mantissenstellen der Dezimal-zahlen muss man geeignete Veranderungen an der Ein- und Ausgabefunktionen vornehmen.

59

Beispiel 2 Ax = b, A regular, aber det(A) = 10−12 und damit sehr schlechte Konditionder Systemmatrix

A =

(

1.441969 1.040807

1.040807 0.751250

)

, b =

(

0.401162

0.289557

)

, Losung x∗ = A−1b =

(

1

−1

)

,

p = (1, 2), P =

(

1 00 1

)


PA = A =

(

1441969

1000000

1040807

1000000

1040807

1000000

751250

1000000

)

=

(

1 01040807

144196911

)(

1441969

1000000

1040807

1000000

0 1

1441969000000

)

=

=

(

1 0

0.7217956835410470 1

)(

1.441969 1.040807

0 0.693496 10−12

)

= LU.





1.441969 1.040807 | 0.4011621.040807 0.751250 | 0.289557

————————————————————————————————————

Loesungsvektor x1.000058 -1.000080

————————————————————————————————————

Ergebnisprotokoll (j/n)?————————————————————————————————————Dimension A(n,m): n=2, m=2

A | b1.441969 1.040807 | 0.4011621.040807 0.751250 | 0.289557

Loesungsvektor x1.000058 -1.000080

Permutationsvektor perm1 2

Zerlegung PA=L\U1.441969 1.0408070.721796 0.000000

————————————————————————————————————

60

Fortsetzung



1.441969 1.0408070.721796 0.000000

C2 = LU mittels laKopiere11.441969 1.0408070.721796 0.000000

C3 = A+C22.883938 2.0816141.762603 0.751250

Normen von b0.401162 0.289557

||b|| u= 0.401162, ||b|| 1= 0.690719, ||b|| 2= 0.494747

Normen von A1.441969 1.0408071.040807 0.751250

||A|| u= 2.482776, ||A|| 1= 2.482776, ||A|| F= 2.193219

Wahlt man im Gauß-Algorithmus gauss4c() beim Test auf die Singularitat der Matrix A dieToleranz zu grob, z. B. ε ≥ 10−12, dann wird wegen u22 = 0.693496 10−12 < ε die Matrix alssingular eingestuft und das Programm abgebrochen. Ein Weiterrechnen in den Resttableausware bei zusatzlicher Spaltenvertauschung durchaus moglich. Dabei konnte man dann auchden Matrixrang ermitteln. Aber die Singularitat als solche bleibt naturlich bestehen.



Dateiname fuer Daten A,b zum LGS Ax=b : Matrix3Kontrollausgabe A,b (j/n)?Dimension A(n,m): n=2, m=2

A | b1.441969 1.040807 | 0.4011621.040807 0.751250 | 0.289557

————————————————————————————————————

Matrix singulaer, Abbruch im Schritt = 2

Ergebnisprotokoll (j/n)?————————————————————————————————————Achtung! Vorzeitiger Abbruch im Schritt t = 2...

61


A =

2 −3 1 −2 4−4 6 −2 5 −6

6 −9 3 −4 102 −4 3 2 −3

−2 5 −3 2 −1

, b =

7−10

1751

, Losung x∗ = A−1b =

14141302

,

p = (3, 5, 4, 2, 1), P =

0 0 1 0 00 0 0 0 10 0 0 1 00 1 0 0 01 0 0 0 0


PA=

6 −9 3 −4 10

−2 5 −3 2 −1

2 −4 3 2 −3

−4 6 −2 5 −6

2 −3 1 −2 4

=

1 0 0 0 0

−1

31 0 0 0

1

3−1

21 0 0

−2

30 0 1 0

1

30 0 −2

71

6 −9 3 −4 10

0 2 −2 2

3

7

3

0 0 1 11

3−31

6

0 0 0 7

3

2

3

0 0 0 0 6

7

=LU.

Eingabe- und Ergebnisfenster (Auswahl)



Dateiname fuer Daten A,b zum LGS Ax=b : Matrix40Kontrollausgabe A,b (j/n)?Dimension A(n,m): n=5, m=5

A | b2.000000 -3.000000 1.000000 -2.000000 4.000000 | 7.000000-4.000000 6.000000 -2.000000 5.000000 -6.000000 | -10.0000006.000000 -9.000000 3.000000 -4.000000 10.000000 | 17.0000002.000000 -4.000000 3.000000 2.000000 -3.000000 | 5.000000-2.000000 5.000000 3.000000 -3.000000 2.000000 | 1.000000

————————————————————————————————————

Loesungsvektor x14.000000 14.000000 13.000000 0.000000 2.000000

————————————————————————————————————...Permutationsvektor perm

3 5 4 2 1

Zerlegung PA=L\U6.000000 -9.000000 3.000000 -4.000000 10.000000-0.333333 2.000000 -2.000000 0.666667 2.3333330.333333 -0.500000 1.000000 3.666667 -5.166667-0.666667 0.000000 0.000000 2.333333 0.6666670.333333 0.000000 0.000000 -0.285714 0.857143

62

Beispiel 4 Ax = b, A singular, rang(A) = 2

A =

1 −2 2−1 2 −1−2 4 1

, b =

103

,

eine der Losungen ist x∗ =

111

,

die einparametrige Losungsmenge ist x∗ =

2c − 1c

1

, c reell.

Abbruch des Gauß-Algorithmus im zweiten Schritt, bis dahin wird folgender “Zwischenstand“erreicht.

p = (3, 2, 1), P =

0 0 10 1 01 0 0

,

L =

1 0 01

2∗ 0

−1

2∗ ∗

, U =

−2 4 1

0 0 −3

2

0 0 5

2

.





1.000000 -2.000000 2.000000 | 1.000000-1.000000 2.000000 -1.000000 | 0.000000-2.000000 4.000000 1.000000 | 3.000000

————————————————————————————————————


63

Fortsetzung


Ergebnisprotokoll (j/n)?————————————————————————————————————Achtung! Vorzeitiger Abbruch im Schritt t = 2Dimension A(n,m): n=3, m=3

A | b1.000000 -2.000000 2.000000 | 1.000000-1.000000 2.000000 -1.000000 | 0.000000-2.000000 4.000000 1.000000 | 3.000000


Zerlegung PA=L\U-2.000000 4.000000 1.0000000.500000 0.000000 -1.500000-0.500000 0.000000 2.500000

————————————————————————————————————


-2.000000 4.000000 1.0000000.500000 0.000000 -1.500000-0.500000 0.000000 2.500000

C2 = LU mittels laKopiere1-2.000000 4.000000 1.0000000.500000 0.000000 -1.500000-0.500000 0.000000 2.500000

C3 = A+C2-1.000000 2.000000 3.000000-0.500000 2.000000 -2.500000-2.500000 4.000000 3.500000

Normen von b1.000000 0.000000 3.000000

||b|| u= 3.000000, ||b|| 1= 4.000000, ||b|| 2= 3.162278

Normen von A1.000000 -2.000000 2.000000-1.000000 2.000000 -1.000000-2.000000 4.000000 1.000000

||A|| u= 7.000000, ||A|| 1= 8.000000, ||A|| F= 6.000000

64


A =

2 2 21 3 02 4 1

, b =

547

, keine Losung x∗.

Abbruch des Gauß-Algorithmus, der ohne Zeilenvertauschung ablauft, im dritten Schritt.Bis dahin wird folgender “Zwischenstand“ erreicht.

L =

1 0 01

21 0

1 1 1

, U =

2 2 2

0 2 −1

0 0 0

.

Die Zerlegung A = LU ist moglich, da erst im letzten Schritt die Komponente u33 = 0 wird.





2.000000 2.000000 2.000000 | 5.0000001.000000 3.000000 0.000000 | 4.0000002.000000 4.000000 1.000000 | 7.000000

————————————————————————————————————


Ergebnisprotokoll (j/n)?————————————————————————————————————Achtung! Vorzeitiger Abbruch im Schritt t = 3Dimension A(n,m): n=3, m=3

A | b2.000000 2.000000 2.000000 | 5.0000001.000000 3.000000 0.000000 | 4.0000002.000000 4.000000 1.000000 | 7.000000


Zerlegung PA=L\U2.000000 2.000000 2.0000000.500000 2.000000 -1.0000001.000000 1.000000 0.000000

65

Wir zeigen einige Modifikationen der bisherigen C-Projekts gauss4c in Bezug auf den Ver-zicht des Befehls getch() zur Zeicheneingabe ohne Bildschirmecho (bei C automatischdurch Einbinden der Bibliotheken vorhanden, bei C++ im Header-File conio.h).Dies betrifft mehrere Funktionen des Projekts. Die Funktion laGlSysEinlesen_ im Quell-text laInOut_c soll dazu als Musterbeispiel ausgewahlt werden.

Bisherige Version:

...


{

FILE *matIn;

int i,j,n,m;



scanf("%s",&s);

matIn = fopen(s,"r");

if(matIn == NULL)

{


getch();

exit(EXIT_FAILURE);

}




for(i=0;i<n;i++)

for(j=0;j<m;j++)


for(i=0;i<n;i++)


fclose(matIn);

}

Mit der scanf-Eingabeanweisung wird eine Zeichenkette (Stream, String) fur den Namen der

Datei eingelesen. Die Eingabe schließt mit der <Enter>-Taste ( ◭yp -Taste) ab. Alle Zeichen

stehen zunachst im Eingabepuffer (Tastaturpuffer), um dann in die Variable s ubertragen zuwerden. Ist der Name unzulassig, kommt die Fehlerausschrift auf dem Monitor. Diese bleibtbis zu einer Tastenbetatigung angezeigt.

Nun sollen zwei Dinge berucksichtigt werden.

1. Damit man sicher ist, dass der Eingabepuffer nach der Namenseingabe auf jeden Fall wie-der leer und damit auch das <Enter>-Zeichen nach der Zeichenketteneingabe geloscht wird,erganzt man die Anweisung (Methode) fflush(stdin) mit stdin aus Header-File stdio.h.

2. Den getch()-Befehl ersetzt man durch die scanf("%c",&ret)-Anweisung mit der char-Variablen ret.

Version (a): im Projekt gauss4ca

...


{

FILE *matIn;

int i,j,n,m;


char ret;


scanf("%s",&s);

fflush(stdin);

66


if(matIn == NULL)

{


scanf("%c",&ret);

exit(EXIT_FAILURE);

}

...

}

Eine weitere Veranderung bedeutet den Verzicht auf die Anweisung fflush(stdin) durchdas Ersetzen mit der scanf()-Anweisung. Dabei werden Dateiname (Zeichenkette) und ab-schließende <Enter>-Taste durch die scanf("%s%c",&s,&ret)-Eingabeanweisung korrektverarbeitet und aus dem Eingabepuffer herausgelesen.

Version (b): im Projekt gauss4cb

...


{

FILE *matIn;

int i,j,n,m;


char ret;


scanf("%s%c",&s,&ret);


if(matIn == NULL)

{


scanf("%c",&ret);

exit(EXIT_FAILURE);

}




for(i=0;i<n;i++)

for(j=0;j<m;j++)


for(i=0;i<n;i++)


fclose(matIn);

}

Ausgehend vom C-Projekt in der Version gauss4cb soll es nun als C++-Projekt dargestelltwerden und funktionieren.Bis auf die Speicherplatzzuweisung fur den Zeigertyp typedef double **laMatrix imHeader-File la_.h klappt alles.Nur die Anweisungen zur Allokation von dynamischen Speicher

if ((A = malloc(sizeof(double*)*n+2*sizeof(int)))==NULL)

...

if ((*(A+i) = malloc(sizeof(double)*m))==NULL)

...

in der Funktion laMatrix laMatrixNeu_(int n, int m) aus dem Quelltext la_.cpp

fuhren zur Fehlermeldung

C:\d\Neundorf\...\projekt4b\la_.cpp [Warning] In function ‘double ** laMatrixNeu_(int, int)’:

23 C:\d\Neundorf\...\projekt4b\la_.cpp ANSI C++ forbids implicit conversion from ‘void *’ in assignment

32 C:\d\Neundorf\...\projekt4b\la_.cpp ANSI C++ forbids implicit conversion from ‘void *’ in assignment

67

C++ kann bei der Allokation von Speicher mit dem Befehl malloc arbeiten. Aber die allo-kierte Variable muss als Zeigertyp mit einem Datentyp der Objekte stehen. Da man hier garnicht uber den Zeiger auf das referenzierte Objekt zugreifen will, sondern nur Speicheradres-sen und damit Adressen von Objekten beliebigen Typs verwaltet, ist ein Zeiger auf void(generischer Zeiger) zu verwenden. Das ist durch den vorangestellten Typ (void*) einfachzu machen.Weitere Veranderungen sind nicht erforderlich.Um spater auf die referenzierten Objekte zugreifen zu konnen, muss der Zeiger dann aberzuerst in einen Zeiger auf den entsprechenden Datentyp des Objekts umgewandelt werden(Dereferenzierung, explizite Konvertierung der Adresse von void* in datentyp*).

Version (c1): im Projekt gauss4cpp1


{

unsigned i;

char ret;

laMatrix A;

if((n<1)||(m<1))

{


scanf("%c",ret);

exit(EXIT_FAILURE);

}

if(((void*)A = malloc(sizeof(double*)*n+2*sizeof(int)))==NULL)

{


scanf("%c",ret);

exit(EXIT_FAILURE);

}

*((int*)A) = n;

*(((int*)A)+1) = m; // auch *((int*)A+1) = m;

A = (double**)(((int*)A) +2); // auch A = (double**)((int*)A +2);

for(i=0;i<n;i++)

if(((void*)(*(A+i)) = malloc(sizeof(double)*m))==NULL)

{


scanf("%c",ret);

exit(EXIT_FAILURE);

}

return(A);

}

Die elegantere C++-Version zur dynamischen Speicherallokation arbeitet mit dem Anwei-sungspaar (Operatoren) new und delete zur Reservierung bzw. Freigabe von Speicher.Man hat jedoch eine besondere Situation bei der Adressierung sowie jeweils Felder vonObjekten. Der Anfangszeiger zeigt zunachst auf eine Objekt (Bereich), in dem die beidenMatrixdimensionen als integer-Zahlen n und m sowie die n Zeiger auf die Matrixzeilen mitjeweils m double-Komponenten stehen.Nachdem die Dimensionen an den ersten beiden “Stellen“ eingetragen wurden, wird der An-fangszeiger um zwei Positionen weitergeruckt. Damit steht er am Beginn des eigentlichenAdressbereiches der Matrix mit ihren double-Elementen. Nun kann jedem Zeilenzeiger seinSpeicherbereich entsprechend der Matrixzeile “zugewiesen“ werden.Der Zugriff auf die Matrixdimensionen geht dabei nicht verloren, denn mit dem Anfangszei-ger kann nach Umtypisierung auf int* auch wieder ruckwarts gezahlt werden.

68

Version (c2): im Projekt gauss4cpp2


{

unsigned i;

char ret;

laMatrix A;

if((n<1)||(m<1))

{


scanf("%c",ret);

exit(EXIT_FAILURE);

}

if(((char*)A = new char[sizeof(double*)*n+2*sizeof(int)])==NULL)

{


scanf("%c",ret);

exit(EXIT_FAILURE);

}

*((int*)A) = n;

*(((int*)A)+1) = m; // auch *((int*)A+1) = m;

A = (double**)(((int*)A) +2); // auch A = (double**)((int*)A +2);

for(i=0;i<n;i++)

if((*(A+i) = new double[m])==NULL)

{


scanf("%c",ret);

exit(EXIT_FAILURE);

}

return(A);

}

/*----------------------------------------------------------------------------*/

void laVernichte_(laMatrix A)

{

unsigned i,n;

n = *(((int*)A)-2);

for(i=0;i<n;i++)

delete [](*(A+i));

delete [](char*)(((int*)A) -2); // auch delete [](char*)((int*)A -2);

}

Schema der dynamischen Speicherverwaltung zu gauss4cpp2

laMatrix

↓ ւ(char*)A = new char[sizeof(double*)*n+2*sizeof(int)]

A →int n

int m ւ*(A+i) = new double[m]

A →double∗ → 1. Matrixzeile mit m double-Elementen

double∗ → 2. Matrixzeile mit m double-Elementen

......

double∗ → n. Matrixzeile mit m double-Elementen

↑

A = (double**)((int*)A +2);

69

4.5 Aufgabe 5

Ein weiteres Projekt zur Behandlung des Algorithmus der Vektoriteration und inversen Vek-toriteration soll folgen.Dabei geht es um die iterative Bestimmung des einzigen betragsgroßten bzw. betragskleinstenEigenwerts einer Matrix sowie des zugehorigen Eigenvektors (partielles Eigenwertproblem).Der zentrale Datentyp ist wieder ein Zeiger auf ein zweidimensionales Feld (Matrix) mitKomponenten vom Typ double. Vektoren werden in dieser Konvention als Feld mit einerSpalte interpretiert. Der Datentyp laMatrix zeigt auf eine Struktur, die als Informationendie Zeilen- und Spaltenmdimension sowie die Matrixelemente zeilenweise enthalt.Die Nummerierung der Zeilen und Spalten beginnt in C/C++ jeweils bei Null.

Die Vektoriteration braucht in jedem Iterationsschritt die Matrix-Vektor-Multiplikation undwegen der Skalierung die Berechnung der Vektornorm.Die inverse Vektoriteration hat anstelle der Matrix-Vektor-Multiplikation die Losung einesLGS. Dabei wird jedoch der bisher verwendete Resttableau-Gauß-Algorithmus mit Spalten-pivotisierung und Zeilenvertauschung

gauss4c(int n,int *t,laMatrix *A,laMatrix b,laMatrix *x,laMatrix *perm,int *sing)

aus Abschnitt 4.4 zerlegt in zwei Prozeduren zur LU -Faktorisierung der Matrix A mit Pivo-strategie (PA = LU, P Permutationsmatrix) sowie zur Losung von Ax = b mittels gestaf-felter LGS gemaß

Ax = b, PAx = LUx = L(Ux) = Pb, Lz = Pb, Ux = z.

Die Matrix A ist Eingabematrix und wird uberschrieben durch die LU -Zerlegung von PA,welche vor der Iteration einmal mit dem Aufwand O(n3) gemacht wird und dann zuruckge-geben wird. Deshalb wird in der inversen Vektoriteration eine Kopie von A angelegt.Der Permutationsvektor von p=perm fur die Zeilenvertauschungen ist eine einspaltige re-elle Ergebnismatrix. Die Komponenten von p sind wertemaßig naturliche Zahlen, die zurVertauschung im Vektor der rechten Seite gebraucht werden. Falls die Faktorisierung nichtdurchfuhrbar ist, so weist der Ergebnisparameter sing6= 0 darauf hin, wobei mit der Großet noch die Stufe (Schritt im Faktorisierungs-Algorithmus) des Abbruchs angezeigt werdenkann.

In jedem Iterationsschritt der inversen Vektoriteration wird dann mit der Faktorisierung furdie jeweils neue rechte Seite eine Vorwarts- und Ruckwartsrechnung mit dem Aufwand O(n2)vorgenommen. Damit ist diese Vorgehensweise naturlich effizienter als die Losung eines LGSmittels Gauß-Algorithmus in jedem Schritt.Die inversen Vektoriteration ist nicht durchfuhrbar, wenn A singular ist. Dann ware derbetragsgroßte Eigenwert von A−1 “unendlich“.

Fur die Iteration braucht man einen nicht verschwindenden Startvektor, der wahlweise ein-gelesen bzw. als Einsvektor oder erster Einheitsvektor definiert wird.

In allen Fallen wird die Iteration gesteuert durch eine maximale Iterationsanzahl maxiterund die Toleranz ε=eps fur den Test der Ubereinstimmung aufeinander folgender Iterier-ter (Betrage des Eigenwerts). Das Vorzeichen des gesuchten Eigenwerts wird anschließendnoch aus dem Vorzeichenvergleich einer geeigneten Komponente der beiden letzten Iterati-onsvektoren bestimmt. Der Iterationsvektor wird skaliert und ist selber Naherung fur denEigenvektor.

70

Zunachst die beiden Funktionen zur LU -Faktorisierung der Matrix A mit Pivostrategie sowiezur Losung von Ax = b mittels gestaffelter LGS.

void gauss_LU(int n, int *t, laMatrix *A, laMatrix *perm, int *sing)

{

double max,s,eps,h;

int i,j,k,index;

char ret;


*sing = 0;

*t = 0;

for(k=0;k<n;k++)

(*perm)[k][0] = k; /* Permutationsvektor */


/* Hinrechnung */

for(k=0;k<n;k++)

{

/* Pivotsuche */

max = 0;

for(i=k;i<n;i++)

{

s = fabs((*A)[i][k]);

if(s>max)

{

max = s;

index = i;

}

}


if(max<eps)

{

*sing = 1;

*t = k+1;

printf("\nMatrix singulaer, Abbruch im Schritt = %d\n",k+1);

scanf("%c",&ret);

exit(EXIT_FAILURE);

}

/* Zeilentausch */

if(index>k)

{

h = (*perm)[k][0]; /* laVertauscheZe_(*perm,k,index); */

(*perm)[k][0] = (*perm)[index][0];

(*perm)[index][0] = h;

laVertauscheZe_(*A,k,index);

}


max = (*A)[k][k];

for(i=k+1;i<n;i++)

{

s = (*A)[i][k]/max;

(*A)[i][k] = s;

for(j=k+1;j<n;j++)

(*A)[i][j] = (*A)[i][j]-s*(*A)[k][j];

}

}

}

/*--------------------------------------------------------------------------------*/

void gestaffelte_LGS(int n, laMatrix LU, laMatrix b, laMatrix *x, laMatrix perm)

{

double h;

int i,j;

laMatrix bp,z;

/* Permutation der rechten Seite */

bp = laMatrixNeu_(n,1);

for (i = 0; i < n; i++)

bp[i][0] = b[(int)perm[i][0]][0]; // mit einfacher Umtypisierung

// bp[i][0] = b[(int)floor(perm[i][0]+1e-15)][0]; // sichere Variante

/* Hinrechnung Lz=bp */

z = laMatrixNeu_(n,1);

z[0][0] = bp[0][0];

71

for(i=1;i<n;i++)

{

h = -bp[i][0];

for(j=0;j<i;j++)

h += LU[i][j]*z[j][0];

z[j][0] = -h;

}

/* Rueckrechnung Ux=z */

for(i=n-1;i>=0;i--)

{

h = -z[i][0];

for(j=n-1;j>i;j--)

h += LU[i][j]*(*x)[j][0];

(*x)[i][0] = -h/LU[i][i];

}

laVernichte_(bp);

laVernichte_(z);

}

Die Vektoriteration und inverse Vektoriteration notieren wir ebenfalls als Funktionen.Dabei sollen die Bibliotheksfunktionen wie

– Initialisieren und Loschen von Feldern (Matrizen und Vektoren),– Kopieren von Feldern,– Normberechnung eines Vektors,– Multiplikation eines Vektors mit einer skalaren Große,– Matrix-Vektor-Multiplikation,– Skalarprodukt zweier Vektoren,

die in la_.h aufgelistet sind, so weit wie moglich Anwendung finden.

void vektorit(int n, laMatrix A, int maxiter, double eps, laMatrix *x, double *la, int *iter)

{

double la_a,la_n,no;

int i,j;

laMatrix y,yn,yh;


*iter = 0;

y = laMatrixNeu_(n,1);

yn = laMatrixNeu_(n,1);

laKopiere_(*x,yn);

no = laNorm_(’2’,yn); // sqrt(laSkalProd_(yn,yn));

la_n = 0;

i = 0;

do

{

la_a = la_n;

laSkMult_(1.0/no,yn);

laKopiere_(yn,y);

laMult_(A,y,&yh);

laKopiere_(yh,yn);

i = i+1;

no = laNorm_(’2’,yn);

la_n = no;

} while((fabs(la_a-la_n)>eps)&&(i<maxiter));

for(j=0;j<n;j++)

if (y[j][0]*yn[j][0]!=0) break;

if(y[j][0]*yn[j][0]<0) la_n = -la_n;


/* Ergebnisse*/

laKopiere_(yn,*x);

*iter = i;

*la = la_n;

laVernichte_(y);

laVernichte_(yn);

}

/*--------------------------------------------------------------------------------*/

72

void inv_vektorit(int n, laMatrix A, int maxiter, double eps, laMatrix *x, double *la, int *iter)

{

double la_a,la_n,no;

int i,j,t,sing;

char ret;

laMatrix y,yn,LU,perm;


*iter = 0;

y = laMatrixNeu_(n,1);

yn = laMatrixNeu_(n,1);

LU = laMatrixNeu_(n,n);

perm = laMatrixNeu_(n,1);

laKopiere_(A,LU);

gauss_LU(n,&t,&LU,&perm,&sing);

if(sing!=0)

{

printf("\nMatrix singulaer -> Abbruch");

scanf("%c",&ret);

exit(EXIT_FAILURE);

}

laKopiere_(*x,yn);

no = laNorm_(’2’,yn); // sqrt(laSkalProd_(yn,yn));

la_n = 0;

i = 0;

do

{

la_a = la_n;


laKopiere_(yn,y);

gestaffelte_LGS(n,LU,y,&yn,perm);

i = i+1;

no = laNorm_(’2’,yn);

la_n = no;

} while((fabs(la_a-la_n)>eps)&&(i<maxiter));

for(j=0;j<n;j++)

if(y[j][0]*yn[j][0]!=0) break;

if(y[j][0]*yn[j][0]<0) la_n = -la_n;


/* Ergebnisse*/

laKopiere_(yn,*x);

*iter = i;

*la = 1.0/la_n;

laVernichte_(y);

laVernichte_(yn);

laVernichte_(LU);

laVernichte_(perm);

}

Die Funktionen fur die Ein- und Ausgabe sind problemangepasst und im Modul laInOut_.c

zusammengefasst. Sie werden alle genutzt und noch durch zwei neue zur Matrixeingabe und-ausgabe erganzt.Dort werden auch die Funktionen laMatrixNeu_, laDim_ aufgerufen.

// laInOut_.c

/* Eine Bibliothek fuer Lineare Algebra, Ein/Ausgaben */

/* Definitionen von laGlSysEinlesen_, laGlSysAusgeben_, laMatAusgeben_, laVektAusgeben_, laVektAusgebenI_ */

#include <stdlib.h>

#include <stdio.h>


#include "la_.h"

...

/*----------------------------------------------------------------------------*/

void laMatEinlesen_(laMatrix *A)

{

73

FILE *matIn;

int i,j,n,m;


printf("Dateiname fuer Daten zur Matrix A : ");

scanf("%s",&s); // auch scanf("%s",s);

matIn = fopen(s,"r"); // matIn = fopen("Matrix.in","r");

if(matIn==NULL)

{


getch();

exit(EXIT_FAILURE);

}



for(i=0;i<n;i++)

for(j=0;j<m;j++)


fclose(matIn);

}

/*----------------------------------------------------------------------------*/

void laMatAusgeben_(laMatrix A)

{

int i,j,n,m;

laDim_(A,&n,&m);

for(i=0;i<n;i++)

{

for(j=0;j<m;j++)

printf(" %10.6lf ",A[i][j]);

printf("\n");

}

}

Die Definitionen befinden sich in verschiedenen Quelltexten sowie im Modul la_.c.

Die Beispiele werden als ASCII-Dateien in der Datenreihenfolge n,m, A(1..n, 1..m), b(1..n), ...bereitgestellt. Der Dateiname fur die Daten A, b zum LGS dient auch fur die Bereitstellungder Matrix A selber. Als Testvarianten nehmen wir regulare und singulare Matrizen(reg.: matrix23, matrix3, matrix31, sing.: matrix58).

matrix23 matrix3 matrix31 matrix58

-------- ------------------- ---------------- --------

3 3 2 2 4 4 3 3

1 -2 2 1.441969 1.040807 1 5 3 9 2 2 2

-1 1 -1 1.040807 0.751250 5 26 16 47 1 3 0

-2 -2 1 0.401162 3 16 14 35 2 4 1

1 0.289557 9 47 35 95 5

-1 Loesung: -12 4

-3 1,-1 -62 7

Loesung: -38 Loesung:

1,1,1, -114 keine

Loesung:

2,-1,3,-2

EW: EW: EW: EW:

1 2.193218999999544 0.031058762099 0

1 4.5595081932092e-13 0.285539766776 1

1 3.409857734575 5

132.273543736549

74

Nun fehlt nur noch der Quelltext (Rahmenprogramm) meinprog_.c mit der main-Funktion.In ihm befinden sich die notwendigen include-Anweisungen, die genannten vier Funktionenzum partiellen Eigenwertproblem, Eingaben, Berechnungen und Ausgaben.

// meinprog_.c

#include<stdlib.h>

#include<stdio.h>

#include<math.h>

#include "la_.h"

void gauss_LU(int n, int *t, laMatrix *A, laMatrix *perm, int *sing)

...

void gestaffelte_LGS(int n, laMatrix LU, laMatrix b, laMatrix *x, laMatrix perm)

...

void vektorit(int n, laMatrix A, int maxiter, double eps, laMatrix *x, double *la, int *iter)

...

void inv_vektorit(int n, laMatrix A, int maxiter, double eps, laMatrix *x, double *la, int *iter)

...

int main()

{

int n,m,i,iter,maxiter;

double eps,la;

laMatrix A,x,x0;

char ein,ret,kont;

printf("Vektoriteration und inverse Vektoriteration fuer Eigenwerte einer Matrix\n");

printf("VI: y(m+1)=Ay(m), |la(A)|_max ~ ||y(m+1)||/||y(m)|| \n");

printf("inv. VI: Ay(m+1)=y(m), |la(A)|_min = 1/|la(A^-1)|_max \n");

printf("Loesung des LGS mittels LU-Faktorisierung und gestaffelter Systeme\n");

printf("-----------------------------------------------------------------------------\n\n");

laMatEinlesen_(&A);

fflush(stdin);

printf("Kontrollausgabe A (j/n)? ");

scanf("%c",&kont,&ret);

if(kont==’j’)

{

laDim_(A,&n,&m);

printf("\nDimension A(n,m): n=%i, m=%i\n",n,m);

printf("\nA\n");

laMatAusgeben_(A);

printf("------------------------------------------------------------------------------\n\n");

scanf("%c",&ret);

}

laDim_(A,&n,&m);

if(n!=m)

{

printf("\nDimension n<>m, keine Berechnung");

scanf("%c",&ret);

return (-1);

}

x = laMatrixNeu_(n,1);

x0 = laMatrixNeu_(n,1);

fflush(stdin);

printf("Eingabe/Belegung des Startvektors\n");

printf("Definition/Einsvektor/1.Einheitsvektor (d/e/i): ");

scanf("%c",&ein,&ret);

switch (ein)

{

case ’d’: for(i=0; i < n; i++)

{ printf("x0[%i] = ",i);

scanf("%lf",&x0[i][0]);

}

break;

case ’e’: for(i=0; i < n; i++) x0[i][0] = 1.0;

break;

default: x0 = laEins_(n,1);

break;

}

printf("Eingabe Toleranz eps>0: ");

scanf("%lf",&eps);

printf("Eingabe maximale Iterationsanzahl maxiter>0: ");

scanf("%i",&maxiter);

75

/* Berechnung */

laKopiere_(x0,x);

vektorit(n,A,maxiter,eps,&x,&la,&iter);

printf("\nVektoriteration\n");

printf("Betragsgroesster Eigenwert la = %19.16e\n",la);

printf("Zugehoeriger skalierter Eigenvektor x\n");

laVektAusgeben_(x);

printf("Iterationsanzahl iter = %i\n",iter);

printf("---------------------------------------------------------------------------------------\n");

laKopiere_(x0,x);

inv_vektorit(n,A,maxiter,eps,&x,&la,&iter);

printf("\nInverse Vektoriteration\n");

printf("Betragskleinster Eigenwert la = %19.16e\n",la);

printf("Zugehoeriger skalierter Eigenvektor x\n");

laVektAusgeben_(x);

printf("Iterationsanzahl iter = %i\n",iter);

printf("---------------------------------------------------------------------------------------\n");

fflush(stdin);

scanf("%c",&ret);

// Speicherplatz freigeben

laVernichte_(A);

laVernichte_(x);

}

Zusammenstellung des C-Projekts vektorit im Arbeitsverzeichnis

C:\d\Neundorf\nwptexte\tech phy\05\dev cpp\projekt5.

Im Projektfenster steht der Baum des Projekts. Die wirklich verwendeten Module sind fetthervorgehoben.

2- –ℵ⋃

≡ vektorit||− la .c||− laAdd_.c||− laEins .c||− laElemMult_.c||− laInOut .c||− laKopiere .c||− laKopiere1_.c||− laKreuzProd_.c||− laMult_.c||− laNorm .c||− laNull .c||− laSkalProd .c||− laSkMult .c||− laSub_.c||− laVertauscheSp_.c||− laVertauscheZe .c|− meinprog .c

Wie in den Projekten vorher kompiliert man das Projekt vektorit und fuhrt es aus.Es folgen Informationen zum Ubersetzen und im Kompilier Log.Im Verzeichnis stehen dann die Projektdatei vektorit.dev, die ausfuhrbare Dateivektorit.exe sowie alle Objektdateien *.o der Projektquellen.


76


A =

1 −2 2−1 1 −1−2 −2 1

, EW: 1, 1, 1

Der einzige Eigenwert λ ist betragsgroßter und betragskleinster sowie mehrfach.Zugehoriger einziger Eigenvektor ist x = (−1, 1, 1)T , denn der Rang der Matrix A−λI ist 2.Die Vektoriteration konvergiert, aber wegen der Vielfachheit sehr langsam.


2 C:\d\Neundorf\nwptexte\tech phy\05\dev cpp\projekt5\vektorit.exe

Vektoriteration und inverse Vektoriteration fuer Eigenwerte einer MatrixVI: y(m+1)=Ay(m), | la(A)| max ∼ ‖y(m+1)‖/‖y(m)‖Inv. VI: Ay(m+1)=y(m), | la(A)| min = 1/|la(A∧ − 1 | maxLoesung des LGS mittels LU-Faktorisierung und gestaffelter Systeme——————————————————————————————————–

Dateiname fuer Daten zur Matrix A : matrix23Kontrollausgabe A (j/n)? jDimension A(n,m): n=3, m=3A

1.000000 -2.000000 2.000000-1.000000 1.000000 -1.000000-2.000000 -2.000000 1.000000

——————————————————————————————————–

Eingabe/Belegung des StartvektorsDefinition/Einsvektor/1.Einheitsvektor (d/e/i): dx0[0] = 1x0[1] = 2x0[2] = 3Eingabe Toleranz eps>0: 1e-12Eingabe maximale Iterationsanzahl maxiter>0: 300

VektoriterationBetragsgroesster Eigenwert la = 1.0067340826915268e+000Zugehoeriger skalierter Eigenvektor x

-0.579284 0.577350 0.575410Iterationsanzahl iter = 300——————————————————————————————————–

Inverse VektoriterationBetragskleinster Eigenwert la = 9.9337741046664496e-001Zugehoeriger skalierter Eigenvektor x

-0.575435 0.577350 0.579259Iterationsanzahl iter = 300——————————————————————————————————–

Ahnlich ist es mit den Startvektoren als Einsvektor oder erster Einheitsvektor.Beginnt man mit dem Startvektor (−1, 1, 1)T (Eigenvektor), dann konvergieren beide Va-rianten nach 2 Schritten und es sind la=1.000000000000 und der skalierte Eigenvektor(−0.577350, 0.577350, 0.577350)T .

77

Beispiel 2 Ax = b, A regular und symmetrisch, det(A) = 10−12

A =

(

1.441969 1.040807

1.040807 0.751250

)

, EW: 2.193218999999544, 4.5595081932092e − 13




Dateiname fuer Daten zur Matrix A : matrix3Kontrollausgabe A (j/n)? nEingabe/Belegung des StartvektorsDefinition/Einsvektor/1.Einheitsvektor (d/e/i): dx0[0] = 1x0[1] = 2Eingabe Toleranz eps>0: 1e-15Eingabe maximale Iterationsanzahl maxiter>0: 10


0.810843 0.585263Iterationsanzahl iter = 3——————————————————————————————————–


-0.585263 0.810843Iterationsanzahl iter = 3——————————————————————————————————–

Ahnlich ist es mit den Startvektoren als Einsvektor oder erster Einheitsvektor.Auch bei noch kleineren Toleranzen erhalt man die gleiche Anzahl von Iterationen.Mit der double-Genauigkeit konnen bei der Berechnung des betragskleinsten Eigenwerts nurungefahr die ersten 4 Dezimalstellen garantiert werden.Noch eine Rechnung mit dem Startvektor (100,−100)T und sonst ohne Veranderungen.

...VektoriterationBetragsgroesster Eigenwert la = 2.1932189999995444e+000Zugehoeriger skalierter Eigenvektor x

0.810843 0.585263Iterationsanzahl iter = 3——————————————————————————————————–


0.585263 -0.810843Iterationsanzahl iter = 4——————————————————————————————————–

78

Beispiel 3 Ax = b, A symmetrisch positiv definit, damit regular

A =

1 5 3 95 26 16 473 16 14 359 47 35 95

,

EW: 0.031058762099, 0.285539766776, 3.409857734575, 132.273543736549

Die betragsgroßten und betragskleinsten Eigenwerte sind gut separiert von den anderen, sodass die Iterationen schnell konvergieren.




Dateiname fuer Daten zur Matrix A : matrix31Kontrollausgabe A (j/n)? jDimension A(n,m): n=4, m=4

A1.000000 5.000000 3.000000 9.0000005.000000 26.000000 16.000000 47.0000003.000000 16.000000 14.000000 35.0000009.000000 47.000000 35.000000 95.000000

——————————————————————————————————–

Eingabe/Belegung des StartvektorsDefinition/Einsvektor/1.Einheitsvektor (d/e/i): dx0[0] = 1x0[1] = 2x0[2] = 3x0[3] = 4Eingabe Toleranz eps>0: 1e-12Eingabe maximale Iterationsanzahl maxiter>0: 20


0.081311 0.424912 0.310068 0.846579Iterationsanzahl iter = 6——————————————————————————————————–


0.980914 -0.140818 0.116582 -0.066234Iterationsanzahl iter = 10——————————————————————————————————–

79

Weitere Berechnungen mit anderen Startvektoren



Dateiname fuer Daten zur Matrix A : matrix31Kontrollausgabe A (j/n)? nEingabe/Belegung des StartvektorsDefinition/Einsvektor/1.Einheitsvektor (d/e/i): eEingabe Toleranz eps>0: 1e-12Eingabe maximale Iterationsanzahl maxiter>0: 20





...Dateiname fuer Daten zur Matrix A : matrix31Kontrollausgabe A (j/n)? nEingabe/Belegung des StartvektorsDefinition/Einsvektor/1.Einheitsvektor (d/e/i): iEingabe Toleranz eps>0: 1e-12Eingabe maximale Iterationsanzahl maxiter>0: 20





80


A =

2 2 21 3 02 4 1

, EW: 0, 1, 5

Die inverse Vektoriteration fur den betragskleinsten Eigenwert ist nicht durchfuhrbar. Schonbei der LU -Faktorisierung der Matrix A wird ihre Singularitat festgestellt und signalisiert.




Dateiname fuer Daten zur Matrix A : matrix58Kontrollausgabe A (j/n)? jDimension A(n,m): n=3, m=3A

2.000000 2.000000 2.0000001.000000 3.000000 0.0000002.000000 4.000000 1.000000

——————————————————————————————————–

Eingabe/Belegung des StartvektorsDefinition/Einsvektor/1.Einheitsvektor (d/e/i): dx0[0] = 1x0[1] = 2x0[2] = 3Eingabe Toleranz eps>0: 1e-12Eingabe maximale Iterationsanzahl maxiter>0: 20


0.666667 0.333333 0.666667Iterationsanzahl iter = 19——————————————————————————————————–


Weitere Berechnungen

Startvektor (1, 1, 1)T , ε = 10−12: la = 5.0000000000002496e+000, iter = 19




Startvektor (1,−2, 3)T , ε = 10−12: la = 5.0000000000001625e+000, iter = 21

Startvektor (1,−2, 3)T , ε = 10−15: la = 5.0000000000000009e+000, iter = 25

81

5 Felder in C++ und C

In den Programmen matsum.cpp, matsum.c (Abschnitt 1.4) sowie den gerade genanntenProjekten gauss4c, gauss4ca, gauss4cb, gauss4cpp1, gauss4cpp2, vektorit zumGauß-Algorithmus / Vektoriteration und speziell in den zahlreichen Definitionen/Modulenwie z. B. laMatrixNeu_, laVernichte_, laGlSysEinlesen_, ist der Umgang mit Matrizenund Vektoren und darauf aufbauenden Objekten schon implementiert worden.Es soll jedoch noch einmal in ubersichtlicher Form die Deklaration, Definition und Verar-beitung von Feldern unter verschiedenen Aspekten dargestellt werden. Das betrifft nicht nurdie Unterschiede zwischen den Sprachdialekten C++ und C, sondern auch die Verwendungvon statischen und dynamischen Strukturen.

Zunachst noch einmal die Zuordnung von einigen Objekten und Standardbezeichnern zuHeader-Files.

C++// fuer Objekte wie

#include <iostream.h> // cin, cout,...

#include <conio.h> // getch, printf, scanf, stdin,...

#include <stdlib.h> // EXIT_FAILURE,...

#include <stdlib.h> // koennen sonst entfallen

#include <stdio.h>

#include <math.h>

C//#include <iostream.h> // muss entfallen

#include <math.h> // falls mathematische Standardfunktionen auftreten

//getch, printf, scanf, ... werden automatisch erkannt

// fuer Objekte wie


#include <stdio.h> // stdin, stderr, stdout,...

#include <conio.h> // stdin,... enthaelt Header-File stdio.h

5.1 Statische Felder

Statische bzw. nichtdynamische Felder (Array) sind hintereinander abgelegte Objekte glei-chen Datentyps. Hier muss die Anzahl der Dimensionen und Elemente in einem Feld beiseiner Deklaration (also zur Programmierzeit) angegeben werden.Die Anwendung des Indizierungsoperators auf den Feldnamen ist nur deshalb moglich, weilder Compiler den Feldnamen in Anweisungen als Zeiger auf das erste Element de Feldesansieht. Die Indizierung wird intern in die Zeigerberechnung umgewandelt.

Eindimensionales Feld: feld[i] ⇔ *(feld+i), i=0,1,2,...

Zeigerarithmetik bedeutet, dass der Wert eines Zeigers bei der Addition oder Subtraktionnicht in Einheiten von Bytes, sondern in Einheiten von der Große der Objekte, auf die derZeiger zeigt, erhoht oder vermindert wird.

Die gezeigten Programme arbeiten mit einem zweidimensionalen Feld. Es sind verschiedeneMoglichkeiten der Deklaration der Anzahl der Elemente in den beiden Dimensionen gezeigt.Von diesem “maximalen“ Feld wird in Abhangigkeit von der eingegebenen Zeilen- und Spal-tenzahl nur ein Teil genutzt. Wurden diese beiden Angaben (oder auch nur eine davon) dievereinbarte obere Grenze ubersteigen, kommt es i. Allg. irgendwann wahrend der Programm-abarbeitung zu einer Fehlermeldung, in C++ eher als in C.a1stat.exe hat Fehler verursacht und wird geschlossen. Starten Sie das Programm neu.

Ein Fehlerprotokoll wird erstellt.

82

C++

// a1stat.cpp

// Felder: statisch

// Berechnen der Zeilensummennorm einer (n,m)-Matrix

#include <iostream.h>

#include <conio.h>

//#define nmax 50

int main()

{

const int nmax=50; // bzw. int nmax=50;

double A[nmax][nmax],sum,max;

int n,m,i,j;

cout<<"Zeilensummennorm einer max. (50,50)-Matrix\n\n";

do

{

cout<<"Eingabe Zeilenzahl n (0<n<=50): ";

cin>>n;

} while ((n<=0)||(n>50));

do

{

cout<<"Eingabe Spaltenzahl m (0<m<=50): ";

cin>>m;

} while ((m<=0)||(m>50));

for(i=0;i<n;i++)

{

for(j=0;j<m;j++)

{

cout<<"A["<<i<<","<<j<<"] = ";

cin>>A[i][j];

}

}

max=0.0;

for(i=0;i<n;i++)

{

sum=0.0;

for(j=0;j<m;j++)

sum=sum+fabs(A[i][j]); // sum+=fabs(A[i][j]);

max=(sum>max?sum:max); // if (sum>max) max=sum;

}

cout<<"\nZeilensummennorm = "<<max;

getch();

return 0;

}

C

// a1stat.c

// Felder: statisch


//#include <iostream.h>

#include <conio.h>

#include <math.h>

//#define nmax 50

int main()

{

const int nmax=50; // bzw. int nmax=50;

double A[nmax][nmax],sum,max;

int n,m,i,j;

char ret;

printf("Zeilensummennorm einer max. (50,50)-Matrix\n\n");

do

{

printf("Eingabe Zeilenzahl n (0<n<=50): ");

83

scanf("%i",&n);

} while((n<=0)||(n>50));

do

{

printf("Eingabe Spaltenzahl m (0<m<=50): ");

scanf("%i",&m);

} while((m<=0)||(m>50));

for(i=0;i<n;i++)

{

for(j=0;j<m;j++)

{

printf("A[%i,%i] = ",i,j);

scanf("%lf",&A[i][j]);

}

}

max=0.0;

for(i=0;i<n;i++)

{

sum=0.0;

for(j=0;j<m;j++)



}

printf("\nZeilensummennorm = %lf",max);

// Tastaturpuffer leeren

fflush(stdin); // stdin in stdio.h bzw. conio.h

scanf("%c",&ret); // anstelle von getch()

return 0;

}

Eingabe- und Ergebnisfenster zu a1stat.c

2 C:\d\Neundorf\nwptexte\tech phy\05\cpp gpp gcc\a1stat.exe

Zeilensummennorm einer max. (50,50)-Matrix

Eingabe Zeilenzahl n (0<n<=50): 2Eingabe Spaltenzahl m (0<m<=50): 5A[0,0] = 1A[0,1] = 2A[0,2] = 3A[0,3] = -1A[0,4] = -2A[1,0] = 4A[1,1] = 5A[1,2] = 6A[1,3] = 0A[1,4] = -1

Zeilensummennorm = 16.000000

84

5.2 Dynamische Felder

Dynamische Felder (Array) sind ebenfalls Objekte gleichen Datentyps. Hier muss die Anzahlder Dimensionen und Elemente in einem Feld erst zur Laufzeit des Programms feststehenund kann sich daher den aktuellen Gegegebenheiten anpassen. Die dynamische Reservierungvon Feldern ist mit der Verwendung von Zeigern verbunden. Die Anwendung des Indizie-rungsoperators auf den Feldnamen ist deshalb moglich, weil der Compiler den Feldnamen inAnweisungen als Zeiger auf das erste Element de Feldes ansieht. Die Indizierung und Zeiger-berechnung sind ineinander umwandelbar.

Die dynamische Speicherallokation mittels den Zeiger reserviert den gewunschten Speicher-platz, der dann ublicherweise auch wieder freizugeben ist, bevor der Zeiger nicht mehr ge-nutzt, ungultig oder umgelenkt wird (Vermeidung von so genannten Speicherlecks).

Speicherallokation Speicherfreigabe

-----------------------------------------------------

C++ new delete

C++/C malloc() free()

In der Zeigerprogrammierung sind u. a. wichtig

– die Deklaration von Zeigern mit dem bestimmten Datentyp,

– generische Zeiger (void) ohne Typbezogenheit,

– explizite Typumwandlung bei Zeigern,

– Zeiger und Felder,

– Zeigerarithmetik und

– Zeiger auf Zeiger.

Die Typumwandlung bei Zeigern ist in den Projekten gauss4c,..., gauss4cpp2 zumGauß-Algorithmus und speziell in der Funktion laMatrixNeu_ schon mehrfach verwendetworden.

Am Beispiel eines zweidimensionalen Felds (rechteckige Matrix) sollen einige Aspekte de-monstriert werden.

Halbdynamische Felder

Zunachst soll speicherplatzsparend dadurch gearbeitet werden, dass bei fester maximalerMatrixzeile mit spmax Elementen (=Spalten) die Anzahl der Zeilen in Abhangigkeit von deraktuell eingegebenen Große n festgelegt wird. Fur die Matrixelemente braucht man dannn*spmax Speicherplatz von gegebenem Elementtyp (hier double). Dazu kommt noch Spei-cherplatz fur etwas Adressverwaltung.

C++Das ist die kurze und elegante Version unter Verwendung eines Zeigervektors (Vektor vonZeigern).Damit liegen die Matrixzeilen im Speicher genau hintereinander und es ist A[i] = *(A+i).Die Matrixelemente konnen mittels A[i][j] oder *(*(A+i)+j) aufgerufen werden.

// a1hdyn1.cpp

// Felder: halbdynamisch, Zeigervektor, new-delete



#include <conio.h>

85

//#define spmax 50

int main()

{

const int spmax=50; // maximale Spaltenzahl

typedef double vektor[spmax]; // Typ des Zeilenvektors

int i,j,n,m;

double max,sum;

cout<<"Zeilensummennorm einer (n,m)-Matrix, m<=50\n\n";

do

{

cout<<"Eingabe Zeilenzahl n (0<n): ";

cin>>n;

} while(n<=0);

do

{

cout<<"Eingabe Spaltenzahl m (0<m<=50): ";

cin>>m;

} while((m<=0)||(m>50));

// Allokation der Matrix auf dem Heap, damit A[i] = *(A+i)

vektor *A = new vektor[n]; // n beliebig

for(i=0;i<n;i++)

for(j=0;j<m;j++)

{

cout<<"A["<<i<<","<<j<<"]= ";

cin>>A[i][j];

}

max=0.0;

for(i=0;i<n;i++)

{

sum=0.0;

for(j=0;j<m;j++)



}


delete []A;

getch();

return 0;

}

CVerwendung des Zeigervektors.Ein Feld von Zeigern (Zeigerfeld) auf den Datentyp T wird als T *ptr[anzahl] deklariert.Dies ist nicht mit einem Zeiger auf ein Feld vom Typ T gemaß T (*ptr)[anzahl] zu ver-wechseln.Da die Deklaration des Zeigerfelds mit einer festen Feldgrenze erfolgen muss, gibt man dorteine maximale obere Schranke an. Die eingegebene Anzahl der Zeilen sollte dann nicht großerals diese Schranke sein. Die Dynamik der Allokation von Speicherplatz fur die Matrixzeilenbesteht dann darin, dass man nur Platz fur die gewunschten Zeilen allokiert. Ist die eingege-bene Zeilenzahl großer als die Schranke, wird es i. Allg. wahrend der Programmabarbeitungzur Fehlermeldung und damit zum Abbruch kommen.

Programmfehler

a1stat.exe hat Fehler verursacht und wird geschlossen. Starten Sie das Programm neu.

Ein Fehlerprotokoll wird erstellt.

Damit ist die dynamische Speicherverwaltung bezuglich der Zeilen nicht vollstandig machbar.

86

// a1hdyn1.c

// Felder: halbdynamisch mit Kontrolle, Zeigervektor, malloc-free


#include <conio.h>

#include <stdlib.h>

#include <math.h>

//#define spmax 50

//#define zemax 50

int main()

{


const int zemax=50; // maximale Zeilenzahl


double *A[zemax]; // Feld (Vektor) von Zeigern

int i,j,n,m;

double max,sum;

char ret;

printf("Zeilensummennorm einer (n,m)-Matrix, n,m<=50\n\n");

do

{

printf("Eingabe Zeilenzahl n (0<n<=50): ");

scanf("%i",&n);

} while((n<=0)||(n>50));

do

{


scanf("%i",&m);

} while(m<=0||m>50);

// Allokation der Matrix auf dem Heap

for(i=0;i<n;i++)

//if((*(A+i) = (double*) malloc(sizeof(vektor)))==NULL) // analog

if((A[i] = (double*) malloc(sizeof(vektor)))==NULL)

{

printf("Zu wenig Speicher\n");

scanf("%c",&ret);

exit(EXIT_FAILURE);

}

for(i=0;i<n;i++)

for(j=0;j<m;j++)

{



}

max=0.0;

for(i=0;i<n;i++)

{

sum=0.0;

for(j=0;j<m;j++)


// sum=sum+fabs(*(*(A+i)+j));


}


fflush(stdin);

for(i=0;i<n;i++) free(*(A+i));

scanf("%c",&ret);

return 0;

}

Nun soll ein Zeiger auf Zeiger genutzt werden.Dabei erkennt man, was notwendig ware, um auch in der Spaltendimension dynamisch zuwerden. Man braucht die Reservierung der Matrixzeilen nur mit der eingelesenen Lange m

anstelle von spmax zu machen.

87

// a1hdyn2.c

// Felder: halbdynamisch mit Kontrolle, Zeiger auf Zeiger, malloc-free


#include <conio.h>

#include <stdlib.h>

#include <math.h>

//#define spmax 50

int main()

{



double **A;

int i,j,n,m;

double max,sum;

char ret;

printf("Zeilensummennorm einer (n,m)-Matrix, m<=50\n\n");

do

{

printf("Eingabe Zeilenzahl n (0<n): ");

scanf("%i",&n);

} while (n<=0);

do

{


scanf("%i",&m);

} while ((m<=0)||(m>50));


if((A = (double**) malloc(n*sizeof(double*)))==NULL)

{


scanf("%c",&ret);

exit(EXIT_FAILURE);

}

for(i=0;i<n;i++)

{

// if((*(A+i) = (double*) malloc(spmax*sizeof(double)))==NULL) // analog

if((*(A+i) = (double*) malloc(sizeof(vektor)))==NULL)

{


scanf("%c",&ret);

exit(EXIT_FAILURE);

}

}

for(i=0;i<n;i++)

for(j=0;j<m;j++)

{



}

max=0.0;

for(i=0;i<n;i++)

{

sum=0.0;

for(j=0;j<m;j++)



}


fflush(stdin);


free(A);

scanf("%c",&ret);

return 0;

}

88

Dynamische Felder

Verwendung des Konzepts Zeiger auf Zeiger fur beide Sprachdialekte.

C++4 Programmversionen

// a1dyn1.cpp

// Felder: dynamisch, Zeiger auf Zeiger, new-delete



#include <conio.h>

int main()

{

typedef double **matrix; // Typ der Matrix

matrix A;

// double **A; // einfache Schreibweise

int i,j,n,m;

double max,sum;

cout<<"Zeilensummennorm einer (n,m)-Matrix\n\n";

do

{


cin>>n;

} while(n<=0);

do

{

cout<<"Eingabe Spaltenzahl m (0<m): ";

cin>>m;

} while(m<=0);


// mit Umtypisierung // einfache Variante a1dyn2.cpp

(char*)A = new char[n*sizeof(double*)]; // A = new double*[n];

for(i=0;i<n;i++) *(A+i) = new double[m];

for(i=0;i<n;i++)

for(j=0;j<m;j++)

{

cout<<"A["<<i<<","<<j<<"]= ";

cin>>A[i][j];

}

max=0.0;

for(i=0;i<n;i++)

{

sum=0.0;

for(j=0;j<m;j++)



}


for(i=0;i<n;i++) delete [](*(A+i));

// einfache Variante a1dyn2.cpp

delete [](char*)A; // delete []A;

getch();

return 0;

}

/*-------------------------------------------------------------------------------------*/

89

// a1dyn3.cpp

// Felder: dynamisch mit Kontrolle, Zeiger auf Zeiger, new-delete



#include <conio.h>

#include <stdlib.h>

int main()

{


matrix A;


int i,j,n,m;

double max,sum;

cout<<"Zeilensummennorm einer (n,m)-Matrix\n\n";

do

{


cin>>n;

} while(n<=0);

do

{

cout<<"Eingabe Spaltenzahl m (0<m): ";

cin>>m;

} while(m<=0);


// Fehlermeldung?

// ev. "Nicht genuegend virtueller Speicher,

// Auslagerungsdatei vergroessern,..."

if ((A = new double*[n])==NULL)

{


getch();

exit(EXIT_FAILURE);

}

for(i=0;i<n;i++)

if((*(A+i) = new double[m])==NULL)

{


getch();

exit(EXIT_FAILURE);

}

for(i=0;i<n;i++)

for(j=0;j<m;j++)

{

cout<<"A["<<i<<","<<j<<"]= ";

cin>>A[i][j];

}

max=0.0;

for(i=0;i<n;i++)

{

sum=0.0;

for(j=0;j<m;j++)



}


for(i=0;i<n;i++) delete [](*(A+i));

delete [](char*)A;

getch();

return 0;

}

/*-------------------------------------------------------------------------------------*/

90

// a1dyn4.cpp

// Felder: dynamisch mit Kontrolle, Zeiger auf Zeiger, malloc-free


//#include <iostream.h>

#include <conio.h>

#include <stdlib.h>

int main()

{


matrix A;


int i,j,n,m;

double max,sum;

char ret;

printf("Zeilensummennorm einer (n,m)-Matrix\n\n");

do

{


scanf("%i",&n);

} while(n<=0);

do

{

printf("Eingabe Spaltenzahl m (0<m): ");

scanf("%i",&m);

} while(m<=0);


if(((void*)A = malloc(n*sizeof(double*)))==NULL)

{


scanf("%c",&ret);

exit(EXIT_FAILURE);

}

for(i=0;i<n;i++)

if(((void*)(*(A+i)) = malloc(m*sizeof(double)))==NULL)

{


scanf("%c",&ret);

exit(EXIT_FAILURE);

}

for(i=0;i<n;i++)

for(j=0;j<m;j++)

{



}

max=0.0;

for(i=0;i<n;i++)

{

sum=0.0;

for(j=0;j<m;j++)



}


fflush(stdin);


free(A);

scanf("%c",&ret);

return 0;

}

Die Analogie zur Version a1dyn4.cpp ohne den generischen Zeiger void* gibt es auch in C.

91

C

// a1dyn4.c

// Felder: dynamisch mit Kontrolle, Zeiger auf Zeiger, malloc-free


#include <conio.h>

#include <math.h>

#include <stdlib.h>

int main()

{


matrix A;


int i,j,n,m;

double max,sum;

char ret;

printf("Zeilensummennorm einer (n,m)-Matrix\n\n");

do

{


scanf("%i",&n);

} while(n<=0);

do

{

printf("Eingabe Spaltenzahl m (0<m): ");

scanf("%i",&m);

} while(m<=0);


// Fehlermeldung?

// ev. "Nicht genuegend virtueller Speicher,

// Auslagerungsdatei vergroessern,..."

if((A = malloc(n*sizeof(double*)))==NULL)

{


scanf("%c",&ret);

exit(EXIT_FAILURE);

}

for(i=0;i<n;i++)

if((*(A+i) = malloc(m*sizeof(double)))==NULL)

{


scanf("%c",&ret);

exit(EXIT_FAILURE);

}

for(i=0;i<n;i++)

for(j=0;j<m;j++)

{



}

max=0.0;

for(i=0;i<n;i++)

{

sum=0.0;

for(j=0;j<m;j++)



}


fflush(stdin);


free(A);

scanf("%c",&ret);

return 0;

}

92

AnhangProgrammbeispiele in C++ und C

Zunachst noch einmal die Zuordnung von einigen Objekten und Standardbezeichnern zuHeader-Files.

C++// fuer Objekte wie

#include <iostream.h> // cin, cout,...

#include <conio.h> // getch, printf, scanf, stdin,...


#include <stdlib.h> // koennen sonst entfallen

#include <stdio.h>

#include <math.h>

C

//#include <iostream.h> // muss entfallen

#include <math.h> // falls mathematische Standardfunktionen auftreten

//getch, printf, scanf, ... werden automatisch erkannt

// fuer Objekte wie


#include <stdio.h> // stdin, stderr, stdout,...

#include <conio.h> // stdin,... enthaelt Header-File stdio.h

In den C++-Programmen verwenden wir nicht die Objekte cin, cout, sondern wie in Cdie Befehle scanf, printf.Neben der zum Teil abweichenden Einbindung von Header-Files in den beiden Sprachdialek-ten, konnen noch andere kleine Unterschiede in der Syntax auftreten. So verlangt z. B. C++bei der Definition von Konstanten zusatzlich die Typangabe. Daruber hinaus gibt es i. Allg.Ubereinstimmung zwischen den Versionen.

(1) Quadratwurzel

// wurzel1.c

// Quadratwurzelberechnung sqrt(x)

#include <stdio.h>

#include <math.h>

int main()

{

int i,n;

double xi,x;

n = 19;

printf("Quadratwurzelberechnung\n");

printf("--------------------------------------------------------------\n");

printf(" i sqrt(i) xi=i,sqrt(xi) x double,sqrt(x)\n");

for(i=1,x=1.0;i<=n;i++,x=x+1.0)

{

xi = i;

printf("%2i %19.16lf %19.16lf %19.16lf\n",i,sqrt(i),sqrt(xi),sqrt(x));

}

getch();

return 0;

}

-------------------------------------------------------------------------------------

// wurzel1.cpp

// Quadratwurzelberechnung sqrt(x)

#include <conio.h>

#include <math.h>

93

int main()

{

int i,n;

double xi,x;

n = 19;

printf("Quadratwurzelberechnung\n");

printf("--------------------------------------------------------------\n");

printf(" i sqrt(i) xi=i,sqrt(xi) x double,sqrt(x)\n");

for(i=1,x=1.0;i<=n;i++,x=x+1.0)

{

xi = i;

printf("%2i %19.16lf %19.16lf %19.16lf\n",i,sqrt(i),sqrt(xi),sqrt(x));

}

getch();

return 0;

}

2 C:\d\Neundorf\nwptexte\tech phy\05\cpp gpp gcc\wurzel1.exe

Quadratwurzelberechnung—————————————————————————————i sqrt(i) xi=i,sqrt(xi) x double,sqrt(x)1 1.0000000000000000 1.0000000000000000 1.00000000000000002 1.4142135623730951 1.4142135623730951 1.41421356237309513 1.7320508075688772 1.7320508075688772 1.73205080756887724 2.0000000000000000 2.0000000000000000 2.00000000000000005 2.2360679774997898 2.2360679774997898 2.23606797749978986 2.4494897427831779 2.4494897427831779 2.44948974278317797 2.6457513110645907 2.6457513110645907 2.64575131106459078 2.8284271247461903 2.8284271247461903 2.82842712474619039 3.0000000000000000 3.0000000000000000 3.0000000000000000

10 3.1622776601683795 3.1622776601683795 3.1622776601683795...

(2) Kreiszahl π

// pi1.c

// Kreiszahl Pi durch Integration des Einheitskreises

#include <conio.h>

#include <math.h>

int main()

{

int i,n;

double x,delta_x,pi,Rl,y;

const double Pi = 3.1415926535897932; // exakt als Vergleichswert

n = 1000;

delta_x = 1.0/n;

/* Flaeche des Viertelkreises */

Rl = 0.0;

for(i=0;i<n;i++)

{

x = i*delta_x;

y = sqrt(1.0-x*x);

Rl = Rl+y*delta_x;

}

pi = 4*Rl;

printf("Kreiszahl Pi durch Integration des Einheitskreises: Rechteckregel links\n\n");

printf("Berechneter Wert (Quadratur) Pi = %17.15lf\n",pi);

printf("Exakter Wert (16 Stellen) Pi = %17.15lf\n",Pi);

// printf("C/C++ : Pi = 17.15%lf",M_PI); // M_PI nicht da

getch();

}

94

----------------------------------------------------------------------------------------------

// pi1.cpp, Kreiszahl Pi durch 1.Quadranten des Einheitskreises

#include <conio.h>

#include <stdlib.h> // kann entfallen

#include <math.h> // kann entfallen

//#define M_PI 3.14159265358979323846

int main()

{

int i,n;

double x,delta_x,pi,Rl,y;


...

}

2 C:\d\Neundorf\nwptexte\tech phy\05\cpp gpp gcc\pi1.exe

Kreiszahl Pi durch Integration des Einheitskreises: Rechteckregel links

Berechneter Wert (Quadratur) Pi = 3.143555466911028Exakter Wert (16 Stellen) Pi = 3.141592653589793

// pi2.c

// Kreiszahl Pi nach Archimedes mittels Umfang von im Kreis einbeschriebenen regelmaessigen Vielecken

#include <conio.h>

#include <math.h>

int main()

{

int j;

double s1,s2,p1,p2,ix;


printf("Pi nach Archimedes\n");

printf(" j s1 Pi1 s2 Pi2\n");

printf("--------------------------------------------------------------------------- \n");

// Initialisierung mit Sechseck

s1 = 1; s2 = 1;

p1 = 3; p2 = 3;

ix = 3;

printf("%2i %17.15lf %17.15lf %17.15lf %17.15lf\n",1,s1,p1,s2,p2);

fflush(stdin);

for(j=2;j<40;j++)

{

ix = 2*ix;

s1 = sqrt(2-sqrt(4-s1*s1));

s2 = s2/sqrt(2+sqrt(4-s2*s2));

p1 = ix*s1;

p2 = ix*s2;

printf("%2i %17.15lf %17.15lf %17.15lf %17.15lf\n",j,s1,p1,s2,p2);

if((j==20)||(j==40)) getch();

}

printf("\nExakter Wert (16 Stellen) Pi = %17.15lf\n",Pi);

//printf("C/C++ : Pi = 17.15%lf",M_PI);

getch();

}

------------------------------------------------------------------------------------------------------

// pi2.cpp

// Kreiszahl Pi nach Archimedes mittels Umfang von im Kreis einbeschriebenen regelmaesssigen Vielecken

#include <conio.h>

#include <math.h>

//#define M_PI 3.14159265358979323846

int main()

{

int j;

double s1,s2,p1,p2,ix;


...

}

95

2 C:\d\Neundorf\nwptexte\tech phy\05\cpp gpp gcc\pi2.exe

Pi nach Archimedesj s1 Pi1 s2 Pi2

————————————————————————————————————-1 1.000000000000000 3.000000000000000 1.000000000000000 3.0000000000000002 0.517638090205042 3.105828541230250 0.517638090205041 3.1058285412302493 0.261052384440103 3.132628613281237 0.261052384440103 3.1326286132812384 0.130806258460286 3.139350203046872 0.130806258460286 3.1393502030468675 0.065438165643553 3.141031950890530 0.065438165643552 3.1410319508905096 0.032723463252972 3.141452472285344 0.032723463252974 3.1414524722854627 0.016362279207873 3.141557607911622 0.016362279207874 3.1415576079118588 0.008181208052471 3.141583892148936 0.008181208052470 3.1415838921483189 0.004090612582340 3.141590463236762 0.004090612582328 3.141590463228050

10 0.002045307360705 3.141592106043048 0.002045307360677 3.14159210599927111 0.001022653813994 3.141592516588155 0.001022653814027 3.14159251669215612 0.000511326923607 3.141592618640789 0.000511326923725 3.14159261936538313 0.000255663463975 3.141592645321216 0.000255663463951 3.14159264503369014 0.000127831731987 3.141592645321216 0.000127831732237 3.14159265145076615 0.000063915865994 3.141592645321216 0.000063915866151 3.14159265305503616 0.000031957932997 3.141592645321216 0.000031957933080 3.14159265345610317 0.000015978971709 3.141593669849427 0.000015978966540 3.14159265355637018 0.000007989482381 3.141592303811738 0.000007989483270 3.14159265358143719 0.000003994762034 3.141608696224804 0.000003994741635 3.14159265358770320 0.000001997367121 3.141586839655041 0.000001997370818 3.14159265358927021 0.000000998711352 3.141674265021758 0.000000998685409 3.14159265358966222 0.000000499355676 3.141674265021758 0.000000499342704 3.14159265358975923 0.000000249788979 3.143072740170040 0.000000249671352 3.14159265358978424 0.000000125559416 3.159806164941135 0.000000124835676 3.14159265358979125 0.000000063220273 3.181980515339464 0.000000062417838 3.14159265358979226 0.000000033320009 3.354101966249685 0.000000031208919 3.14159265358979327 0.000000021073424 4.242640687119286 0.000000015604460 3.14159265358979328 0.000000014901161 6.000000000000000 0.000000007802230 3.14159265358979329 0.000000000000000 0.000000000000000 0.000000003901115 3.14159265358979330 0.000000000000000 0.000000000000000 0.000000001950557 3.14159265358979331 0.000000000000000 0.000000000000000 0.000000000975279 3.14159265358979332 0.000000000000000 0.000000000000000 0.000000000487639 3.14159265358979333 0.000000000000000 0.000000000000000 0.000000000243820 3.14159265358979334 0.000000000000000 0.000000000000000 0.000000000121910 3.14159265358979335 0.000000000000000 0.000000000000000 0.000000000060955 3.14159265358979336 0.000000000000000 0.000000000000000 0.000000000030477 3.14159265358979337 0.000000000000000 0.000000000000000 0.000000000015239 3.14159265358979338 0.000000000000000 0.000000000000000 0.000000000007619 3.14159265358979339 0.000000000000000 0.000000000000000 0.000000000003810 3.14159265358979340 0.000000000000000 0.000000000000000 0.000000000001905 3.14159265358979341 0.000000000000000 0.000000000000000 0.000000000000952 3.141592653589793

Exakter Wert (16 Stellen) Pi = 3.141592653589793

96

(3) Primzahlen

// prim1.c

// Die ersten n Primzahlen

#include <math.h>

int main()

{

const n = 100; // auch int const n = 100;

int i,j,k,prim;

double is;

i = 1;

k = 0;

printf("Die ersten n Primzahlen\n");

do

{

i++;

is = sqrt(i);

prim = 1;

for(j=2;j<=is;j++)

if((i%j)==0) // Modulo = Rest der Division

{

prim = 0;

break;

}

if(prim!=0)

{

k++;

printf("%3i.te Primzahl = %i\n",k,i);

}

} while(k<n);

getch();

}

-----------------------------------------------------------------------------------------

// prim1.cpp

// Die ersten n Primzahlen

#include <conio.h>

int main()

{

int const n = 100; // auch const int n = 100;

int i,j,k,prim;

double is;

...

}

2 C:\d\Neundorf\nwptexte\tech phy\05\cpp gpp gcc\prim1.exe

Die ersten n Primzahlen1. Primzahl = 22. Primzahl = 33. Primzahl = 54. Primzahl = 75. Primzahl = 11

...10. Primzahl = 2911. Primzahl = 31...26. Primzahl = 101...99. Primzahl = 523

100. Primzahl = 541

97

// prim2.c

// Die ersten n Primzahlen als Vektor

#include <math.h>

int main()

{

const n = 100;

int i,j,k,prim,vprim[n];

i = 2;

k = 0;

printf("Die ersten n Primzahlen als Vektor\n");

vprim[0] = 2;

printf("%3i. Primzahl = %i\n",k+1,i);

do

{

i++;

prim = 1;

for(j=0;j<k;j++)

if((i%vprim[j])==0)

{

prim = 0;

break;

}

if(prim!=0)

{

++k;

vprim[k] = i;

printf("%3i. Primzahl = %i\n",k+1,i);

}

} while(k<n-1);

getch();

printf("\nAusgabe des Vektors der Primzahlen\n");

for(k=0;k<n;k++) printf("vprim[%3i] = %i\n",k,vprim[k]);

getch();

}

---------------------------------------------------------------------------

// prim3.c

// Die ersten n Primzahlen als Vektor

#include <math.h>

int main()

{

const n = 10;

int i,j,k,prim,is,vprim[n+1];

i = vprim[0] = 1;

k = 2;

vprim[1] = 2;

vprim[2] = 3;

printf("Die ersten n Primzahlen als Vektor\n");

printf("%3i. Primzahl = %i\n",1,vprim[1]);

do

{

i = i+2;

is = (int)sqrt(i); // Typumwandlung

prim = 1;

j = 2;

while((vprim[j]<=is)&&(prim!=0))

{

prim = i%vprim[j];

j++;

}

if(prim!=0) printf("%3i. Primzahl = %i\n",k++,vprim[k]=i);

} while(k<=n);

getch();

printf("\nAusgabe des Vektors der Primzahlen\n");

for(k=1;k<=n;k++) printf("vprim[%3i] = %i\n",k,vprim[k]);

getch();

}

98


Die ersten n Primzahlen1. Primzahl = 22. Primzahl = 33. Primzahl = 5

...10. Primzahl = 29

Ausgabe des Vektors der Primzahlenvprim[ 1] = 2vprim[ 2] = 3vprim[ 3] = 5...vprim[ 10] = 29

// prim4.c

// Alle Primzahlen <= n, Sieb des Eratosthenes

#include <math.h>

int main()

{

const nmax = 1000;

int i,k,n,ns,vprim[nmax+1];

printf("Sieb des Eratosthenes: alle Primzahlen <= n\n");

do

{

printf("n (2<=n<=%i) = ",nmax);

scanf("%i",&n);

} while((n<2)||(n>nmax));

vprim[0] = 0;

vprim[1] = 1;

for(i=2;i<=n;i++) vprim[i]=i;

ns = (int)sqrt(n);

for(i=2;i<=ns;i++)

if(vprim[i]!=0)

{

k = 2*i;

while(k<=n)

{

vprim[k] = 0;

k = k+i;

}

}

printf("\nPrimzahlen <= %i\n",n);

for(i=2;i<=n;i++)

{

if(vprim[i]!=0) printf(" %i",i);

if(i%30==0) getch();

}

getch();

}


Sieb des Eratosthenes: alle Primzahlen <= nn (2<=n<=1000) = 100

Primzahlen <= 1002 3 5 7 11 13 17 19 31 37 41 43 47 53 59 61 71 7379 83 89 97

99

// prim5.c

// Primzahltest mittels ganzzahliger Division

#include <math.h>

int main()

{

int n,ns,teiler,kleinster_teiler,prim;

printf("Primzahltest mittels ganzzahliger Division\n");

do

{

printf("Zahl n (n>=2) = ");

scanf("%i",&n);

} while(n<=1);

prim = 1;

teiler = 2;

kleinster_teiler = 2;

ns = (int)sqrt(n);

while(teiler<=ns)

{

if((n%teiler==0)&&(n>3))

{

prim = 0;

kleinster_teiler = teiler;

break;

}

teiler++;

}

if(prim)

printf("\n%i ist Primzahl",n);

else

{

printf("\n%i ist keine Primzahl",n);

printf("\nKleinster Teiler = %i",kleinster_teiler);

}

getch();

}

----------------------------------------------------------------------------------------------

// prim5.cpp

// Primzahltest mittels ganzzahliger Division

#include <conio.h>

int main()

{

int n,ns,teiler,kleinster_teiler,prim;

...

}


Primzahltest mittels ganzzahliger DivisionZahl n (n>=2) = 93

93 ist keine PrimzahlKleinster Teiler = 3—————————————————————————————————————

Primzahltest mittels ganzzahliger DivisionZahl n (n>=2) = 107

107 ist Primzahl

100

(4) Summenberechnungen

// summe1_gauss.c

// Summe 1+2+...+n=n(n+1)/2 Gauss-Formel

main()

{

int i,n,s,g;

n = 100;

s = 0;

for(i=1;i<=n;i=i+1)

s = s+i; // s+=i;

g = n*(n+1)/2; // Ergebnis ist ganzzahlig

printf("Gauss-Formel 1+2+...+n=n(n+1)/2\n");

printf("n=%i\n",n);

printf("g=n(n+1)/2 = %i\n",g);

printf("s=0+1+2+...+n = %i\n",s);

getch();

}

----------------------------------------------------------------------

// summe1_gauss.cpp

// Summe 1+2+...+n=n(n+1)/2 Gauss-Formel

#include <conio.h>

main()

{

int i,n,s,g;

...

}

2 C:\d\Neundorf\nwptexte\tech phy\05\cpp gpp gcc\summe1 gauss.exe

Gauss-Formel 1+2+...+n=n(n+1)/2n=100g=n(n+1)/2 = 5050s=0+1+2+...+n = 5050

// summe2_geom.c

// Summe 2=1+1/2+1/2^2+1/2^3+1/2^4+.... geometrische Reihe

int main()

{

int i,n;

double s,a;

n = 16;

a = 1.0;

s = 1.0;

for(i=1;i<=n;i=i+1)

{

a = a/2.0;

s = s+a;

}

printf("Geometrische Reihe 2=1+1/2+1/2^2+1/2^3+1/2^4+....\n");

printf("Partialsumme s(n)=s(%i)=%14.12lf\n",n,s);

getch();

return 0;

}

2 C:\d\Neundorf\nwptexte\tech phy\05\cpp gpp gcc\summe2 geom.exe

Geometrische Reihe 2=1+1/2+1/2∧2+1/2∧3+1/2∧4+....Partialsumme s(n)=s(16)=1.999984741211

101

// summe3_sin.c

// Summe sin(x)=x-x^3/3!+x^5/5!-x^7/7!+x^9/9!-... Potenzreihe

#include <math.h>

int main()

{

int i,n;

double x,x2,s,a;

x = 1.5;

x2 = x*x;

n = 15;

a = x;

s = x;

for(i=2;i<=n;i=i+2)

{

a = -a*x2/(i*(i+1));

s = s+a;

}

printf("Potenzreihe sin(x)=x-x^3/3!+x^5/5!-x^7/7!+x^9/9!-...\n");

printf("sin(%lf)=%17.15lf\n",x,sin(x));

printf("Partialsumme s(n)=s(%i)=%17.15lf\n",n,s);

getch();

return 0;

}

2 C:\d\Neundorf\nwptexte\tech phy\05\cpp gpp gcc\summe3 sin.exe

Potenzreihe sin(x)=x-x∧3/3!+x∧5/5!-x∧7/7!+x∧9/9!-...sin(1.500000)=0.997494986604054Partialsumme s(n)=s(15)=0.997494986601303

// summe4_sin.c


// nach Transformation von x auf das Intervall [-Pi,Pi]

#include <math.h>

int main()

{

const double Pi = 3.1415926535897932;

int i,n,Q,h,vz;

double xe,x,x2,sinx,xt,s,a;


printf("nach Transformation von x auf das Intervall [-Pi,Pi]\n\n");

n = 100;

vz = 1;

printf("x = ");

scanf("%lf",&xe);

/*Transformation in den Bereich von -Pi bis +Pi*/

x = xe;

if(x<0)

{

x = -x;

vz = -1;

}

Q = (int)(x/(0.5*Pi));

h = Q%2; // Rest der ganzzahligen Division

xt = (x-2.0*Pi*(Q/4+(Q/2)%2))*(1-2*h*h)+Pi*h;

if(xt>Pi) xt = xt-2*Pi;

/*Berechnung von sin(x)*/

x2 =xt*xt;

for(i=0,sinx=0.0,s=a=xt;(sinx!=s)&&(i<n);i=i+1)

s=(sinx=s)+(a=-a*x2/((2*i+2)*(2*i+3)));

sinx = vz*sinx;

102

printf("\nErgebnisse\n");

printf("x =%19.15lf\n",xe);

printf("xt=%19.15lf\n\n",xt);

printf("sin(%19.15lf)=%19.15lf\n",xe,sin(xe));

printf("Partialsumme s(n)=s(%i)=%19.15lf\n",n,sinx);

getch();

return 0;

}


Potenzreihe sin(x)=x-x∧3/3!+x∧5/5!-x∧7/7!+x∧9/9!-...nach Transformation von x auf das Intervall [-Pi,Pi]

x = 0.5

Ergebnissex = 0.500000000000000xt = 0.500000000000000

sin( 0.500000000000000)= 0.479425538604203Partialsumme s(n)=s(100)= 0.479425538604203

————————————————————————————————————-

x = 2.2

Ergebnissex = 2.200000000000000xt = 0.941592653589793


————————————————————————————————————-

x = 3.5

Ergebnissex = 3.500000000000000xt = -2.783185307179586

sin( 3.500000000000000)= -0.350783227689620Partialsumme s(n)=s(100)= -0.350783227689620

————————————————————————————————————-

x = 6

Ergebnissex = 6.000000000000000xt = -2.858407346410207


————————————————————————————————————-

x = 100

Ergebnissex =100.000000000000000xt = -2.610627738716413


103

// summe5_sin.c


// nach Transformation von x auf das Intervall [0,Pi/2]

#include <math.h>

int main()

{

const double Pi = 3.1415926535897932;

int i,n,vz,ix,vorz;

double xe,x,x2,sinx,s,a;


printf("nach Transformation von x auf das Intervall [0,Pi/2]\n\n");

n = 100;

vz = 1;

vorz = 1;

printf("x = ");

scanf("%lf",&xe);

/*Transformation in den Bereich von 0 bis 2*Pi*/

x = xe;

if(x<0)

{

x = -x;

vz = -1;

}

ix = (int)(x/(2*Pi)); // wie oft ist 2*Pi in x?

x = x-ix*2*Pi;

/* Transformation des Argumentes in [0,Pi/2), Symmetrien ausnutzen */

if(x>=1.5*Pi) // 4. Quadrant

{

vorz = -1;

x = 2*Pi-x;

}

else

if(x>=Pi) // 3. Quadrant

{

vorz = -1;

x = x-Pi;

}

else

if(x>=0.5*Pi) // 2. Quadrant

x = Pi-x;

/*Berechnung von sin(x)*/

x2 =x*x;

for(i=0,sinx=0.0,s=a=x;(sinx!=s)&&(i<n);i=i+1)

s = (sinx=s)+(a=-a*x2/((2*i+2)*(2*i+3)));

sinx = vz*vorz*sinx;

printf("\nErgebnisse\n");

printf("x =%19.15lf\n",xe);

printf("xt=%19.15lf\n\n",x);

printf("sin(%19.15lf)=%19.15lf\n",xe,sin(xe));

printf("Partialsumme s(n)=s(%i)=%19.15lf\n",n,sinx);

getch();

return 0;

}

-----------------------------------------------------------------------------------

// summe5_sin.cpp


// nach Transformation von x auf das Intervall [0,Pi/2]

#include <conio.h>

int main()

{

const double Pi = 3.1415926535897932;

int i,n,vz,ix,vorz;

double xe,x,x2,sinx,s,a;

...

}

104


Potenzreihe sin(x)=x-x∧3/3!+x∧5/5!-x∧7/7!+x∧9/9!-...nach Transformation von x auf das Intervall [-Pi,Pi]

x = 0.5

Ergebnissex = 0.500000000000000xt = 0.500000000000000


————————————————————————————————————-x = 2.2

Ergebnissex = 2.200000000000000xt = 0.941592653589793


————————————————————————————————————-x = 3.5

Ergebnissex = 3.500000000000000xt = 0.358407346410207


————————————————————————————————————-x = 6

Ergebnissex = 6.000000000000000xt = 0.283185307179586


————————————————————————————————————-x = 100

Ergebnissex =100.000000000000000xt = 0.503964914873373


————————————————————————————————————-x = -99

Ergebnissex = -99.000000000000000xt = 1.530964914873373

sin( -99.000000000000000)= 0.999206834186354Partialsumme s(n)=s(100)= 0.999206834186353

105

(5) Integration, Quadraturformeln

// integrate1.c

// Zusammengesetzte Newton-Cotes-Formel: Trapezregel, Simpson-Regel

#include <math.h>

// Integrand

double f(double x)

{

return exp(x);

}

// Zusammengesetzte Trapezregel, Integrand global

double trapez_z(double a,double b,int n)

{

int i;

double h,s1;

h = (b-a)/n;

s1 = 0.5*(f(a)+f(b));

for(i=1;i<=n-1;i++) s1 = s1+f(a+i*h);

return s1*h;

}

// Zusammengesetzte Simpson-Regel, Integrand global

double simpson_z(double a,double b,int n)

{

int i;

double h,s1;

h = (b-a)/(2*n);

s1 = f(a)+f(b);

i = 2;

while(i<=(2*n-2))

{

s1 = s1+2*f(a+i*h);

i = i+2;

}

i = 1;

while(i<=(2*n-1))

{

s1 = s1+4*f(a+i*h);

i = i+2;

}

return s1*h/3;

}

int main()

{

int n;

double a,b,T1,T2,S;

a = 0.0; // linke und rechte Intervallgrenze

b = 1.0;

n = 10; // Anzahl der Integrationsintervalle

printf("Bestimmtes Integral Int(f(x),x=a..b)\n");

printf("Numerische Integration mittels Newton-Cotes-Quadraturformeln\n");

printf("Zusammengesetzte Trapezregel, Simpson-Regel\n\n");

printf("Int(exp(x),x=0..1) = exp(1)-1 = %18.15lf\n",exp(1.0)-1.0);

printf("Anzahl der Teilintervalle n = %i\n",n);

T1 = trapez_z(a,b,n);

printf("T_z(n) = %12.9lf\n",T1);

T2 = trapez_z(a,b,2*n);

printf("T_z(2n) = %12.9lf\n",T2);

S = simpson_z(a,b,n);

printf("S_z(n) = %12.9lf\n",S);

printf("\nKontrolle S_z(n)=(4T_z(2n)-T_z(n))/3 \n");

printf("(4T_z(2n)-T_z(n))/3 = %12.9lf\n",(4*T2-T1)/3);

getch();

}

106

2 C:\d\Neundorf\nwptexte\tech phy\05\cpp gpp gcc\integrate1.exe

Bestimmtes Integral Int(f(x),x=a..b)Numerische Integration mittels Newton-Cotes-QuadraturformelnZusammengesetzte Trapezregel, Simpson-Regel

Int(exp(x),x=0..1) = exp(1)-1 = 1.718281828459045Anzahl der Teilintervalle n = 10T z(n) = 1.719713491T z(2n) = 1.718639789S z (n) = 1.718281888

Kontrolle S z(n)=(4T z(2n)-T z(n))/3(4T z(2n)-T z(n))/3 = 1.718281888

// integrate2.c

// Duale adaptive Trapezregel mit Grob- und Feinrechnung

#include <math.h>

// globale Groessen

int anz;

double hmin;

// Definition des Integranden

double f(double x)

{

return 1.0/(1e-4+x*x);

// andere Integranden

// return sin(x);

// return sqrt(sqrt(sqrt(fabs(x))));

}

// Bestimmtes Integral mit dualer adaptiver Trapezregel

void Integral(double x1,double x2,double f1,double f2,double I,double h,double eps,double *s)

// hmin,anz,f global

{

double I1,I2,I12,xm,fm;

xm = (x1+x2)*0.5;

fm = f(xm);

anz++;

h = 0.5*h; if(h<hmin) hmin = h;

I1 = 0.5*h*(f1+fm);

I2 = 0.5*h*(fm+f2);

I12 = I1+I2;

if(fabs(I-I12)<eps*(fabs(I)+eps)) (*s) = (*s)+I12;

else

{

Integral(x1,xm,f1,fm,I1,h,eps,s); // Integral(x1,xm,f1,fm,I1,h,eps,&*s);

Integral(xm,x2,fm,f2,I2,h,eps,s);

}

}

int main()

{

double a,b,eps,fa,fb,h,s,I;


printf("Duale adaptive Trapezregel mit Grob- und Feinrechnung\n\n");

printf("Untere Grenze a = "); scanf("%lf",&a);

printf("Obere Grenze b = "); scanf("%lf",&b);

printf("Genauigkeit eps = "); scanf("%lf",&eps);

printf("\n");

107

// Vorbereitung des rekursiven Prozeduraufrufes

hmin = b-a;

fa = f(a);

fb = f(b);

anz = 2;

h = b-a;

s = 0.0;

I = 0.5*h*(fa+fb);

Integral(a,b,fa,fb,I,h,eps,&s);

printf("Integralwert s = %15.12lf\n",s);

printf("Anzahl der FW-Auswertungen anz = %i\n",anz);

printf("Minimaler Knotenabstand hmin = %11.9lf\n",hmin);

getch();

}

Beispiel

1∫

−1

dx

10−4 + x2= 200 arctan(100) = 312.159 332 021 646 276 20.


Bestimmtes Integral Int(f(x),x=a..b)Duale adaptive Trapezregel mit Grob- und Feinrechnung

Untere Grenze a = -1Obere Grenze b = 1Genauigkeit eps = 1e-6

Integralwert s = 312.159350642377Anzahl der FW-Auswertungen anz = 16097Minimaler Knotenabstand hmin = 0.000007629————————————————————————————————————-Genauigkeit eps = 1e-9




Integralwert s = 312.159332021827Anzahl der FW-Auswertungen anz = 156574921Minimaler Knotenabstand hmin = 0.000000001

Die letzten Rechnungen mit sehr kleiner Toleranz dauern im Sekundenbereich.Der Fehler wird auf Grund der extrem vielen Schritte dabei wieder zunehmen (Fehlerakku-mulation).

108

// integrate3.c

// Zusammengesetzte Trapezregel mit Grob- und Feinrechnung sowie Intervallhalbierung

#include <math.h>

// Prototypen

double integrate1(double(*funktion)(double),double,double,double);

double integrate2(double(*funktion)(double),double,double,double);

double myf(double);

int main()

{

double a,b,eps,I;

a = 0.0;

b = 1.0;

eps = 1e-6;


printf("Zusammengesetzte Trapezregel mit Grob- und Feinrechnung sowie Intervallhalbierung\n\n");

I = integrate1(myf,a,b,eps);

printf("While-Schleife: I1 = %12.9lf\n",I);

I = integrate2(myf,a,b,eps);

printf("Do-While-Schleife: I2 = %12.9lf\n",I);

getch();

}

/*---------------------------------------------------------------------------*/

double integrate1(double(*funktion)(double x),double a,double b,double eps)

{

double Ialt,Ineu,xi,dx,diff;

int i,N;

/* Erste Naeherung, einfache Trapezregel */

Ialt = (b-a)*0.5*((*funktion)(a)+(*funktion)(b));

/* Zweite Naeherung */

N = 2;

dx = (b-a)/N;

xi = a+dx;

Ineu = dx*0.5*((*funktion)(a)+2*(*funktion)(xi)+(*funktion)(b));

/* Iteration mit Grob- und Feinrechnung */

diff = fabs(Ineu-Ialt);

while(diff>eps) // While-Schleife

{

Ialt = Ineu;

N = N*2;

dx = (b-a)/N;

Ineu = 0.0;

for(i=0;i<N;i++)

{

xi = a+i*dx;

Ineu = Ineu+dx*0.5*((*funktion)(xi)+(*funktion)(xi+dx));

}

diff = fabs(Ineu-Ialt);

}

return Ineu;

}

double integrate2(double(*funktion)(double x),double a,double b,double eps)

{

double Ialt,Ineu,xi,dx,diff;

int i,N;

/* Erste Naeherung, einfache Trapezregel ohne Intervallbreite*/

Ineu = 0.5*((*funktion)(a)+(*funktion)(b));

N = 1;

/* Iteration mit Grob- und Feinrechnung */

do // Do-While-Schleife

{

Ialt = Ineu;

N = N*2;

dx = (b-a)/N;

for(i=1;i<N;i=i+2)

{

109

xi = a+i*dx;

Ineu = Ineu+(*funktion)(xi); // rekursive Berechnung

// damit weniger FW-Berechnungen

}

diff = fabs(Ineu*dx-Ialt*2.0*dx);

}

while(diff>eps);

return Ineu*dx;

}

/*---------------------------------------------------------------------------*/

double myf(double x)

{

return x*x; // Normalparabel

}

Beispiel

b∫

a

x2 dx =1

3(b3 − a3), a = 0, b = 1, ε = 10−6.


Bestimmtes Integral Int(f(x),x=a..b)Zusammengesetzte Trapezregel mit Grob- und Feinrechnung sowie Intervallhalbierung

While-Schleife: I1 = 0.333333492Do-While-Schleife: I2 = 0.333333492

// integrate4.c

// Vergleich verschiedener Integrationsmethoden fuer Integrand = Polynom

// - Zusammengesetzte Trapezregel

// - Monte-Carlo-Methode

// - exaktes Verfahren mit Stammfunktion

#include <stdlib.h>

#include <math.h>

#include <time.h>

#define maxgrad 4 // maximaler Grad des Polynoms

// Prototypen

double direkt(double,double);

double Trapezregel(double(*funktion)(double),double,double,double);

double Monte_Carlo(double(*funktion)(double),double,double,double);

double polynom(double);

/* Globale Variablen */

double ai[maxgrad+1]; // Koeffizienten des Polynoms

int grad; // Grad des Polynoms

/*---------------------------------------------------------------------------*/

int main()

{

double a,b,eps,Id,It,Im;

int i,NN,td,tt,tm;

time_t ta,te;

printf("Berechnung des bestimmten Integrals eines Polynoms Int(pn(x),x=a..b)\n");

printf("Das Integral wird NN mal berechnet, um Laufzeitunterschiede zu ermitteln.\n\n");

printf("Eingaben\n");

printf("Untere Grenze a = "); scanf("%lf",&a);

printf("Obere Grenze b = "); scanf("%lf",&b);

printf("Genauigkeit eps = "); scanf("%lf",&eps);

printf("\nPolynom\n");

do

{

printf("Grad des Polynoms (max. %i. Grades) : ",maxgrad);

110

scanf("%i",&grad);

} while ((grad<0)||(grad>maxgrad));

printf("Polynomkoeffizienten a[i] von pn(x)=a[0]+a[1]x+...+a[n]x^n\n");

for(i=0;i<=grad;i++)

{

printf("a[%i] = ",i);

scanf("%lf",&ai[i]);

}

for(i=grad+1;i<=maxgrad;i++) ai[i] = 0.0;

// mit Zeitmessung, N hinreichend gross waehlen fuer Zeit in sec

NN = 100;

printf("\nBerechnungen laufen\n");

printf("Anzahl der Durchlaeufe NN = %i\n\n",NN);

ta = time(NULL);

for(i=1;i<=NN;i++) Id = direkt(a,b);

te = time(NULL);

td = te-ta;

ta = time(NULL);

for(i=1;i<=NN;i++) It = Trapezregel(polynom,a,b,eps);

te = time(NULL);

tt = te-ta;

ta = time(NULL);

for(i=1;i<=NN;i++) Im = Monte_Carlo(polynom,a,b,eps);

te = time(NULL);

tm = te-ta;

printf("Integrationsmethoden\n");

printf("Direkt Id = %12.9lf Zeit in sec = %ld\n",Id,td);

printf("Tapezregel It = %12.9lf Zeit in sec = %ld\n",It,tt);

printf("Monte-Carlo Im = %12.9lf Zeit in sec = %ld\n",Im,tm);

getch();

}

/*---------------------------------------------------------------------------*/

double direkt(double a,double b)

{

int i;

double I;

I = 0.0;

for(i=0;i<=grad;i++) I = I+ai[i]/(i+1)*(pow(b,i+1)-pow(a,i+1));

return I;

}

/*---------------------------------------------------------------------------*/

double Trapezregel(double (*funktion)(double x),double a,double b,double eps)

{

int i,N;

double Ialt,Ineu,xi,dx;

// Erste Naeherung, einfache Trapezregel

Ineu = 0.5*((*funktion)(a)+(*funktion)(b));

N = 1;

// Iteration

do

{

Ialt = Ineu;

N = N*2;

dx = (b-a)/N;

for(i=1;i<N;i=i+2)

{

xi = a+i*dx;

Ineu = Ineu+(*funktion)(xi);

}

} while(fabs(Ineu*dx-Ialt*2.0*dx)>eps);

return Ineu*dx;

}

/*---------------------------------------------------------------------------*/

double Monte_Carlo(double(*funktion)(double x),double a,double b,double eps)

{

const int N = 100;

111

int i,anzp;

double fmittel,fmittelalt,xr,r;

//srand48(1357); // Zufallszahlengenerator initialisieren

srand(1357);

r = (double)RAND_MAX;

fmittel = 0.0;

for(i=1;i<=N;i++)

{

//xr = drand48();

xr = a+(b-a)*rand()/r;

fmittel = fmittel+(*funktion)(xr);

}

anzp = N;

do

{

fmittelalt = fmittel;

for(i=1;i<=N;i++)

{

//xr = drand48();

xr = a+(b-a)*rand()/r;

fmittel = fmittel+(*funktion)(xr);

}

anzp = anzp+N;

// printf("I: %lf AP: %i\n", fmittel/anzp*(b-a), anzp);

} while(fabs(fmittel/anzp-fmittelalt/(anzp-N))*(b-a)>eps);

return fmittel/anzp*(b-a);

}

/*---------------------------------------------------------------------------*/

double polynom(double x)

{

int i;

double p;

p = 0.0;

for(i=0;i<=grad;i++) p = p+ai[i]*pow(x,i);

return p;

}

----------------------------------------------------------------------------------------

// integrate4.cpp

// Vergleich verschiedener Integrationsmethoden fuer Integrand = Polynom

// - Zusammengesetzte Trapezregel

// - Monte-Carlo-Methode

// - exaktes Verfahren mit Stammfunktion

#include <stdlib.h>

#include <conio.h>

#include <math.h>

#include <time.h>

#define maxgrad 4 // maximaler Grad des Polynoms

// Prototypen

double direkt(double,double);

double Trapezregel(double(*funktion)(double),double,double,double);

double Monte_Carlo(double(*funktion)(double),double,double,double);

double polynom(double);

/* Globale Variablen */

double ai[maxgrad+1]; // Koeffizienten des Polynoms

int grad; // Grad des Polynoms

/*---------------------------------------------------------------------------*/

int main()

{

double a,b,eps,Id,It,Im;

int i,NN,td,tt,tm;

time_t ta,te;

...

}

...

112

Beispiel

b∫

a

pn(x) dx =

b∫

a

(a0 + a1x + ... + anxn) dx

=[

a0x +1

2a1x

2 + ... +1

n + 1anx

n+1

]b

a

=n∑

i=0

1

i + 1ai (b

i+1 − ai+1).


Berechnung des bestimmten Integrals eines Polynoms Int(pn(x),x=a..b)Das Integral wird NN mal berechnet, um Laufzeitunterschiede zu ermitteln.

EingabenUntere Grenze a = 0Obere Grenze b = 1Genauigkeit eps = 1e-6

PolynomGrad des Polynoms (max. 4. Grades) : 3Polynomkoeffizienten a[i] von pn(x)=a[0]+a[1]x+...+a[n]x∧na[0] = 1a[1] = 1a[2] = 1a[3] = 1

Berechnungen laufenAnzahl der Durchlaeufe NN = 100

IntegrationsmethodenDirekt (exakt) Id = 2.083333333 Zeit in sec = 0 (NN=200: 0, NN=400: 0, NN=600: 0)Tapezregel It = 2.083333433 Zeit in sec = 0 ( 0, 1, 1)Monte-Carlo Im = 2.072050397 Zeit in sec = 1 ( 3, 6, 8)———————————————————————————————————————–

...Genauigkeit eps = 1e-8...Anzahl der Durchlaeufe NN = 100

IntegrationsmethodenDirekt (exakt) Id = 2.083333333 Zeit in sec = 0 (NN=200: 0, NN=400: 0, NN=600: 0)Tapezregel It = 2.083333335 Zeit in sec = 1 ( 3, 6, 9)Monte-Carlo Im = 2.080377829 Zeit in sec = 24 ( 46, 95, 140)

Das Verfahren von Monte-Carlo ist im Vergleich mit Abstand am langsamsten undsehr ungenau.

113

dev cpp - Startseite TU Ilmenau · Hinweise zur Nutzung der Programmieroberﬂ¨ache Dev-C++, der...

Documents

Transcript of dev cpp - Startseite TU Ilmenau · Hinweise zur Nutzung der Programmieroberﬂ¨ache Dev-C++, der...