Commit 7c5530c6 authored by Sander Mathijs van Veen's avatar Sander Mathijs van Veen

fixed merge conflict in main.c

parents 3f68304d 5d96e4ca
...@@ -37,11 +37,13 @@ int Euler(double t0, double t1, double dt, double *y0, double *y1, int N, ...@@ -37,11 +37,13 @@ int Euler(double t0, double t1, double dt, double *y0, double *y1, int N,
puts("Error calculating slope."); puts("Error calculating slope.");
return -1; return -1;
} }
printf("%.12e: %.12e,%.12e\n", t, y1[0], y1[1]);
for( i = 0; i < N; i++ ) for( i = 0; i < N; i++ )
y1[i] = y1[i] + dy[i] * dt; y1[i] = y1[i] + dy[i] * dt;
// Value can be printed/logged here // Value can be printed/logged here
printf("%.12e: %.12e,%.12e\n", t, y1[0], y1[1]);
logger(t, N, y1); logger(t, N, y1);
} }
...@@ -130,6 +132,7 @@ int RungeKutta4(double t0, double t1, double dt, double *y0, double *y1, int N, ...@@ -130,6 +132,7 @@ int RungeKutta4(double t0, double t1, double dt, double *y0, double *y1, int N,
} }
// Value can be printed/logged here // Value can be printed/logged here
printf("%.12e: %.12e, %.12e\n", t, y1[0], y1[1]);
logger(t, N, y1); logger(t, N, y1);
} }
......
...@@ -24,4 +24,5 @@ void logger(double t, int N, double *y) { ...@@ -24,4 +24,5 @@ void logger(double t, int N, double *y) {
void logger_close() { void logger_close() {
if( logfile && logfile != stdout && fclose(logfile) ) if( logfile && logfile != stdout && fclose(logfile) )
printf("Failed to close logfile. Error no: %d", errno); printf("Failed to close logfile. Error no: %d", errno);
logfile = NULL;
} }
...@@ -5,7 +5,7 @@ ...@@ -5,7 +5,7 @@
#include "models.h" #include "models.h"
#include "logger.h" #include "logger.h"
#define DEFAULT_DT .001 #define DEFAULT_DT .01
#define COMPARE_DT 1.0 #define COMPARE_DT 1.0
#define INTEGRATE(method, t0, t1, dt, N, FN) { \ #define INTEGRATE(method, t0, t1, dt, N, FN) { \
...@@ -28,26 +28,38 @@ int main(int argc, char **argv) { ...@@ -28,26 +28,38 @@ int main(int argc, char **argv) {
args_forced_osc[3] = {1.0, 0.1, 2.0}, // {k, D, w} 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_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_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; *args_gilpin = NULL;
switch( atoi(argv[1]) ) { switch( atoi(argv[1]) ) {
case 0: case 0:
y0[0] = 10.0;
y0[1] = 20.0;
logger_open("report/osc_euler"); logger_open("report/osc_euler");
puts("Integrating Euler method..."); puts("Integrating Euler method...");
INTEGRATE(Euler, .0, 40.0, COMPARE_DT, 2, osc); INTEGRATE(Euler, .0, 100.0, COMPARE_DT, 2, osc);
logger_close(); logger_close();
logger_open("report/osc_rk2"); logger_open("report/osc_rk2");
puts("Integrating Runge-Kutta 2 method..."); 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(); logger_close();
logger_open("report/osc_rk4"); logger_open("report/osc_rk4");
puts("Integrating Runge-Kutta 4 method..."); 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; break;
case 1: INTEGRATE(RungeKutta4, .0, 20.0, DEFAULT_DT, 2, forced_osc); break; case 1: INTEGRATE(RungeKutta4, .0, 30.0, DEFAULT_DT, 2, forced_osc); break;
case 2: INTEGRATE(RungeKutta4, .0, 20.0, DEFAULT_DT, 2, lottka_volt); break; case 2:
case 3: INTEGRATE(RungeKutta4, .0, 20.0, DEFAULT_DT, 2, modified_pp); break; 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: case 4:
y0[0] = .0; // = v(0) y0[0] = .0; // = v(0)
y0[1] = 3.0; // = s(0) y0[1] = 3.0; // = s(0)
...@@ -60,7 +72,7 @@ int main(int argc, char **argv) { ...@@ -60,7 +72,7 @@ int main(int argc, char **argv) {
s_buf[i] = y0[1]; 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(v_buf);
free(s_buf); free(s_buf);
break; break;
......
...@@ -19,7 +19,7 @@ for i in range(len(lines)): ...@@ -19,7 +19,7 @@ for i in range(len(lines)):
figure(1) figure(1)
clf() clf()
for i in range(len(y)): for i in range(len(y)):
plot(x, y[i], 'x') plot(x, y[i], '-')
if len(argv) == 2: if len(argv) == 2:
savefig(argv[1]) savefig(argv[1])
show() show()
...@@ -25,7 +25,8 @@ ...@@ -25,7 +25,8 @@
De opdracht bestaat uit twee delen: het eerste deel betreft het implementeren 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 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} \section{Implementatie}
...@@ -78,7 +79,7 @@ De vergelijking van de drie integratiemethoden kan worden uitgevoerd met: ...@@ -78,7 +79,7 @@ De vergelijking van de drie integratiemethoden kan worden uitgevoerd met:
$ sh ./compare_osc.sh $ sh ./compare_osc.sh
\end{verbatim} \end{verbatim}
In de volgende grafieken staan hiervan de resultaten (blauw
\begin{figure}[H] \begin{figure}[H]
\includegraphics[scale=.5]{osc_euler.pdf} \includegraphics[scale=.5]{osc_euler.pdf}
...@@ -99,15 +100,81 @@ In de volgende grafieken staan hiervan de resultaten (blauw ...@@ -99,15 +100,81 @@ In de volgende grafieken staan hiervan de resultaten (blauw
\subsection{Lottka-Voltera} \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} \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} \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} \subsection{Gilpin's Model}
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment