hallo

abgabe6nochmal.cpp

@@ -0,0 +1,141 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <ctype.h>
+#include <math.h>
+#include <complex>
+>
+#include<vector>
+using namespace std;
+#ifndef M_PI
+#define M_PI           3.14159265358979323846
+#endif
+typedef struct {
+	long N;
+	double b, a, beta, alpha;
+	double tmin, tmax,h;
+	int m;
+    int sign;
+} fft_param;
+
+complex<double> expo(fft_param *fp,int k) {
+    int sign = fp->sign;
+	return exp(sign*2 * M_PI * 1i / ((double) k));
+}
+
+vector<complex<double>> fourierTrafo(fft_param *fp, vector<complex<double>> x) {
+	long N = x.size();
+    int sign=fp->sign;
+	vector<complex<double>> gerPkt(N / 2);
+	vector<complex<double>> ungerPkt(N / 2);
+    
+	if (N == 1) {
+		
+		return x;
+	} else {
+		for (int k = 0; k < N / 2; k++) {
+			gerPkt[k] = x[2 * k];
+			ungerPkt[k] = x[2 * k + 1];
+		}
+		
+		vector<complex<double>> ftger(N / 2) = fourierTrafo(fp, gerPkt);   //???????????????????????
+		vector<complex<double>> ftunger(N / 2) = fourierTrafo(fp, ungerPkt); //???????????????????
+		vector<complex<double>> out(N);
+        
+		for (int k = 0; k < N / 2; k++) {
+			out[k] = (ftger[k] + ftunger[k] * expo(fp,k))/(double)N;
+			out[k + N / 2] = (ftger[k] - ftunger[k] *expo(fp,k))/(double)N;
+		}
+		
+		return out;
+    }
+}
+vector<complex<double>> test(fft_param *fp){
+    double h=fp->h;
+    double tmin=fp->tmin;
+    double tmax=fp->tmax;
+    double t=tmin;
+    int N=fp->N;
+    
+    vector<complex<double>> x(N);
+    
+    for(int k=0;k<N;k++){
+        x(k)=sin(tmin+k*h);
+    }
+    vector<complex<double>> out = fourierTrafo(fp, x);
+    fft_param n= *fp;
+    n.sign=-1;
+    vector<complex<double>> doppelout = fourierTrafo(&n, out);
+    FILE *fp2 = fopen("testFFT.csv" ,'w');
+    for(int k=0;k<N;k++){
+        fprintf(fp2,"%f\t%f\t%f\n",x[k], doppelout[k],(tmin+k*h));
+    }
+    
+    return x;
+}   
+vector<complex<double>> makeS(fft_param *fp){
+    double h=fp->h;
+    double tmax=fp->tmax;
+    double tmin=0;
+    double t;
+    double tl=fp->tl;
+    double alpha=f->alpha;
+    int N=fp->N;
+    vector<complex<double>> s(N);
+    for(int k=0;k<N;k++){
+        t=tmin+k*h;
+        s[k]=sin(2*M_PI*alpha*(t-tl));
+    }
+    return s;
+}
+vector<complex<double>> makeE(fft_param *fp){
+    double h=fp->h;
+    double tmax=fp->tmax;
+    double tmin=0;
+    double t;
+    double beta=f->beta;
+    double tl=fp->tl;
+    double a=fp->a;
+    double b=fp->b;
+    int N=fp->N;
+    vector<complex<double>> e(N);
+    vector<complex<double>> s=makeS(fp);
+    for(int k=0;k<N;k++){
+        t=tmin+k*h;
+       e[k]=b*sin(2*M_PI*beta*t)+a*rand(1)+s[k];
+    }
+    return e;
+}
+vector<complex<double>> aufgabe3 (fft_param *fp){
+    double h=fp->h;
+    double tmin=fp->tmin;
+    double tmax=fp->tmax;
+    double t=tmin;
+    int N=fp->N;
+    vector<complex<double>> e(N)=makeE;
+    vector<complex<double>> s(N)=makeS(fp);
+    vector<comlex<double>> Fe=fourierTrafo(fp,e);
+    fft_param fps=*fp;
+    fps.sign=-1;
+    vector<comlex<double>> Fs=fourierTrafo(&fps,s);
+    vector<comlex<double>> FFes(N);
+    
+    for(int k=0;k<N;k++){
+        FFes[k]=Fe[k]*conj(Fs[k]);
+        
+    }
+    vector<comlex<double>> FFFes(N)=fourierTrafo(fp, FFes);
+    return FFFes;
+    
+}
+double main(){
+    fft_param fp;
+    fp.N= pow(2,13);
+    fp.tmin=0;
+    fp.tmax = 100;
+    fp.h=(fp.tmax-fp.tmin)/fp.N;
+    fp.sign=1;
+    test(&fp);
+    
+    
+}

