ModSim ass3: worked on report.

modsim/ass3/main.c

@@ -5,7 +5,7 @@
 #include "models.h"
 #include "logger.h"
 
-#define DEFAULT_DT .001
+#define DEFAULT_DT .01
 #define COMPARE_DT 1.0
 
 #define INTEGRATE(method, t0, t1, dt, N, FN) { \
@@ -28,26 +28,38 @@ int main(int argc, char **argv) {
            args_forced_osc[3]  = {1.0, 0.1, 2.0},              // {k, D, w}
            args_lottka_volt[4] = {.5, 1.0, .1, .1},            // {a, b, c, d}
            args_modified_pp[5] = {.5, 1.0, .1, .1, 100.0},     // {a, b, c, d, y_m}
-           args_pid_control[5] = {-2.0, -1.0, .1, 10.0, .001}, // {P, D, a, M, d}
+           args_pid_control[5] = {-2.0, -1.0, .1, 10.0, 0.1}, // {P, D, a, M, d}
            *args_gilpin = NULL;
 
     switch( atoi(argv[1]) ) {
         case 0:
+            y0[0] = 10.0;
+            y0[1] = 20.0;
             logfile = fopen("report/osc_euler", "w");
             puts("Integrating Euler method...");
-            INTEGRATE(Euler, .0, 40.0, COMPARE_DT, 2, osc);
+            INTEGRATE(Euler, .0, 100.0, COMPARE_DT, 2, osc);
             logger_close();
             logfile = fopen("report/osc_rk2", "w");
             puts("Integrating Runge-Kutta 2 method...");
-            INTEGRATE(RungeKutta2, .0, 40.0, COMPARE_DT, 2, osc);
+            INTEGRATE(RungeKutta2, .0, 100.0, COMPARE_DT, 2, osc);
             logger_close();
             logfile = fopen("report/osc_rk4", "w");
             puts("Integrating Runge-Kutta 4 method...");
-            INTEGRATE(RungeKutta4, .0, 40.0, COMPARE_DT, 2, osc);
+            INTEGRATE(RungeKutta4, .0, 100.0, COMPARE_DT, 2, osc);
         break;
-        case 1: INTEGRATE(RungeKutta4, .0, 20.0, DEFAULT_DT, 2, forced_osc); break;
-        case 2: INTEGRATE(RungeKutta4, .0, 20.0, DEFAULT_DT, 2, lottka_volt); break;
-        case 3: INTEGRATE(RungeKutta4, .0, 20.0, DEFAULT_DT, 2, modified_pp); break;
+        case 1: INTEGRATE(RungeKutta4, .0, 30.0, DEFAULT_DT, 2, forced_osc); break;
+        case 2:
+            y0[0] = 0.5;
+            y0[1] = 1.0;
+            //y0[0] = 50.0;
+            //y0[1] = 100.0;
+            INTEGRATE(RungeKutta4, .0, 400.0, DEFAULT_DT, 2, lottka_volt); break;
+        case 3:
+            y0[0] = 0.5;
+            y0[1] = 1.0;
+            //y0[0] = 50.0;
+            //y0[1] = 100.0;
+            INTEGRATE(RungeKutta4, .0, 100.0, DEFAULT_DT, 2, modified_pp); break;
         case 4:
             y0[0] = .0;  // = v(0)
             y0[1] = 3.0; // = s(0)
@@ -60,7 +72,7 @@ int main(int argc, char **argv) {
                 s_buf[i] = y0[1];
             }
 
-            INTEGRATE(RungeKutta4, .0, 100.0, DEFAULT_DT, 2, pid_control);
+            INTEGRATE(RungeKutta4, .0, 150.0, DEFAULT_DT, 2, pid_control);
             free(v_buf);
             free(s_buf);
         break;

modsim/ass3/plot.py

@@ -19,7 +19,7 @@ for i in range(len(lines)):
 figure(1)
 clf()
 for i in range(len(y)):
-    plot(x, y[i], 'x')
+    plot(x, y[i], '-')
 if len(argv) == 2:
     savefig(argv[1])
 show()

modsim/ass3/report/report.tex

@@ -25,7 +25,8 @@
 
 De opdracht bestaat uit twee delen: het eerste deel betreft het implementeren
 van de Euler- en Runge-Kutta-methoden, in het tweede deel worden deze methoden
-vergeleken en gebruikt om enkele ODE's te benaderen.
+vergeleken en gebruikt om enkele ODE's te benaderen. De opdrahchten zelf zullen
+we niet heel uitgebreid bespreken, alleen de antwoorden die we erop hebben.
 
 \section{Implementatie}
 
@@ -78,7 +79,7 @@ De vergelijking van de drie integratiemethoden kan worden uitgevoerd met:
 $ sh ./compare_osc.sh
 \end{verbatim}
 
-In de volgende grafieken staan hiervan de resultaten (blauw
+
 
 \begin{figure}[H]
     \includegraphics[scale=.5]{osc_euler.pdf}
@@ -99,15 +100,81 @@ In de volgende grafieken staan hiervan de resultaten (blauw
 
 \subsection{Lottka-Voltera}
 
+We hebben de Lotka-Volterra simulatie met de opgegeven constanten uitgevoerd
+voor verschillende startwaarden. De ene keer met startwaarden op het ``stabiele
+punt'' en de andere keer met startwaarden daar ver vandaan. Het stabiele punt
+ligt op $y = \frac{a}{b} = 0.5$ en $x = \frac{c}{d} = 1.0$.
+
+\begin{figure}[H]
+    \includegraphics[scale=.5]{lot_vol_close.pdf}
+    \caption{$x(0) = 0.5$, $y(0) = 1.0$}
+\end{figure}
+
+\begin{figure}[H]
+    \includegraphics[scale=.5]{lot_vol_far.pdf}
+    \caption{$x(0) = 100$, $y(0) = 200$}
+\end{figure}
 
+We zien dat de grafiek bij de hoge waardes een kleiner periode heeft, dus vaker
+fluctueert. Bij het stabiele punt zijn er grotere periodes, waardoor de prooi
+sneller een groot aantal krijgt.
 
 \subsection{Modified predator-prey model}
 
+We gebruiken exact dezelfde startwaarden als bij Lotka-Volterra en nu krijgen
+het volgende resultaat:
+
+\begin{figure}[H]
+    \includegraphics[scale=.5]{mod_pp_close.pdf}
+    \caption{$x(0) = 0.5$, $y(0) = 1.0$}
+\end{figure}
+
+\begin{figure}[H]
+    \includegraphics[scale=.5]{mod_pp_far.pdf}
+    \caption{$x(0) = 100$, $y(0) = 200$}
+\end{figure}
 
+We zien dat in beide gevallen de grafiek niet boven de 100 uitkomt (namelijk
+omdat $y_m = 100$, dit is de maximumwaarde). De tweede grafiek vertoont na
+verloop van tijd min of meer hetzelfde gedrag als de eerste, alleen wordt in het
+eerste gedeelte de stabiele positie ``opgezocht'' waarna naar dezelfde
+evenwichtsstand wordt bewogen.
+
+Hieruit concluderen we dat de startwaarden waarschijnlijk weinig uitmaken, omdat
+na verloop van de tijd de grafiek toch wel een evenwichtsstand bereikt.
 
 \subsection{PID controller}
 
+De constante P is de mate waarin de afstand tot de doelpositie ($s$) meetelt
+voor de snelheid van de beweging. Als $s(0)$ positief gekozen wordt moet P
+negatief zijn (bijv. $3 > 0$, als $s(0) = 3$, er moet naar $s = 0$ worden
+bewogen dus een negatieve beweging), en andersom.
+
+D is de mate waarin de bewegingssnelheid wordt beperkt, men wil waarschijnlijk niet dat de arm
+een bepaalde snelheid overschrijdt. Omdat D een beperkende factor is, is deze
+ook negatief.
+
+$d$ is de delay voor het bewegen voor de arm, we hebben voor verschillende waarden
+hiervan de grafiek geplot. We hebben hierbij gekozen voor $P = -3$, $D = -1$, $a
+= 0.1$ en $m = 10$ als constanten.
+
+\begin{figure}[H]
+    \includegraphics[scale=.5]{pid_big.pdf}
+    \caption{$d = 1.0$}
+\end{figure}
+
+\begin{figure}[H]
+    \includegraphics[scale=.5]{pid_medium.pdf}
+    \caption{$d = 0.5$}
+\end{figure}
+
+\begin{figure}[H]
+    \includegraphics[scale=.5]{pid_small.pdf}
+    \caption{$d = 0.1$}
+\end{figure}
 
+We zien dat hoe kleiner $d$, hoe vlugger de grafiek een evenwichtsstand bereikt.
+Dit is logisch, want de arm reageert sneller als de delay kleiner is.
 
 \subsection{Gilpin's Model}