ModSim: Started on report of assignment 2

modsim/ass2/Makefile

@@ -2,7 +2,7 @@ CC=clang
 CFLAGS=-Wall -Wextra -pedantic -std=c99 -D_GNU_SOURCE -g
 LFLAGS=-lm
 
-all: q1 q2 q3 q4 q5 q7
+all: q1 q2 q3 q4 q5 q7 report.pdf
 
 q%: q%.o
 	$(CC) $(CFLAGS) $(LFLAGS) -o $@ $^
@@ -10,6 +10,10 @@ q%: q%.o
 %.o: %.c
 	$(CC) $(CFLAGS) -o $@ -c $^
 
+%.pdf: %.tex
+	pdflatex $^
+	pdflatex $^
+
 .PHONY: clean
 clean:
 	for i in `seq 7`; do \

modsim/ass2/bisec.h

 #include "func_ptr.h"
 
-#define EPSILON 1e-11
-
-double bisec(func_ptr f, double left, double right, int *steps) {
+double bisec(func_ptr f, double left, double right, double epsilon, unsigned int *steps) {
     int i;
     double mid;
 
-    for( i = 1; fabs(right - left) > 2 * EPSILON; i++ ) {
+    for( i = 1; fabs(right - left) > 2 * epsilon; i++ ) {
         mid = (right + left) / 2;
 
         if( f(left) * f(mid) < 0 )

modsim/ass2/plot.py

+import matplotlib.pyplot as plt
+
+data = [[5 , 1e-1],
+        [8 , 1e-2],
+        [11, 1e-3],
+        [15, 1e-4],
+        [18, 1e-5],
+        [21, 1e-6],
+        [25, 1e-7],
+        [28, 1e-8],
+        [31, 1e-9]]
+
+plt.plot(zip(*data)[1], zip(*data)[0], 'o')
+plt.xscale('log')
+plt.xlabel('epsilon')
+plt.ylabel('steps')
+plt.grid(True)
+plt.savefig('plot.pdf')

modsim/ass2/q1.c

@@ -3,9 +3,8 @@
 #include <math.h>
 #include "func_ptr.h"
 
-#define H 1e-3
-#define TABLE_LINE(func, x) (printf("%-24s%.12f\t%.12f\n", #x, \
-            slope_right(func, x, H), slope_central(func, x, H)))
+#define TABLE_LINE(func, x, h) (printf("%-24s%.12f\t%.12f\n", #x, \
+            slope_right(func, x, h), slope_central(func, x, h)))
 
 double slope_right(func_ptr func, double x, double h) {
     return (func(x + h) - func(x)) / h;
@@ -17,9 +16,9 @@ double slope_central(func_ptr func, double x, double h) {
 
 int main(void) {
     puts("x\t\t\tright\t\tcentral");
-    TABLE_LINE(&sin, M_PI / 3);
-    TABLE_LINE(&sin, 100 * M_PI + M_PI / 3);
-    TABLE_LINE(&sin, 1e12 * M_PI + M_PI / 3);
+    TABLE_LINE(&sin, M_PI / 3, 1e-5);
+    TABLE_LINE(&sin, 100 * M_PI + M_PI / 3, 1e-4);
+    TABLE_LINE(&sin, 1e12 * M_PI + M_PI / 3, 1e-3);
 
     return 0;
 }

modsim/ass2/q2.c

@@ -7,11 +7,21 @@ double func(double x) {
     return x * sin(x) - 1;
 }
 
-int main(void) {
-    int steps;
+#define BISEC(exponent) (bisec(&func, 0, 2, pow(10, -1.0 * exponent), &steps))
 
-    printf("zero point: %.20f\n", bisec(&func, 0, 2, &steps));
-    printf("Steps: %d\n", steps);
+int main(int argc, char *argv[]) {
+    unsigned int steps, i, begin, end;
+
+    if( argc != 3 ) {
+        printf("Usage: %s BEGIN END", argv[0]);
+        return 1;
+    }
+
+    begin = atoi(argv[1]);
+    end = atoi(argv[2]);
+
+    for( i = begin; i <= end; i++ )
+        printf("zero point: %.30f for epsilon = 1e-%d (%d steps)\n", BISEC(i), i, steps);
 
     return 0;
 }

modsim/ass2/q3.c

@@ -2,16 +2,19 @@
 #include <stdio.h>
 #include <math.h>
 #include "bisec.h"
+#include "newton_raphson.h"
+
+#define EPSILON 1e-11
 
 double quad(double x) {
     return x * x - 2;
 }
 
 int main(void) {
-    int steps;
+    unsigned int steps;
 
-    printf("Square root of 2: %.11f\n", bisec(&quad, 1.4, 1.5, &steps));
-    printf("Steps: %d\n", steps);
+    printf("Square root of 2 using bisection method: %.11f (%d steps)\n",
+           bisec(&quad, 1.4, 1.5, EPSILON, &steps), steps);
 
     return 0;
 }

modsim/ass2/report.tex

+\documentclass[10pt,a4paper]{article}
+\usepackage{float,url,graphicx}
+
+\usepackage[dutch]{babel}
+
+\title{Modelleren, Simuleren \& Contin\"ue Wiskunde \\
+       Assignment 2: Differentiation, roots and Integration}
+\author{Tadde\"us Kroes (6054129) \and Sander van Veen (6167969)}
+
+\begin{document}
+
+\maketitle
+
+\section{Differentiatie} % {{{
+\label{sec:Differentiatie}
+
+We kunnen op verschillende manieren de richtingsco\"effici\"ent van een functie
+berekenen. Wij hebben twee manieren getest: het gebruik van een interval rechts
+van een gegeven punt, en een interval om het gegeven punt heen (``right-hand''
+en ``central'' differencing). Voor de berekening gebruiken we respectievelijk de
+volgende functies:
+
+\begin{equation}
+    \label{eq:right}
+    \frac{\sin(x + h) - \sin(x)}{h}
+\end{equation}
+\begin{equation}
+    \label{eq:central}
+    \frac{\sin(x + h) - \sin(x - h)}{2h}
+\end{equation}
+
+Als het interval rechts van een gegeven punt ligt nemen we geen informatie mee
+over de grafiek links van het punt, terwijl we dit wel doen bij centrale
+differentiatie. We verwachten daarom een klein verschil in de resultaten van
+beide methodes. Om de accuratie van de resultaten te testen, berekenen we ook
+handmatig het resultaat met behulp van $\sin(x)' = \cos(x)$. Hiermee berekenen
+we dat $[\sin(\frac{\pi}{3})]' = \cos(\frac{\pi}{3}) = \frac{1}{2}$. De periode
+van $\sin(x)$ id $2\pi$, dus zijn ook beide andere richtingsco\"effici\"enten
+gelijk aan $\frac{1}{2}$ (immers zijn $100\pi$ en $10^{12}\pi$ veelvouden van
+$2\pi$). De resultaten van ons programma staan in de volgende tabel:
+
+\begin{table}[h]
+\begin{tabular}{l|lll}
+    $x$                          & right-hand       & central          & $h$ \\
+    \hline
+    $\frac{\pi}{3}$              & $0.499995669867$ & $0.499999999992$ & $10^{-5}$ \\
+    $100\pi + \frac{\pi}{3}$     & $0.499956697770$ & $0.499999999041$ & $10^{-4}$ \\
+    $10^{12}\pi + \frac{\pi}{3}$ & $0.487956179030$ & $0.488369107296$ & $10^{-3}$ \\
+\end{tabular}
+\end{table}
+
+We hebben ook kleinere waardes voor $h$ gebruikt, maar dan werd het resultaat
+voor $x = 10^{12}\pi + \frac{\pi}{3}$ gelijk aan $0$. Dit is te verklaren door
+het feit dat de grootte $h$ significant genoeg moet zijn om $x$ te veranderen.
+Bij een grote $|x|$ en een kleine $h$ kan het dus voorkomen dat voor het begin
+en het eind van het interval dezelfde waarde wordt gebruikt (namelijk $x$), dit
+levert een richtingsco\"effici\"ent van $0$ op. Een oplossing voor dit probleem
+is om $h$ dynamisch te maken, door $h$ te berekenen aan de hand van $x$,
+bijvoorbeeld door de minst significante mantissa van $x$ te gebruiken. \\ De
+kleinere waarden voor $h$ geven wel betere resultaten voor de kleinere waarden
+van $x$. We zien wel dat right-hand differencing een kleinere
+richtingsco\"effici\"ent geeft dan central differencing. Dit kan worden
+verklaard door het feit dat $\sin(x)$ op deze waarden van $x$ een afnemend stijgende
+lijn is. Aan de rechterkant van $x$ stijgt de lijn dan minder hard dan gemiddeld
+over de linker- en rechterkant.
+
+% }}}
+
+\section{Bisection} % {{{
+\label{sec:Bisection}
+
+We hebben een programma geschreven dat de bisection methode toepast om het
+nulpunt van $f(x)$ te benaderen op het interval $[0, 2]$. De benadering is
+gelijk aan $1.114157142 \dots$, interessanter is echter de relatie tussen de
+nauwkeurigheid en het aantal stappen van de berekening:
+
+\begin{figure}[H]
+\includegraphics[scale=.5]{plot}
+\caption{Het verband tussen de nauwkeurigheid en het aantal stappen van de
+berekening is logaritmisch.}
+\end{figure}
+
+\noindent Uit de grafiek kunnen we concluderen dat de groei van het aantal
+benodigde stappen voor de berekening logaritmisch (met basis 10) afneemt met de
+beoogde nauwkeurigheid.
+
+% }}}
+
+\section{Newton-Raphson} % {{{
+\label{sec:Newton-Raphson}
+
+
+
+% }}}
+
+\end{document}