fuenfteabgabe.c

@@ -1,60 +1,90 @@
 #include <stdio.h>
 #include <stdlib.h>
+#include <iostream>
 #include <math.h>
+#include <string.h>
+using namespace std;
 typedef struct {
-	double N;
-	double d;
+	double N; //Parameter anzahl der Solitonen
+	double d; //Schrittweiten von Zeit und Ort
 	double h;
-	double xmin, xmax, tmin, tmax;
-	double j;
-	double n;
-	int jmax;
+	double xmin, xmax, tmin, tmax; //Intervalle
+	double x; //momentanes x und t
+	double t;
+	int jmax; //Maximale Anzahl der Schritte
 	int nmax;
 
 } soli_param;
-double *u0;
-double *uj1;
+char dateiname[256]="plotu1";
+double *u0; //Pointer auf denen u(n=0,j) gespeichert wird
+double *uj1; //Pointer auf denne u(n=1,j) gespeichert wird
+//u0 mit Glg (3) berechnet
 void ustart(soli_param *sp) {
-	double N = sp->N;
+	double N = sp->N; //Variablen aus Struct holen
 	double jmax = sp->jmax;
 	double h = sp->h;
 	double x = sp->xmin;
-    //FILE *fp2=fopen("cosh.csv","w");
+    
 	for (int j = 0; j <= jmax; j++) {
 
 		u0[j] = -N * (N + 1) / (pow(cosh(x), 2));
-		
-		//fprintf(fp2,"%f\t%f\n", x,u0[j]);
-        x = x + h;
+		x=x+h;
 
 	}
 
 }
