Commit b8727620 authored by Victor-M-Gomes's avatar Victor-M-Gomes

Commit exo1 module2

parent 0660017b
......@@ -4,7 +4,7 @@
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en">
<head>
<!-- 2020-04-27 lun. 20:01 -->
<!-- 2020-04-28 mar. 11:29 -->
<meta http-equiv="Content-Type" content="text/html;charset=utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1" />
<title>On the computation of pi</title>
......@@ -244,67 +244,98 @@
<body>
<div id="content">
<h1 class="title">On the computation of pi</h1>
<ol class="org-ol">
<li>Asking the math library
<div id="table-of-contents">
<h2>Table of Contents</h2>
<div id="text-table-of-contents">
<ul>
<li><a href="#orgec675a6">1. 1 Asking the math library</a></li>
<li><a href="#org9cc90ed">2. 2 * Buffon&rsquo;s needle</a></li>
<li><a href="#org8060f5c">3. 3. Using a surface fraction argument</a></li>
</ul>
</div>
</div>
<div id="outline-container-orgec675a6" class="outline-2">
<h2 id="orgec675a6"><span class="section-number-2">1</span> 1 Asking the math library</h2>
<div class="outline-text-2" id="text-1">
<p>
My computer tells me that &pi; is approximatively
#begin<sub>src</sub> python :results output :exports both
from math import *
</p>
<div class="org-src-container">
<pre class="src src-python"><span style="color: #51afef;">from</span> math <span style="color: #51afef;">import</span> *
pi
#end<sub>src</sub></li>
</pre>
</div>
</div>
</div>
<li>* Buffon&rsquo;s needle
<div id="outline-container-org9cc90ed" class="outline-2">
<h2 id="org9cc90ed"><span class="section-number-2">2</span> 2 * Buffon&rsquo;s needle</h2>
<div class="outline-text-2" id="text-2">
<p>
Applying the method of <a href="https://en.wikipedia.org/wiki/Buffon%27s_needle_problem">Buffon&rsquo;s needle</a>, we get the <b>approximation</b>
#begin<sub>src</sub> python :results output :exports both
import numpy as np
np.random.seed(seed=42)
N = 10000
x = np.random.uniform(size=N, low=0, high=1)
theta = np.random.uniform(size=N, low)0, high=pi/2)
2/(sum((x+np.sin(theta))&gt;1)/N)
#end<sub>src</sub></li>
<li><p>
Using a surface fraction argument
A method that is easier to understand and does not make use of the <code>sin</code> function is based on the fact that if \(X \approx U(0,1)\) and \(Y \approx U(0,1)\), then \(P[X^2 + Y^2 \leq 1] = \pi/4\) (see <a href="https://en.wikipedia.org/wiki/Monte_Carlo_method">&ldquo;Monte Carlo method on Wikipedia&rdquo;</a>). The following code uses this approach:
#begin<sub>src</sub> python :results output :exports both
import matplotlib.pyplot as plt
</p>
<div class="org-src-container">
<pre class="src src-python"><span style="color: #51afef;">import</span> numpy <span style="color: #51afef;">as</span> np
np.random.seed(seed=<span style="color: #da8548; font-weight: bold;">42</span>)
<span style="color: #dcaeea;">N</span> = <span style="color: #da8548; font-weight: bold;">10000</span>
<span style="color: #dcaeea;">x</span> = np.random.uniform(size=N, low=<span style="color: #da8548; font-weight: bold;">0</span>, high=<span style="color: #da8548; font-weight: bold;">1</span>)
<span style="color: #dcaeea;">theta</span> = np.random.uniform(size=N, low=<span style="color: #da8548; font-weight: bold;">0</span>, high=pi/<span style="color: #da8548; font-weight: bold;">2</span>)
<span style="color: #da8548; font-weight: bold;">2</span>/(<span style="color: #c678dd;">sum</span>((x+np.sin(theta))&gt;<span style="color: #da8548; font-weight: bold;">1</span>)/N)
</pre>
</div>
</div>
</div>
<div id="outline-container-org8060f5c" class="outline-2">
<h2 id="org8060f5c"><span class="section-number-2">3</span> 3. Using a surface fraction argument</h2>
<div class="outline-text-2" id="text-3">
<p>
np.random.seed(seed=42)
N = 1000
x = np.random.uniform(size=N, low=0, high=1)
y = np.random.uniform(size=N, low=0, high=1)
A method that is easier to understand and does not make use of the <code>sin</code> function is based on the fact that if \(X \approx U(0,1)\) and \(Y \approx U(0,1)\), then \(P[X^2 + Y^2 \leq 1] = \pi/4\) (see <a href="https://en.wikipedia.org/wiki/Monte_Carlo_method">&ldquo;Monte Carlo method on Wikipedia&rdquo;</a>). The following code uses this approach:
</p>
<div class="org-src-container">
<pre class="src src-python"><span style="color: #51afef;">import</span> matplotlib.pyplot <span style="color: #51afef;">as</span> plt
<p>
accept = (x*x+y*y) &lt;= 1
reject = np.logical<sub>not</sub>(accept)
</p>
np.random.seed(seed=<span style="color: #da8548; font-weight: bold;">42</span>)
<span style="color: #dcaeea;">N</span> = <span style="color: #da8548; font-weight: bold;">1000</span>
<span style="color: #dcaeea;">x</span> = np.random.uniform(size=N, low=<span style="color: #da8548; font-weight: bold;">0</span>, high=<span style="color: #da8548; font-weight: bold;">1</span>)
<span style="color: #dcaeea;">y</span> = np.random.uniform(size=N, low=<span style="color: #da8548; font-weight: bold;">0</span>, high=<span style="color: #da8548; font-weight: bold;">1</span>)
<p>
fig, ax = plt.subplots(1)
ax.scatter(x[accept], y[accept], c=&rsquo;b&rsquo;, alpha=0.2, edgecolor=None)
ax.scatter(x[reject], y[reject], c=&rsquo;r&rsquo;, alpha=0.2, edgecolor=None)
ax.set<sub>aspect</sub>(&rsquo;equal&rsquo;)
</p>
<span style="color: #dcaeea;">accept</span> = (x*x+y*y) &lt;= <span style="color: #da8548; font-weight: bold;">1</span>
<span style="color: #dcaeea;">reject</span> = np.logical_not(accept)
<span style="color: #dcaeea;">fig</span>, <span style="color: #dcaeea;">ax</span> = plt.subplots(<span style="color: #da8548; font-weight: bold;">1</span>)
ax.scatter(x[accept], y[accept], c=<span style="color: #98be65;">'b'</span>, alpha=<span style="color: #da8548; font-weight: bold;">0.2</span>, edgecolor=<span style="color: #a9a1e1;">None</span>)
ax.scatter(x[reject], y[reject], c=<span style="color: #98be65;">'r'</span>, alpha=<span style="color: #da8548; font-weight: bold;">0.2</span>, edgecolor=<span style="color: #a9a1e1;">None</span>)
ax.set_aspect(<span style="color: #98be65;">'equal'</span>)
plt.savefig(<span style="color: #98be65;">'matplot_lib_filename'</span>)
<span style="color: #51afef;">print</span>(<span style="color: #98be65;">'matplot_lib_filename'</span>)
</pre>
</div>
<pre class="example">
matplot_lib_filename
</pre>
<p>
plt.savefig(matplot<sub>lib</sub><sub>filename</sub>)
print(matplot<sub>lib</sub><sub>filename</sub>)
#end<sub>src</sub>
</p>
<p>
It is then straightforward to obtain a (not really good) approximation to &pi; by counting how many times, on average, \(X^2 + Y^2\) is smaller than 1:
#begin<sub>src</sub> python :results output :exports both
4*np.mean(accept)
#end<sub>src</sub>
</p></li>
</ol>
</p>
<div class="org-src-container">
<pre class="src src-python"><span style="color: #da8548; font-weight: bold;">4</span>*np.mean(accept)
</pre>
</div>
<pre class="example">
3.112
</pre>
</div>
</div>
</div>
<div id="postamble" class="status">
<p class="author">Author: Victor Martins Gomes</p>
<p class="date">Created: 2020-04-27 lun. 20:01</p>
<p class="date">Created: 2020-04-28 mar. 11:29</p>
</div>
</body>
</html>
#+TITLE: On the computation of pi
#+LANGUAGE: en
# #+PROPERTY: header-args :eval never-export
1. Asking the math library
My computer tells me that \pi is approximatively
#begin_src python :results output :exports both
from math import *
pi
#end_src
2. * Buffon's needle
Applying the method of [[https://en.wikipedia.org/wiki/Buffon%27s_needle_problem][Buffon's needle]], we get the *approximation*
#begin_src python :results output :exports both
import numpy as np
np.random.seed(seed=42)
N = 10000
x = np.random.uniform(size=N, low=0, high=1)
theta = np.random.uniform(size=N, low)0, high=pi/2)
2/(sum((x+np.sin(theta))>1)/N)
#end_src
3. Using a surface fraction argument
A method that is easier to understand and does not make use of the =sin= function is based on the fact that if $X \approx U(0,1)$ and $Y \approx U(0,1)$, then $P[X^2 + Y^2 \leq 1] = \pi/4$ (see [[https://en.wikipedia.org/wiki/Monte_Carlo_method]["Monte Carlo method on Wikipedia"]]). The following code uses this approach:
#begin_src python :results output :exports both
import matplotlib.pyplot as plt
np.random.seed(seed=42)
N = 1000
x = np.random.uniform(size=N, low=0, high=1)
y = np.random.uniform(size=N, low=0, high=1)
accept = (x*x+y*y) <= 1
reject = np.logical_not(accept)
fig, ax = plt.subplots(1)
ax.scatter(x[accept], y[accept], c='b', alpha=0.2, edgecolor=None)
ax.scatter(x[reject], y[reject], c='r', alpha=0.2, edgecolor=None)
ax.set_aspect('equal')
plt.savefig(matplot_lib_filename)
print(matplot_lib_filename)
#end_src
It is then straightforward to obtain a (not really good) approximation to \pi by counting how many times, on average, $X^2 + Y^2$ is smaller than 1:
#begin_src python :results output :exports both
4*np.mean(accept)
#end_src
# #+PROPERTY: header-args :session :exports both
* 1 Asking the math library
My computer tells me that \pi is approximatively
#+begin_src python :results value :session *python* :exports both
from math import *
pi
#+end_src
#+RESULTS:
: 3.141592653589793
* 2 * Buffon's needle
Applying the method of [[https://en.wikipedia.org/wiki/Buffon%27s_needle_problem][Buffon's needle]], we get the *approximation*
#+begin_src python :results value :session :*python* :exports both
import numpy as np
np.random.seed(seed=42)
N = 10000
x = np.random.uniform(size=N, low=0, high=1)
theta = np.random.uniform(size=N, low=0, high=pi/2)
2/(sum((x+np.sin(theta))>1)/N)
#+end_src
#+RESULTS:
: 3.128911138923655
* 3. Using a surface fraction argument
A method that is easier to understand and does not make use of the =sin= function is based on the fact that if $X \approx U(0,1)$ and $Y \approx U(0,1)$, then $P[X^2 + Y^2 \leq 1] = \pi/4$ (see [[https://en.wikipedia.org/wiki/Monte_Carlo_method]["Monte Carlo method on Wikipedia"]]). The following code uses this approach:
#+begin_src python :results output file :var matplot_lib_filename="figure_pi_mc2.png" :session *python* :exports both
import matplotlib.pyplot as plt
np.random.seed(seed=42)
N = 1000
x = np.random.uniform(size=N, low=0, high=1)
y = np.random.uniform(size=N, low=0, high=1)
accept = (x*x+y*y) <= 1
reject = np.logical_not(accept)
fig, ax = plt.subplots(1)
ax.scatter(x[accept], y[accept], c='b', alpha=0.2, edgecolor=None)
ax.scatter(x[reject], y[reject], c='r', alpha=0.2, edgecolor=None)
ax.set_aspect('equal')
plt.savefig(matplot_lib_filename)
print(matplot_lib_filename)
#+end_src
#+RESULTS:
[[file:Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/tmp/babel-ekv1Vl/python-2ZIXqL", line 4, in <module>
np.random.seed(seed=42)
NameError: name 'np' is not defined]]
File "<stdin>", line 1, in <module>
File "/tmp/babel-ekv1Vl/python-xlXPMX", line 4, in <module>
np.random.seed(seed=42)
NameError: name 'np' is not defined]]
It is then straightforward to obtain a (not really good) approximation to \pi by counting how many times, on average, $X^2 + Y^2$ is smaller than 1:
#+begin_src python :results output :session *python* :exports both
4*np.mean(accept)
#+end_src
#+RESULTS:
: 3.112
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