ModSim: finished assignment 2.4

modsim/ass2/bisection.c

@@ -4,7 +4,7 @@ double bisec(func_ptr f, double left, double right,
         double epsilon, unsigned int *steps, unsigned int max_steps) {
     double mid = left, fmid;
 
-    for( *steps = 0; (fabs(right - left) > 2 * epsilon)
+    for( ; (fabs(right - left) > 2 * epsilon)
             && (!max_steps || (*steps < max_steps)); (*steps)++ ) {
         mid = (right + left) / 2;
 

modsim/ass2/integration.py

+import matplotlib.pyplot as plt
+from numpy import arange, zeros
+
+t = arange(-3.0, 3.0, 0.01)
+s = t**3 - 3*t - 2
+
+plt.plot(t, s)
+plt.grid(True)
+plt.xlabel('x')
+plt.ylabel('f(x)')
+
+plt.savefig('integration.pdf')
+

modsim/ass2/newton_raphson.c

@@ -4,7 +4,7 @@ double newton_raphson(func_ptr f, func_ptr df, double x_start,
         double epsilon, unsigned int *steps, unsigned int max_steps) {
     double x = x_start, last_x = x - 1;
 
-    for( *steps = 0; fabs(x - last_x) >= epsilon; (*steps)++ ) {
+    for( ; fabs(x - last_x) >= epsilon; (*steps)++ ) {
         if( *steps == max_steps )
             return NAN;
 

modsim/ass2/q3.c

@@ -26,15 +26,18 @@ int main(int argc, char *argv[]) {
     end = atoi(argv[2]);
 
     for( i = begin; i <= end; i++) {
+        steps = 0;
         epsilon = pow(10, -1.0 * i);
         tmp= bisec(&f, 1, 2, epsilon, &steps, 100000);
         printf("Sqrt(2) using bisection:      %.20f (%d steps; epsilon=%.0e)\n",
             tmp, steps, epsilon);
 
+        steps = 0;
         tmp = newton_raphson(&f, &df, 1.4, epsilon, &steps, 100000);
         printf("Sqrt(2) using Newton-Raphson: %.20f (%d steps; epsilon=%.0e)\n",
             tmp, steps, epsilon);
 
+        steps = 0;
         tmp = regula_falsi(&f, 1, 2, epsilon, &steps, 100000);
         printf("Sqrt(2) using Regula Falsi:   %.20f (%d steps; epsilon=%.0e)\n",
             tmp, steps, epsilon);

modsim/ass2/q4.c

@@ -37,6 +37,7 @@ double df3(double x) {
 }
 
 #define PRINT_ZERO(f, df, x_start) { \
+    steps = start_steps; \
     if( !isnan(root = newton_raphson(&f, &df, x_start, EPSILON, &steps, MAX_STEPS)) ) \
         printf(#f ": %.11f (%d steps from start point %.11f)\n", root, steps, (double)x_start); \
     else \
@@ -44,13 +45,14 @@ double df3(double x) {
 }
 
 #define PRINT_ZERO_RANGE(f, df, method, lnbd, ubnd) { \
+    start_steps = 0; \
     x_start = method(&f, lnbd, ubnd, EPSILON, &start_steps, max_start_steps); \
     printf("Did first %d steps using " #method "\n", start_steps); \
     PRINT_ZERO(f, df, x_start); \
 }
 
 int main(int argc, char *argv[]) {
-    unsigned int steps, start_steps, max_start_steps;
+    unsigned int steps, start_steps = 0, max_start_steps;
     double root, x_start;
 
     if( argc != 2 ) {
@@ -67,11 +69,13 @@ int main(int argc, char *argv[]) {
 
     printf("\nNow using Bisection to determine the first %d steps.\n\n", max_start_steps);
 
+    PRINT_ZERO_RANGE(f2, df2, bisec, -10, 0);
     PRINT_ZERO_RANGE(f2, df2, bisec, 0, 10);
     PRINT_ZERO_RANGE(f3, df3, bisec, 0, 10);
 
     printf("\nNow using Regula Falsi to determine the first %d steps.\n\n", max_start_steps);
 
+    PRINT_ZERO_RANGE(f2, df2, regula_falsi, -10, 0);
     PRINT_ZERO_RANGE(f2, df2, regula_falsi, 0, 10);
     PRINT_ZERO_RANGE(f3, df3, regula_falsi, 0, 10);
 

modsim/ass2/regula_falsi.c

@@ -5,11 +5,13 @@ double regula_falsi(func_ptr f, double s, double t, double e,
     int side = 0;
     double r, fr, fs = f(s), ft = f(t);
 
-    for( *steps = 0; *steps < m; (*steps)++ ) {
+    for( ; *steps < m; (*steps)++ ) {
         r = (fs * t - ft * s) / (fs - ft);
 
-        if( fabs(t - s) < e * fabs(t + s) )
+        if( fabs(t - s) < e * fabs(t + s) ) {
+            (*steps)++;
             break;
+        }
 
         fr = f(r);
 
@@ -29,9 +31,11 @@ double regula_falsi(func_ptr f, double s, double t, double e,
                 ft /= 2;
 
             side = 1;
-        } else
+        } else {
+            (*steps)++;
             break;
         }
+    }
 
     return r;
 }

modsim/ass2/report.tex

@@ -127,12 +127,98 @@ aantal stappen is afgerond op een geheel getal.
 
 % }}}
 
+\pagebreak
+
 \section{Newton-Raphson} % {{{
 \label{sec:Newton-Raphson}
 
-% TODO: alleen voor snijlijnen en niet voor raaklijnen op de x-as (bisection,
-% regula falsi). Starting values: plotten. En voeg output van q4 toe en
-% vergelijk het aantal benodigde stappen.
+We hebben een programma geschreven dat de nulpunten van de volgende drie
+functies bepaald:
+
+\begin{itemize}
+    \item $f_1(x) = x^2 - x + 2$
+    \item $f_2(x) = x^3 - 3x - 2$
+    \item $f_3(x) = (x^2 - 1)(x - 4)$
+\end{itemize}
+
+Als parameter voor het programma dient een geheel positief getal op te worden
+gegeven. Dat getal wordt gebruikt voor de bisection en regula falsi methode om
+het totaal aantal stappen om de nulpunten te vinden te reduceren. Hieronder is
+een voorbeeld weergave van het programma met een maximum van \texttt{5} stappen
+voor de bisection en regula falsi methode:
+
+\begin{verbatim}
+$ ./q4 5
+f1: could not find a root after 100000 steps
+f2: -1.00000000011 (39 steps from start point -10.00000000000)
+f2: 2.00000000000 (10 steps from start point 10.00000000000)
+f3: 4.00000000000 (2 steps from start point 0.00000000000)
+
+Now using Bisection to determine the first 5 steps.
+
+Did first 1 steps using bisec
+f2: -1.00000000010 (38 steps from start point -5.00000000000)
+Did first 5 steps using bisec
+f2: 2.00000000000 (10 steps from start point 2.18750000000)
+Did first 5 steps using bisec
+f3: 4.00000000000 (9 steps from start point 4.06250000000)
+
+Now using Regula Falsi to determine the first 5 steps.
+
+Did first 5 steps using regula_falsi
+f2: -0.99999999985 (37 steps from start point 0.39763149048)
+Did first 5 steps using regula_falsi
+f2: -0.99999999992 (38 steps from start point 0.37773413585)
+Did first 5 steps using regula_falsi
+f3: 4.00000000000 (30 steps from start point 1.02650759729)
+\end{verbatim}
+
+De functie $f_1$ heeft geen nulpunten en na 100000 stappen zal onze
+implementatie van de Newton-Raphson methode daarom ook aangeven dat er geen
+nulpunten zijn gevonden.
+
+Het aantal benodigde stappen voor het benaderen van de nulpunten van $f_2$ is
+bij de Newton-Raphson methode gelijk aan $39 + 10 = 49$. Als we de bisection
+methode toepassen, kost het benaderen van de nulpunten van $f_2$ een totaal $10
++ 38$ stappen. Regula falsi in combinatie met de Newton-Raphson methode heeft
+voor de functies $f_2$ en $f_3$ veel meer stappen nodig (resp. $37 + 38 = 75$ en
+$30$ stappen).
+
+Op het moment dat we het maximum stappen van de regula falsi methode verhogen
+naar \texttt{10}, dan ziet het resultaat er heel anders uit:
+
+\begin{verbatim}
+$ ./q4 10
+f1: could not find a root after 100000 steps
+
+ .. snip ..
+
+Now using Regula Falsi to determine the first 10 steps.
+
+Did first 7 steps using regula_falsi
+f2: 2.00000000000 (13 steps from start point 2.62300283682)
+Did first 10 steps using regula_falsi
+f2: 2.00000000000 (14 steps from start point 1.99601959065)
+Did first 10 steps using regula_falsi
+f3: 4.00000000000 (14 steps from start point 4.00986504483)
+\end{verbatim}
+
+In dit geval is het aantal benodigde stappen voor regula falsi in combinatie met
+de Newton-Raphson methode gereduceerd tot slechts 14 stappen voor $f_3$. Merk op
+dat de regula falsi methode voor $f_2$ het tweede nulpunt ($x = -1$) niet kan
+vinden.
+
+\begin{figure}[H]
+\centering
+\includegraphics[width=12cm]{integration}
+\caption{Plot van $f(x) = x^3 - 3x - 2$ voor $x$ in $[-3, 3]$.}
+\end{figure}
+
+De bisection en regula falsi methodes werken goed voor snijlijnen en niet voor
+raaklijnen met de x-as. De functie $f_2$ is in bovenstaande figuur geplot en
+uit het figuur valt goed op te maken dat er een raaklijn op de x-as is voor $x =
+-1$. De output van ons programma laat zien dat deze raaklijn op $x = -1$ niet
+wordt gevonden, aangezien het nooit de x-as snijdt.
 
 % }}}