+//analytische Lösungen berechnen
+double u_analytic(soli_param *sp){
+    double x=sp->x; //Variablen aus Struct holen
+    double t=sp->t;
+    double N= sp->N;
+    //für verschiedene N
+    if (N==1){
+        return (-2/pow(cosh(x-4*t),2));
+    }
+    else if (N==2){
+        return (-12*(3+4*cosh(2*x-8*t)+cosh(4*x-64*t))/(pow((3*cosh(x-28*t)+cosh(3*x-36*t)),2)));
+    }
+    else{
+        printf("Keine analytische Lösung für dieses N");
+        return 0;
+    }
+}
+//uj1 nach Glg (8) berechnen
 void ujeins(soli_param *sp) {
-	double jmax = sp->jmax;
+	double jmax = sp->jmax; //Variablen aus Struct holen
 	double d = sp->d;
 	double h = sp->h;
     double x=sp->xmin;
+    //für alle j
 	for (int j = 0; j < jmax; j++) {
+        //Randbed: u0(j-2) fällt weg wenn j=1 ist 
 		if (j == 1) {
 			uj1[j] = (d * 6 * u0[j] * (u0[j + 1] - u0[j - 1]) / (2 * h))
 					- (d * (u0[j + 2] + 2 * u0[j - 1]-2*u0[j+1]) / (2 * pow(h, 3)))
 					+ u0[j];
 
-		} else if (j == 0) {
+		}
+		//Randbed: u0(j-2) und u0(j-1) fällt weg wenn j=0 ist 
+		else if (j == 0) {
 			uj1[j] = (d * 6 * u0[j] * (u0[j + 1]) / (2 * h))
 					- (d * (u0[j + 2]-2*u0[j+1]) / (2 * pow(h, 3))) + u0[j];
 
-		} else if (j == jmax - 2) {
+		}
+		//Randbed: u0(j+2) fällt weg wenn j=jmax-2 ist 
+		else if (j == jmax - 2) {
 			uj1[j] = (d * 6 * u0[j] * (u0[j + 1] - u0[j - 1]) / (2 * h))
 					- (d * (-2*u0[j+1]+2 * u0[j - 1] - u0[j - 2]) / (2 * pow(h, 3)))
 					+ u0[j];
 
-		} else if (j == jmax-1) {
+		}
+		//Randbed: u0(j+2) und u0(j+1) fällt weg wenn j=jmax -1 ist 
+		else if (j == jmax-1) {
 			uj1[j] = (d * 6 * u0[j] * (-u0[j - 1]) / (2 * h))
 					- (d * (2 * u0[j - 1] - u0[j - 2]) / (2 * pow(h, 3)))
 					+ u0[j];
-		} else {
+		}
+		else {
 			uj1[j] = ((d * 6 * u0[j] * (u0[j + 1] - u0[j - 1])) / (2 * h))
 					- (d * (u0[j + 2] -2*u0[j+1]+ 2 * u0[j - 1] - u0[j - 2])
 							/ (2 * pow(h, 3))) + u0[j];
@@ -69,109 +99,105 @@ void nextu(soli_param *sp) {
 	int nmax = sp->nmax;
 	double d = sp->d;
 	double h = sp->h;
-	double * unm = u0;
-	double * un = uj1;
-	double unp[jmax];
+	double * unm = u0; //unm=u(n-1) wird auf u0 gesetzt 
+	double * un = uj1; //analog dazu auf uj1
+	double unp[jmax]; //unp=u(n+1) ist also im ersten schritt in der Schleife u(n=2)
     double t=sp->tmin;
     double x=sp->xmin;
     int n;
-    
+    //Berechnen für alle n bis tmax erreicht ist
 	for ( n = 2; n < nmax; n++) {
         
 		for (int j = 0; j < jmax; j++) {
+            //Randbedingungen werden genauso benutzt wie in ujeins
 			if (j == 1) {
-				/*unp[j] = 2*d*(
-								(un[j + 1] + un[j] + un[j - 1])
-										* (un[j + 1] - un[j - 1]) / h
-										- (un[j + 2] - 2 * un[j + 1]
-												+ 2 * un[j - 1])
-												/ (2 * pow(h, 3))) + unm[j];*/
+				
                 unp[j]=(2*(d/h) *(un[j+1]+un[j]+un[j-1])*(un[j+1]-un[j-1]))-(d/((h*h*h))*(un[j+2]-2*un[j+1]+2*un[j-1]) ) +unm[j];  
 
 			} else if (j == 0) {
-				/*unp[j] = 2
-						* d*(
-								(un[j + 1] + un[j]) * (un[j + 1]) / h
-										- (un[j + 2] - 2 * un[j + 1])
-												/ (2 * pow(h, 3))) + unm[j];*/
                 unp[j]=(2*(d/h) *(un[j+1]+un[j])*(un[j+1]))-(d/(h*h*h)*(un[j+2]-2*un[j+1]) ) +unm[j];  
 
 			} else if (j == jmax - 2) {
-				/*unp[j] = 2
-						* d*(
-								(un[j + 1] + un[j] + un[j - 1])
-										* (un[j + 1] - un[j - 1]) / h
-										- (-2 * un[j + 1] + 2 * un[j - 1]
-												- un[j - 2]) / (2 * pow(h, 3)))
-						+ unm[j];*/
                 unp[j]=(2*(d/h) *(un[j+1]+un[j]+un[j-1])*(un[j+1]-un[j-1]))-(d/((h*h*h))*(-2*un[j+1]+2*un[j-1]-un[j-2]) ) +unm[j];  
 
 			} else if (j == jmax - 1) {
-				/*unp[j] = 2
-						* d*(
-								(un[j] + un[j - 1]) * (-un[j - 1]) / h
-										- (2 * un[j - 1] - un[j - 2])
-												/ (2 * pow(h, 3))) + unm[j];*/
                 unp[j]=(2*(d/h) *(un[j]+un[j-1])*(-un[j-1]))-(d/((h*h*h))*(2*un[j-1]-un[j-2]) ) +unm[j];  
 
 			} else {
-				/*unp[j] = 2
-						* d*(
-								(un[j + 1] + un[j] + un[j - 1])
-										* (un[j + 1] - un[j - 1]) / h
-										- (un[j + 2] - 2 * un[j + 1]
-												+ 2 * un[j - 1] - un[j - 2])
-												/ (2 * pow(h, 3))) + unm[j];*/
-                                                
-                 unp[j]=(2*(d/h) *(un[j+1]+un[j]+un[j-1])*(un[j+1]-un[j-1]))-(d/((h*h*h))*(un[j+2]-2*un[j+1]+2*un[j-1]-un[j-2]) ) +unm[j];                     
+				unp[j]=(2*(d/h) *(un[j+1]+un[j]+un[j-1])*(un[j+1]-un[j-1]))-(d/((h*h*h))*(un[j+2]-2*un[j+1]+2*un[j-1]-un[j-2]) ) +unm[j];                     
 			}
-			//unm[j]=un[j];
-            //un[j]=unp[j];
-			
-
 		}
+		//Variablen fuer den nächsten Schritt tauschen (komponentenweise, da Pointer)
 		for (int j = 0; j < jmax; j++) {
     		  unm[j]=un[j];
     		  un[j]=unp[j];
             }
 		
-		
-        
-        
 
 	}
-	FILE *fp= fopen("numplot8.csv", "w");
-	printf("%f\t%f\n", d*n, d*nmax);
+    //Speichern der Werte
+    soli_param diffp= *sp; //Neue Variablen in struct um neues x und t u_analytic mitzgeben
+    diffp.t=sp->tmax;
+    diffp.x=sp->xmin;
+    double diff;
+	FILE *fp= fopen("numplot0_1_2.csv", "w"); //Datei öffnen
 	for(int j=0;j<jmax;j++){
-		
-		fprintf(fp, "%f\t%f\n",x,un[j]);
-		//printf( "%f\t%f\n",x,un[j]);
-        x=x+h;
-
+		diff= fabs(un[j]-u_analytic(&diffp)); //Differenz berechnen
+		fprintf(fp, "%f\t%f\t%f\t%f\n",diffp.x, un[j],u_analytic(&diffp),diff); //In datei ausgeben
+        printf("%f\t%f\t%f\t%f\n",diffp.x, un[j],u_analytic(&diffp),diff);
+        diffp.x=diffp.x+h; //neues x berechnen
 	}
 	fclose(fp);
 }
-
+//analytische Lösung fuer H10.1 speichern
+void plotu_analytic(soli_param * sp){
+	FILE *fp =fopen(dateiname, "w"); //Dateiöffnen
+	double f; //Variable zum speichern des Wertes
+	double x=sp->xmin;
+    double xmax=sp->xmax;
+    double xmin=x;
+    soli_param ua= *sp; //neuer struct zum ändern der Parameter
+    ua.x=x;
+    ua.t=sp->t;
+	while(ua.x<=xmax){
+		f=u_analytic(&ua); //rechnen
+		fprintf(fp, "%f\t%f\n", ua.x,f); //in Datei schreiben
+		ua.x=ua.x+0.1;
+
+	} 
+}	
 
 int main(void) {
-    soli_param sp;
+    soli_param sp; //definieren der Variablen im struct
     sp.N = 1;
     sp.h=0.1;
-    sp.d=0.0001;
+    sp.d=pow(sp.h,3)/10;
     sp.tmax=1;
     sp.tmin=0;
     sp.xmin=-30;
     sp.xmax=30;
     sp.jmax=(sp.xmax-sp.xmin)/sp.h;
     sp.nmax=(sp.tmax-sp.tmin)/sp.d;
-    //sp.nmax=500;
-    printf("%d\n", sp.nmax);
-    u0 = (double*) malloc(sp.jmax * sizeof(double));
-    ustart(&sp);
+    //H10.1:
+    /*sp.t=-1;
+    sp.xmin=-20;
+    sp.xmax=20;
+	int i=0;
+	while(sp.t<=1){
+		sprintf(&dateiname[0],"%s%d","plotu1",i);
+		plotu_analytic(&sp);
+		sp.t+=0.1;
+		i++;
+	}*/
+	//H10.2
+	sp.N=2;
+    sp.xmin=-30;
+    sp.xmax=30;
+    u0 = (double*) malloc(sp.jmax * sizeof(double)); //speicher allocieren
+    ustart(&sp); //ustart aufrufen um u0 zu berechnen
     uj1 = (double*) malloc(sp.jmax * sizeof(double));
-    ujeins(&sp);
-    nextu(&sp);
-    free(u0);
+    ujeins(&sp); //ujeins aufrufen um uj1 zu berechnen
+    nextu(&sp); //berechnet u für t=1
+    free(u0); //speicher befreien
     free(uj1);
-
 }