140 likes | 221 Views
pendulum.c. Thomas Baumgartner, Mat.nr. 0425162 Alexander Gross, Mat.nr. 0425230. q (t) durch Integration mit CVODE berechnen Phasenraumpunkte (A, w ) bestimmen, für die das Pendel stehen bleibt. Aufgabenstellung. m. q. L. Reduktion auf 1. Ordnung. N_Vector y. N_Vector ydot.
E N D
pendulum.c Thomas Baumgartner, Mat.nr. 0425162 Alexander Gross, Mat.nr. 0425230
q(t) durch Integration mit CVODE berechnen Phasenraumpunkte (A, w) bestimmen, für die das Pendel stehen bleibt Aufgabenstellung m q L
Reduktion auf 1. Ordnung N_Vector y N_Vector ydot
Rechte Seite der Gleichung static void f(integer N, real t, N_Vector y, N_Vector ydot, void *f_data) { double* v_y; double* v_ydot; v_y = N_VDATA(y); v_ydot = N_VDATA(ydot); v_ydot[0] = v_y[1]; v_ydot[1] = (1/L)*(g-A*pow(w,2)*sin(w*t))*sin(v_y[0]); }
(A, w) bestimmen Parameter: A1, A2, w1, w2, Ares, wres, q0 Ausgabedatei: phasespace.dat Betriebsarten q(t) berechnen • Parameter: A, w, q0 • Ausgabedatei: theta.dat
Ablauf: q(t) berechnen • Programm mit gewünschten Parametern aufrufen • Integration von t = 0 bis t = 10s mittels linearem Solver CVSPGMR und Adams-Methode • Ausgabe in Datei theta.dat für jeden berechneten Zeitpunkt
Code: q(t) berechnen cvode_mem = CVodeMalloc(NEQ, f, T0, y, ADAMS, FUNCTIONAL, SS, &reltol, &abstol, fdata, NULL, FALSE, NULL, NULL, NULL); if(cvode_mem == NULL) { printf("CVodeMalloc failed.\n"); return(1); } while(t <= Tend){ flag = CVode(cvode_mem, TSTEP, y, &t, ONE_STEP); if (flag != SUCCESS) { printf("CVode failed, flag = %d.\n", flag); return(1); } fprintf(fout, "%e %e %e\n", t, v_y[0], v_y[1]); ++timeSteps; }
Beispiel: q(t) berechnen ./pendulum 0.2 10 0.02 Calculating Theta(t) with the following parameters: A = 0.150000 w = 70.000000 Th0 = 0.020000 Integration finished. 9994 time steps calculated.
Ablauf: (A, w) bestimmen • Programm mit gewünschten Parametern aufrufen • Schleifen für die Wertebereiche von A und w werden durchlaufen • Integration von t = 0 bis t = 10s für jedes (A, w)-Wertepaar (CVSPGMR mit Adams-Methode)
Ablauf: (A, w) bestimmen • Stabilitätskriterium 1: • Stabilitätskriterium 2: • Ausgabe in Datei phasespace.dat für jedes Wertepaar (A, w) mit stabiler Lösung
Code: (A, w) bestimmen while((t <= Tend)&&(flagStable == 1)){ . . . if (fabs(v_y[0] - Th0) > ThTOL){ flagStable = 0; } } if (flagStable == 1){ . . . if (fabs(ThAvg - Th0) <= ThAvgTol){ printf("Stable solution found: A = %f, w = %f\n", A, w); fprintf(fout, "%e %e\n", A, w); ++stablePoints; } }
Beispiel: (A, w) bestimmen ./pendulum 0 1 0 100 0.01 0.5 0.02 Searching phase space for stable solutions with the following parameters: A = 0.000000 ... 1.000000 w = 0.000000 ... 100.000000 Ares = 0.010000 wres = 0.500000 Th0 = 0.020000 Stable solution found: A = 0.060000, w = 87.000000 Stable solution found: A = 0.060000, w = 94.000000 . . . Phase space analysis finished. 1388 stable solutions found.