-
Notifications
You must be signed in to change notification settings - Fork 1
/
carminati-esercizi-07.html
388 lines (384 loc) · 54.9 KB
/
carminati-esercizi-07.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml" lang="it-IT" xml:lang="it-IT">
<head>
<meta charset="utf-8" />
<meta name="generator" content="pandoc-markdown-css-theme" />
<meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes" />
<meta name="author" content="Leonardo Carminati" />
<meta name="author" content="Maurizio Tomasi" />
<title>Lezione 7: Quadratura numerica</title>
<link rel="stylesheet" href="css/theme.css" />
<link rel="stylesheet" href="css/skylighting-solarized-theme.css" />
<script defer="" src="https://cdn.jsdelivr.net/npm/[email protected]/dist/katex.min.js"></script>
<script>document.addEventListener("DOMContentLoaded", function () {
var mathElements = document.getElementsByClassName("math");
var macros = [];
for (var i = 0; i < mathElements.length; i++) {
var texText = mathElements[i].firstChild;
if (mathElements[i].tagName == "SPAN") {
katex.render(texText.data, mathElements[i], {
displayMode: mathElements[i].classList.contains('display'),
throwOnError: false,
macros: macros,
fleqn: false
});
}}});
</script>
<link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/[email protected]/dist/katex.min.css" />
<script src="./fmtinstall.js"></script>
</head>
<body>
<header>
<h1 class="title">Lezione 7: Quadratura numerica</h1>
<blockquote class="metadata">
<p class="author">
Leonardo CarminatiMaurizio Tomasi
</p>
<p class="date before-toc"><time datetime="A.A. 2024−2025">A.A. 2024−2025</time></p>
</blockquote>
</header>
<nav id="TOC" role="doc-toc">
<strong>Contents</strong><label for="contents">⊕</label>
<input type="checkbox" id="contents">
<ul>
<li><a href="#esercizi-per-oggi" id="toc-esercizi-per-oggi">Esercizi per oggi</a></li>
<li><a href="#esercizio-7.0" id="toc-esercizio-7.0">Esercizio 7.0 - Integrazione con il metodo <em>midpoint</em> a numero di passi fissato</a>
<ul>
<li><a href="#il-metodo-del-mid-point" id="toc-il-metodo-del-mid-point">Il metodo del mid-point</a></li>
<li><a href="#cenni-sullimplementazione" id="toc-cenni-sullimplementazione">Cenni sull’implementazione</a></li>
<li><a href="#e-i-test" id="toc-e-i-test">E i test?</a></li>
</ul></li>
<li><a href="#esercizio-7.1" id="toc-esercizio-7.1">Esercizio 7.1 - Integrazione alla Simpson (da consegnare)</a>
<ul>
<li><a href="#il-metodo-simpson" id="toc-il-metodo-simpson">Il metodo Simpson</a></li>
<li><a href="#test" id="toc-test">Test</a></li>
</ul></li>
<li><a href="#esercizio-7.2" id="toc-esercizio-7.2">Esercizio 7.2 - Integrazione con la formula dei trapezi con precisione fissata (da consegnare)</a>
<ul>
<li><a href="#il-metodo-dei-trapezi" id="toc-il-metodo-dei-trapezi">Il metodo dei trapezi</a></li>
<li><a href="#stima-runtime-dellerrore" id="toc-stima-runtime-dellerrore">Stima runtime dell’errore</a></li>
<li><a href="#cenni-sullimplementazione-1" id="toc-cenni-sullimplementazione-1">Cenni sull’implementazione</a></li>
<li><a href="#test-1" id="toc-test-1">Test</a></li>
</ul></li>
<li><a href="#esercizio-7.3" id="toc-esercizio-7.3">Esercizio 7.3 - Integrazione di una funzione Gaussiana (facoltativo)</a></li>
<li><a href="#errori-comuni" id="toc-errori-comuni">Errori comuni</a></li>
</ul>
</nav>
<main>
<h1 id="esercizi-per-oggi">Esercizi per oggi</h1>
<p>[La pagina con la spiegazione originale degli esercizi si trova qui: <a href="https://labtnds.docs.cern.ch/Lezione7/Lezione7/" class="uri">https://labtnds.docs.cern.ch/Lezione7/Lezione7/</a>.]</p>
<p>In questa lezione implementeremo alcuni algoritmi di <em>quadratura numerica</em>, cioè algoritmi per il calcolo di integrali definiti di funzioni in un intervallo chiuso e limitato. I metodi numerici che studieremo in questa sessione si possono rendere necessari in varie occasioni: quando non sappiamo valutare analiticamente l’integrale in esame, quando la soluzione analitica è molto complicata ed il calcolo numerico è molto più semplice, oppure quando la funzione è conosciuta in un numeri finito di punti.</p>
<h1 id="esercizio-7.0">Esercizio 7.0 - Integrazione con il metodo <em>midpoint</em> a numero di passi fissato</h1>
<p>Implementare un codice per il calcolo della funzione <span class="math inline">x \sin(x)</span> su <span class="math inline">[0, \pi/2]</span> con il metodo midpoint.</p>
<ol type="1">
<li><p>Per prima cosa, costruiamo un programma di test che calcoli l’integrale utilizzando un numero di intervalli fissato e passato da riga di comando.</p></li>
<li><p>Per controllare la precisione ottenuta con un numero di passi fissato, stampiamo una tabella con la differenza tra il risultato numerico ed il valore esatto (ottenuto analiticamente) in funzione del numero di passi (o della lunghezza del passo <span class="math inline">h</span>). In aggiunta alla tabella si può rappresentare l’andamento dell’errore in funzione della lunghezza del passo <span class="math inline">h</span> con un grafico, usando <a href="https://github.com/ziotom78/gplotpp">gplot++</a> oppure un <code>TGraph</code> di ROOT.</p></li>
</ol>
<h2 id="il-metodo-del-mid-point">Il metodo del mid-point</h2>
<p>Ricordiamo che in questo metodo l’approssimazione dell’integrale è definita dalla formula</p>
<p><span class="math display">
\int_a^b f(x)\,\mathrm{d}x = h \cdot \bigl(f(x_0) + f(x_1) + \ldots + f(x_{N - 1})\bigl),
</span> dove <span class="math display">
h = \frac{b - a}N
</span> e <span class="math display">
x_k = a + \left(k + \frac12 \right) h,\qquad k = 0, 1, \ldots, N - 1.
</span></p>
<p>La formula fornisce un’accuratezza dell’integrale di <span class="math inline">O(h^2)</span>. Notate che questo metodo non richiede il calcolo della funzione negli estremi di integrazione.</p>
<h2 id="cenni-sullimplementazione">Cenni sull’implementazione</h2>
<p>L’algoritmo può essere implementato con uno schema del tipo <em>classe madre astratta</em>/<em>classe derivata</em>, che abbiamo già utilizzato nella lezione precedente. Si potrebbe pensare ad una organizzazione di questo tipo:</p>
<ul>
<li><p>Prepariamo una <em>classe astratta</em> <code>Integral</code>, che rappresenta un generico algoritmo di integrazione. Il suo costruttore accetta come parametri gli estremi <span class="math inline">a</span> e <span class="math inline">b</span> di integrazione, e controlla il loro ordinamento: se <span class="math inline">a > b</span>, il segno dell’integrale deve essere scambiato.</p></li>
<li><p>La classe implementa (tra le varie cose) un metodo virtuale puro</p>
<div class="sourceCode" id="cb1"><pre class="sourceCode cpp"><code class="sourceCode cpp"><span id="cb1-1"><a href="#cb1-1" aria-hidden="true" tabindex="-1"></a><span class="kw">virtual</span> <span class="dt">double</span> calculate<span class="op">(</span><span class="dt">int</span> step<span class="op">,</span> <span class="at">const</span> FunzioneBase <span class="op">&</span> f<span class="op">)</span> <span class="op">=</span> <span class="dv">0</span><span class="op">;</span></span></code></pre></div></li>
<li><p>Il metodo <code>calculate</code> è definito come <code>private</code> e non può essere invocato all’esterno della classe. Esso viene invocato dal metodo pubblico (stavolta sì!) <code>integrate</code>, che si preoccupa di verificare gli estremi <code>a</code> e <code>b</code> e li scambia se <code>a > b</code>.</p></li>
<li><p>La suddivisione dell’implementazione in <code>calculate</code> (privato) e <code>integrate</code> (pubblico) permette alle classi derivate di assumere nell’implementazione di <code>calculate</code> sempre che <code>a < b</code>, perché se ne occupa il metodo <code>integrate</code> che è sempre lo stesso (non è <code>virtual</code>).</p></li>
<li><p>Il metodo <em>midpoint</em> viene concretamente implementato in una classe derivata <code>Midpoint</code> sovrascrivendo il metodo privato <code>calculate</code>.</p></li>
</ul>
<div class="sourceCode" id="cb2"><pre class="sourceCode cpp"><code class="sourceCode cpp"><span id="cb2-1"><a href="#cb2-1" aria-hidden="true" tabindex="-1"></a><span class="pp">#pragma once</span></span>
<span id="cb2-2"><a href="#cb2-2" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb2-3"><a href="#cb2-3" aria-hidden="true" tabindex="-1"></a><span class="pp">#include </span><span class="im">"funzioni.h"</span></span>
<span id="cb2-4"><a href="#cb2-4" aria-hidden="true" tabindex="-1"></a><span class="pp">#include </span><span class="im"><cstdlib></span><span class="pp"> </span><span class="co">// Per std::exit()</span></span>
<span id="cb2-5"><a href="#cb2-5" aria-hidden="true" tabindex="-1"></a><span class="pp">#include </span><span class="im"><iostream></span><span class="pp"> </span><span class="co">// Per std::cerr</span></span>
<span id="cb2-6"><a href="#cb2-6" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb2-7"><a href="#cb2-7" aria-hidden="true" tabindex="-1"></a><span class="kw">using</span> <span class="kw">namespace</span> std<span class="op">;</span></span>
<span id="cb2-8"><a href="#cb2-8" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb2-9"><a href="#cb2-9" aria-hidden="true" tabindex="-1"></a><span class="kw">class</span> Integral <span class="op">{</span></span>
<span id="cb2-10"><a href="#cb2-10" aria-hidden="true" tabindex="-1"></a><span class="kw">public</span><span class="op">:</span></span>
<span id="cb2-11"><a href="#cb2-11" aria-hidden="true" tabindex="-1"></a> Integral<span class="op">()</span> <span class="op">:</span> <span class="va">m_a</span><span class="op">{},</span> <span class="va">m_b</span><span class="op">{},</span> <span class="va">m_sign</span><span class="op">{},</span> <span class="va">m_h</span><span class="op">{}</span> <span class="op">{}</span></span>
<span id="cb2-12"><a href="#cb2-12" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb2-13"><a href="#cb2-13" aria-hidden="true" tabindex="-1"></a> <span class="dt">double</span> integrate<span class="op">(</span><span class="dt">double</span> a<span class="op">,</span> <span class="dt">double</span> b<span class="op">,</span> <span class="dt">int</span> nstep<span class="op">,</span> FunzioneBase <span class="op">&</span>f<span class="op">)</span> <span class="op">{</span></span>
<span id="cb2-14"><a href="#cb2-14" aria-hidden="true" tabindex="-1"></a> <span class="co">// Questo metodo fa molto poco: imposta `a` e</span></span>
<span id="cb2-15"><a href="#cb2-15" aria-hidden="true" tabindex="-1"></a> <span class="co">// `b`, verifica che a < b, e poi invoca il</span></span>
<span id="cb2-16"><a href="#cb2-16" aria-hidden="true" tabindex="-1"></a> <span class="co">// metodo virtuale `calculate`, che va ridefinito</span></span>
<span id="cb2-17"><a href="#cb2-17" aria-hidden="true" tabindex="-1"></a> <span class="co">// nelle classi derivate</span></span>
<span id="cb2-18"><a href="#cb2-18" aria-hidden="true" tabindex="-1"></a> checkInterval<span class="op">(</span>a<span class="op">,</span> b<span class="op">);</span></span>
<span id="cb2-19"><a href="#cb2-19" aria-hidden="true" tabindex="-1"></a> <span class="cf">return</span> <span class="va">m_sign</span> <span class="op">*</span> calculate<span class="op">(</span>nstep<span class="op">,</span> f<span class="op">);</span></span>
<span id="cb2-20"><a href="#cb2-20" aria-hidden="true" tabindex="-1"></a> <span class="op">}</span></span>
<span id="cb2-21"><a href="#cb2-21" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb2-22"><a href="#cb2-22" aria-hidden="true" tabindex="-1"></a> <span class="dt">double</span> getA<span class="op">()</span> <span class="at">const</span> <span class="op">{</span> <span class="cf">return</span> <span class="va">m_a</span><span class="op">;</span> <span class="op">}</span></span>
<span id="cb2-23"><a href="#cb2-23" aria-hidden="true" tabindex="-1"></a> <span class="dt">double</span> getB<span class="op">()</span> <span class="at">const</span> <span class="op">{</span> <span class="cf">return</span> <span class="va">m_b</span><span class="op">;</span> <span class="op">}</span></span>
<span id="cb2-24"><a href="#cb2-24" aria-hidden="true" tabindex="-1"></a> <span class="dt">double</span> getSign<span class="op">()</span> <span class="at">const</span> <span class="op">{</span> <span class="cf">return</span> <span class="va">m_sign</span><span class="op">;</span> <span class="op">}</span></span>
<span id="cb2-25"><a href="#cb2-25" aria-hidden="true" tabindex="-1"></a> <span class="dt">double</span> getH<span class="op">()</span> <span class="at">const</span> <span class="op">{</span> <span class="cf">return</span> <span class="va">m_h</span><span class="op">;</span> <span class="op">}</span></span>
<span id="cb2-26"><a href="#cb2-26" aria-hidden="true" tabindex="-1"></a> <span class="dt">void</span> setH<span class="op">(</span><span class="dt">double</span> h<span class="op">)</span> <span class="op">{</span> <span class="va">m_h</span> <span class="op">=</span> h<span class="op">;</span> <span class="op">}</span></span>
<span id="cb2-27"><a href="#cb2-27" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb2-28"><a href="#cb2-28" aria-hidden="true" tabindex="-1"></a><span class="kw">private</span><span class="op">:</span></span>
<span id="cb2-29"><a href="#cb2-29" aria-hidden="true" tabindex="-1"></a> <span class="co">// Questa è la funzione da ridefinire con `override` nelle classi derivate</span></span>
<span id="cb2-30"><a href="#cb2-30" aria-hidden="true" tabindex="-1"></a> <span class="co">// Essa usa come estremi di integrazione m_a ed m_b, ed è *sempre*</span></span>
<span id="cb2-31"><a href="#cb2-31" aria-hidden="true" tabindex="-1"></a> <span class="co">// garantito che m_a < m_b (perché se ne occupa `Integral::integrate`)</span></span>
<span id="cb2-32"><a href="#cb2-32" aria-hidden="true" tabindex="-1"></a> <span class="co">// (Notate che un metodo si può ridefinire nelle classi derivate anche</span></span>
<span id="cb2-33"><a href="#cb2-33" aria-hidden="true" tabindex="-1"></a> <span class="co">// se è `private`!)</span></span>
<span id="cb2-34"><a href="#cb2-34" aria-hidden="true" tabindex="-1"></a> <span class="kw">virtual</span> <span class="dt">double</span> calculate<span class="op">(</span><span class="dt">int</span> nstep<span class="op">,</span> FunzioneBase <span class="op">&</span>f<span class="op">)</span> <span class="op">=</span> <span class="dv">0</span><span class="op">;</span></span>
<span id="cb2-35"><a href="#cb2-35" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb2-36"><a href="#cb2-36" aria-hidden="true" tabindex="-1"></a> <span class="dt">void</span> checkInterval<span class="op">(</span><span class="dt">double</span> a<span class="op">,</span> <span class="dt">double</span> b<span class="op">)</span> <span class="op">{</span></span>
<span id="cb2-37"><a href="#cb2-37" aria-hidden="true" tabindex="-1"></a> <span class="va">m_a</span> <span class="op">=</span> <span class="bu">std::</span>min<span class="op">(</span>a<span class="op">,</span> b<span class="op">);</span></span>
<span id="cb2-38"><a href="#cb2-38" aria-hidden="true" tabindex="-1"></a> <span class="va">m_b</span> <span class="op">=</span> <span class="bu">std::</span>max<span class="op">(</span>a<span class="op">,</span> b<span class="op">);</span></span>
<span id="cb2-39"><a href="#cb2-39" aria-hidden="true" tabindex="-1"></a> <span class="va">m_sign</span> <span class="op">=</span> <span class="op">(</span>a <span class="op"><</span> b<span class="op">)</span> <span class="op">?</span> <span class="dv">1</span> <span class="op">:</span> <span class="op">-</span><span class="dv">1</span><span class="op">;</span></span>
<span id="cb2-40"><a href="#cb2-40" aria-hidden="true" tabindex="-1"></a> <span class="op">}</span></span>
<span id="cb2-41"><a href="#cb2-41" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb2-42"><a href="#cb2-42" aria-hidden="true" tabindex="-1"></a> <span class="dt">double</span> <span class="va">m_a</span><span class="op">,</span> <span class="va">m_b</span><span class="op">;</span></span>
<span id="cb2-43"><a href="#cb2-43" aria-hidden="true" tabindex="-1"></a> <span class="dt">double</span> <span class="va">m_sign</span><span class="op">;</span></span>
<span id="cb2-44"><a href="#cb2-44" aria-hidden="true" tabindex="-1"></a> <span class="dt">double</span> <span class="va">m_h</span><span class="op">;</span></span>
<span id="cb2-45"><a href="#cb2-45" aria-hidden="true" tabindex="-1"></a><span class="op">};</span></span>
<span id="cb2-46"><a href="#cb2-46" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb2-47"><a href="#cb2-47" aria-hidden="true" tabindex="-1"></a><span class="co">// Classe derivata, implementa il metodo mid-point</span></span>
<span id="cb2-48"><a href="#cb2-48" aria-hidden="true" tabindex="-1"></a><span class="kw">class</span> Midpoint <span class="op">:</span> <span class="kw">public</span> Integral <span class="op">{</span></span>
<span id="cb2-49"><a href="#cb2-49" aria-hidden="true" tabindex="-1"></a><span class="kw">public</span><span class="op">:</span></span>
<span id="cb2-50"><a href="#cb2-50" aria-hidden="true" tabindex="-1"></a> Midpoint<span class="op">()</span> <span class="op">:</span> Integral<span class="op">()</span> <span class="op">{}</span></span>
<span id="cb2-51"><a href="#cb2-51" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb2-52"><a href="#cb2-52" aria-hidden="true" tabindex="-1"></a><span class="kw">private</span><span class="op">:</span></span>
<span id="cb2-53"><a href="#cb2-53" aria-hidden="true" tabindex="-1"></a> <span class="dt">double</span> calculate<span class="op">(</span><span class="dt">int</span> nstep<span class="op">,</span> FunzioneBase <span class="op">&</span>f<span class="op">)</span> <span class="kw">override</span> <span class="op">{</span></span>
<span id="cb2-54"><a href="#cb2-54" aria-hidden="true" tabindex="-1"></a> <span class="co">// Ricordare: in quest'implementazione possiamo</span></span>
<span id="cb2-55"><a href="#cb2-55" aria-hidden="true" tabindex="-1"></a> <span class="co">// assumere che m_a < m_b, perché il controllo</span></span>
<span id="cb2-56"><a href="#cb2-56" aria-hidden="true" tabindex="-1"></a> <span class="co">// viene fatto da `Integral::integrate`</span></span>
<span id="cb2-57"><a href="#cb2-57" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb2-58"><a href="#cb2-58" aria-hidden="true" tabindex="-1"></a> <span class="cf">if</span> <span class="op">(</span>nstep <span class="op"><</span> <span class="dv">0</span><span class="op">)</span> <span class="op">{</span></span>
<span id="cb2-59"><a href="#cb2-59" aria-hidden="true" tabindex="-1"></a> cerr <span class="op"><<</span> <span class="st">"Errore, il numero di passi non può essere negativo!</span><span class="sc">\n</span><span class="st">"</span><span class="op">;</span></span>
<span id="cb2-60"><a href="#cb2-60" aria-hidden="true" tabindex="-1"></a> exit<span class="op">(</span><span class="dv">1</span><span class="op">);</span></span>
<span id="cb2-61"><a href="#cb2-61" aria-hidden="true" tabindex="-1"></a> <span class="op">}</span></span>
<span id="cb2-62"><a href="#cb2-62" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb2-63"><a href="#cb2-63" aria-hidden="true" tabindex="-1"></a> <span class="va">m_h</span> <span class="op">=</span> <span class="op">(</span>getB<span class="op">()</span> <span class="op">-</span> getA<span class="op">())</span> <span class="op">/</span> nstep<span class="op">;</span></span>
<span id="cb2-64"><a href="#cb2-64" aria-hidden="true" tabindex="-1"></a> <span class="dt">double</span> sum<span class="op">{};</span></span>
<span id="cb2-65"><a href="#cb2-65" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb2-66"><a href="#cb2-66" aria-hidden="true" tabindex="-1"></a> <span class="cf">for</span> <span class="op">(</span><span class="dt">int</span> i<span class="op">{};</span> i <span class="op"><</span> nstep<span class="op">;</span> i<span class="op">++)</span> <span class="op">{</span></span>
<span id="cb2-67"><a href="#cb2-67" aria-hidden="true" tabindex="-1"></a> sum <span class="op">+=</span> f<span class="op">.</span>Eval<span class="op">(</span><span class="va">m_a</span> <span class="op">+</span> <span class="op">(</span>i <span class="op">+</span> <span class="fl">0.5</span><span class="op">)</span> <span class="op">*</span> <span class="va">m_h</span><span class="op">);</span></span>
<span id="cb2-68"><a href="#cb2-68" aria-hidden="true" tabindex="-1"></a> <span class="op">}</span></span>
<span id="cb2-69"><a href="#cb2-69" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb2-70"><a href="#cb2-70" aria-hidden="true" tabindex="-1"></a> <span class="co">// Non c'è bisogno di moltiplicare per m_sign,</span></span>
<span id="cb2-71"><a href="#cb2-71" aria-hidden="true" tabindex="-1"></a> <span class="co">// questo verrà fatto dalla funzione (pubblica)</span></span>
<span id="cb2-72"><a href="#cb2-72" aria-hidden="true" tabindex="-1"></a> <span class="co">// `Integral::integrate` della classe primitiva</span></span>
<span id="cb2-73"><a href="#cb2-73" aria-hidden="true" tabindex="-1"></a> <span class="cf">return</span> sum <span class="op">*</span> <span class="va">m_h</span><span class="op">;</span></span>
<span id="cb2-74"><a href="#cb2-74" aria-hidden="true" tabindex="-1"></a> <span class="op">}</span></span>
<span id="cb2-75"><a href="#cb2-75" aria-hidden="true" tabindex="-1"></a><span class="op">};</span></span></code></pre></div>
<p>(Attenzione: il codice sopra ha un problema di accessibilità e non compila. Cercate di capire da soli in che modo sistemarlo, e se non riuscite, chiedete al docente o all’assistente).</p>
<p>Notate in che modo il codice implementa il calcolo: il metodo pubblico è <code>Integral::integrate</code>, che <strong>non</strong> è virtuale: esso si preoccupa di invocare <code>Integral::checkInterval</code> (privato) per verificare gli estremi di integrazione, e poi invoca il metodo privato <code>Integral::calculate</code> che fa il conto vero e proprio:</p>
<div class="sourceCode" id="cb3"><pre class="sourceCode cpp"><code class="sourceCode cpp"><span id="cb3-1"><a href="#cb3-1" aria-hidden="true" tabindex="-1"></a><span class="dt">double</span> integrate<span class="op">(</span><span class="dt">double</span> a<span class="op">,</span> <span class="dt">double</span> b<span class="op">,</span> <span class="dt">int</span> nstep<span class="op">,</span> FunzioneBase <span class="op">&</span>f<span class="op">)</span> <span class="op">{</span></span>
<span id="cb3-2"><a href="#cb3-2" aria-hidden="true" tabindex="-1"></a> checkInterval<span class="op">(</span>a<span class="op">,</span> b<span class="op">);</span></span>
<span id="cb3-3"><a href="#cb3-3" aria-hidden="true" tabindex="-1"></a> <span class="cf">return</span> <span class="va">m_sign</span> <span class="op">*</span> calculate<span class="op">(</span>nstep<span class="op">,</span> f<span class="op">);</span></span>
<span id="cb3-4"><a href="#cb3-4" aria-hidden="true" tabindex="-1"></a><span class="op">}</span></span></code></pre></div>
<p>In questo modo il metodo <code>calculate</code> può usare gli estremi restituiti da <code>getA()</code>/<code>getB()</code> senza doversi preoccupare del caso <span class="math inline">a > b</span>. Se vi stupisce che il metodo <code>calculate</code>, pur essendo dichiarato <code>private</code>, possa essere ridefinito nella classe derivata <code>Midpoint</code>, considerate che <code>private</code> indica chi può <strong>chiamare</strong> il metodo (solo la classe <code>Integrate</code> e non le sue derivate), ma non pone restrizioni su chi possa <strong>ridefinirlo</strong>.</p>
<p>Viene ora fornito un codice per verificare il funzionamento di quanto implementato finora, che usa la libreria <code>fmtlib</code>. Come al solito, potete installarla usando lo script <a href="./install_fmt_library"><code>install_fmt_library</code></a>: scaricatelo nella directory dell’esercizio ed eseguitelo con <code>sh install_fmt_library</code>, oppure eseguite direttamente questo comando:</p>
<p><input type="text" value="curl https://ziotom78.github.io/tnds-tomasi-notebooks/install_fmt_library | sh" id="installFmtCommand" readonly="1" size="60"><button onclick='copyFmtInstallationScript("installFmtCommand")'>Copia</button></p>
<p>In alternativa, scaricate questo <a href="./fmtlib.zip">file zip</a> nella directory dell’esercizio e decomprimetelo. Le istruzioni dettagliate sono qui: <a href="index.html#fmtinstall">index.html#fmtinstall</a>.</p>
<p>Ecco il codice di esempio:</p>
<div class="sourceCode" id="cb4"><pre class="sourceCode cpp"><code class="sourceCode cpp"><span id="cb4-1"><a href="#cb4-1" aria-hidden="true" tabindex="-1"></a><span class="pp">#include </span><span class="im">"integral.h"</span></span>
<span id="cb4-2"><a href="#cb4-2" aria-hidden="true" tabindex="-1"></a><span class="pp">#include </span><span class="im">"funzioni.h"</span></span>
<span id="cb4-3"><a href="#cb4-3" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb4-4"><a href="#cb4-4" aria-hidden="true" tabindex="-1"></a><span class="pp">#include </span><span class="im"><cmath></span></span>
<span id="cb4-5"><a href="#cb4-5" aria-hidden="true" tabindex="-1"></a><span class="pp">#include </span><span class="im"><iostream></span></span>
<span id="cb4-6"><a href="#cb4-6" aria-hidden="true" tabindex="-1"></a><span class="pp">#include </span><span class="im"><numbers></span></span>
<span id="cb4-7"><a href="#cb4-7" aria-hidden="true" tabindex="-1"></a><span class="pp">#include </span><span class="im"><string></span></span>
<span id="cb4-8"><a href="#cb4-8" aria-hidden="true" tabindex="-1"></a><span class="pp">#include </span><span class="im">"fmtlib.h"</span></span>
<span id="cb4-9"><a href="#cb4-9" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb4-10"><a href="#cb4-10" aria-hidden="true" tabindex="-1"></a><span class="kw">using</span> <span class="kw">namespace</span> std<span class="op">;</span></span>
<span id="cb4-11"><a href="#cb4-11" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb4-12"><a href="#cb4-12" aria-hidden="true" tabindex="-1"></a><span class="dt">int</span> main <span class="op">(</span><span class="dt">int</span> argc<span class="op">,</span> <span class="dt">char</span><span class="op">*</span> argv<span class="op">[])</span> <span class="op">{</span></span>
<span id="cb4-13"><a href="#cb4-13" aria-hidden="true" tabindex="-1"></a> <span class="cf">if</span> <span class="op">(</span>argc <span class="op">!=</span> <span class="dv">2</span><span class="op">)</span> <span class="op">{</span></span>
<span id="cb4-14"><a href="#cb4-14" aria-hidden="true" tabindex="-1"></a> fmt<span class="op">::</span>println<span class="op">(</span>stderr<span class="op">,</span> <span class="st">"Usage: </span><span class="sc">{}</span><span class="st"> <NSTEP>"</span><span class="op">,</span> argv<span class="op">[</span><span class="dv">0</span><span class="op">]);</span></span>
<span id="cb4-15"><a href="#cb4-15" aria-hidden="true" tabindex="-1"></a> <span class="cf">return</span> <span class="dv">1</span><span class="op">;</span></span>
<span id="cb4-16"><a href="#cb4-16" aria-hidden="true" tabindex="-1"></a> <span class="op">}</span></span>
<span id="cb4-17"><a href="#cb4-17" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb4-18"><a href="#cb4-18" aria-hidden="true" tabindex="-1"></a> <span class="dt">int</span> nstep<span class="op">{</span>stoi<span class="op">(</span>argv<span class="op">[</span><span class="dv">1</span><span class="op">])};</span></span>
<span id="cb4-19"><a href="#cb4-19" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb4-20"><a href="#cb4-20" aria-hidden="true" tabindex="-1"></a> XSinX f<span class="op">{};</span></span>
<span id="cb4-21"><a href="#cb4-21" aria-hidden="true" tabindex="-1"></a> Midpoint myInt<span class="op">{};</span></span>
<span id="cb4-22"><a href="#cb4-22" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb4-23"><a href="#cb4-23" aria-hidden="true" tabindex="-1"></a> <span class="dt">double</span> I<span class="op">{</span>myInt<span class="op">.</span>integrate<span class="op">(</span><span class="dv">0</span><span class="op">,</span> numbers<span class="op">::</span>pi <span class="op">/</span> <span class="dv">2</span><span class="op">,</span> nstep<span class="op">,</span> f<span class="op">)};</span></span>
<span id="cb4-24"><a href="#cb4-24" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb4-25"><a href="#cb4-25" aria-hidden="true" tabindex="-1"></a> fmt<span class="op">::</span>println<span class="op">(</span><span class="st">"Passi: </span><span class="sc">{}</span><span class="st">, I = </span><span class="sc">{}</span><span class="st">"</span><span class="op">,</span> nstep<span class="op">,</span> I<span class="op">);</span></span>
<span id="cb4-26"><a href="#cb4-26" aria-hidden="true" tabindex="-1"></a><span class="op">}</span></span></code></pre></div>
<p>dove abbiamo utilizzato una classe <code>XSinX</code> che eredita dalla classe astratta <code>FunzioneBase</code>.</p>
<p>Per creare i grafici, potete ovviamente usare ROOT, oppure <a href="https://github.com/ziotom78/gplotpp">gplot++</a>. In quest’ultimo caso, scaricate il file <a href="https://raw.githubusercontent.com/ziotom78/gplotpp/master/gplot%2B%2B.h">gplot++.h</a> (facendo click col tasto destro sul link) e scrivete un codice del genere:</p>
<div class="sourceCode" id="cb5"><pre class="sourceCode cpp"><code class="sourceCode cpp"><span id="cb5-1"><a href="#cb5-1" aria-hidden="true" tabindex="-1"></a><span class="bu">std::</span>vector<span class="op"><</span><span class="dt">int</span><span class="op">></span> steps<span class="op">{</span><span class="dv">10</span><span class="op">,</span> <span class="dv">50</span><span class="op">,</span> <span class="dv">100</span><span class="op">,</span> <span class="dv">500</span><span class="op">,</span> <span class="dv">1000</span><span class="op">};</span></span>
<span id="cb5-2"><a href="#cb5-2" aria-hidden="true" tabindex="-1"></a><span class="bu">std::</span>vector<span class="op"><</span><span class="dt">double</span><span class="op">></span> step_sizes<span class="op">(</span>ssize<span class="op">(</span>steps<span class="op">));</span></span>
<span id="cb5-3"><a href="#cb5-3" aria-hidden="true" tabindex="-1"></a><span class="bu">std::</span>vector<span class="op"><</span><span class="dt">double</span><span class="op">></span> errors<span class="op">(</span>ssize<span class="op">(</span>steps<span class="op">));</span></span>
<span id="cb5-4"><a href="#cb5-4" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb5-5"><a href="#cb5-5" aria-hidden="true" tabindex="-1"></a><span class="co">// Calcola gli errori e stampa una tabella usando "fmtlib.h"</span></span>
<span id="cb5-6"><a href="#cb5-6" aria-hidden="true" tabindex="-1"></a><span class="dt">double</span> true_value<span class="op">{</span><span class="dv">1</span><span class="op">};</span></span>
<span id="cb5-7"><a href="#cb5-7" aria-hidden="true" tabindex="-1"></a>fmt<span class="op">::</span>println<span class="op">(</span><span class="st">"Passi Intervallo h Errore"</span><span class="op">);</span></span>
<span id="cb5-8"><a href="#cb5-8" aria-hidden="true" tabindex="-1"></a><span class="cf">for</span> <span class="op">(</span><span class="dt">size_t</span> i<span class="op">{};</span> i <span class="op"><</span> ssize<span class="op">(</span>steps<span class="op">);</span> <span class="op">++</span>i<span class="op">)</span> <span class="op">{</span></span>
<span id="cb5-9"><a href="#cb5-9" aria-hidden="true" tabindex="-1"></a> <span class="dt">double</span> estimated_value<span class="op">{</span>myInt<span class="op">.</span>integrate<span class="op">(</span>steps<span class="op">[</span>i<span class="op">],</span> f<span class="op">)};</span></span>
<span id="cb5-10"><a href="#cb5-10" aria-hidden="true" tabindex="-1"></a> errors<span class="op">[</span>i<span class="op">]</span> <span class="op">=</span> fabs<span class="op">(</span>estimated_value <span class="op">-</span> true_value<span class="op">);</span></span>
<span id="cb5-11"><a href="#cb5-11" aria-hidden="true" tabindex="-1"></a> step_sizes<span class="op">[</span>i<span class="op">]</span> <span class="op">=</span> myInt<span class="op">.</span>GetH<span class="op">();</span></span>
<span id="cb5-12"><a href="#cb5-12" aria-hidden="true" tabindex="-1"></a> fmt<span class="op">::</span>println<span class="op">(</span><span class="st">"</span><span class="sc">{:12d}</span><span class="st"> </span><span class="sc">{:14.8e}</span><span class="st"> </span><span class="sc">{:20.8e}</span><span class="st">"</span><span class="op">,</span> steps<span class="op">[</span>i<span class="op">],</span> step_sizes<span class="op">[</span>i<span class="op">],</span> errors<span class="op">[</span>i<span class="op">]);</span></span>
<span id="cb5-13"><a href="#cb5-13" aria-hidden="true" tabindex="-1"></a><span class="op">}</span></span>
<span id="cb5-14"><a href="#cb5-14" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb5-15"><a href="#cb5-15" aria-hidden="true" tabindex="-1"></a><span class="co">// Crea un plot</span></span>
<span id="cb5-16"><a href="#cb5-16" aria-hidden="true" tabindex="-1"></a>Gnuplot plt<span class="op">{};</span></span>
<span id="cb5-17"><a href="#cb5-17" aria-hidden="true" tabindex="-1"></a><span class="at">const</span> <span class="bu">std::</span>string output_file_name<span class="op">{</span><span class="st">"midpoint-error.png"</span><span class="op">};</span></span>
<span id="cb5-18"><a href="#cb5-18" aria-hidden="true" tabindex="-1"></a>plt<span class="op">.</span>redirect_to_png<span class="op">(</span>output_file_name<span class="op">,</span> <span class="st">"800,600"</span><span class="op">);</span></span>
<span id="cb5-19"><a href="#cb5-19" aria-hidden="true" tabindex="-1"></a>plt<span class="op">.</span>set_logscale<span class="op">(</span>Gnuplot<span class="op">::</span>AxisScale<span class="op">::</span>LOGXY<span class="op">);</span></span>
<span id="cb5-20"><a href="#cb5-20" aria-hidden="true" tabindex="-1"></a>plt<span class="op">.</span>plot<span class="op">(</span>step_sizes<span class="op">,</span> errors<span class="op">);</span></span>
<span id="cb5-21"><a href="#cb5-21" aria-hidden="true" tabindex="-1"></a>plt<span class="op">.</span>set_xlabel<span class="op">(</span><span class="st">"Passo di integrazione h"</span><span class="op">);</span></span>
<span id="cb5-22"><a href="#cb5-22" aria-hidden="true" tabindex="-1"></a>plt<span class="op">.</span>set_ylabel<span class="op">(</span><span class="st">"Errore"</span><span class="op">);</span></span>
<span id="cb5-23"><a href="#cb5-23" aria-hidden="true" tabindex="-1"></a>plt<span class="op">.</span>show<span class="op">();</span></span>
<span id="cb5-24"><a href="#cb5-24" aria-hidden="true" tabindex="-1"></a><span class="co">// È sempre consigliato fornire un messaggio all'utente</span></span>
<span id="cb5-25"><a href="#cb5-25" aria-hidden="true" tabindex="-1"></a><span class="co">// per comunicare che è stato salvato un plot. Includete</span></span>
<span id="cb5-26"><a href="#cb5-26" aria-hidden="true" tabindex="-1"></a><span class="co">// sempre il nome del file nel messaggio!</span></span>
<span id="cb5-27"><a href="#cb5-27" aria-hidden="true" tabindex="-1"></a>fmt<span class="op">::</span>println<span class="op">(</span><span class="st">"Plot saved in '</span><span class="sc">{}</span><span class="st">'"</span><span class="op">,</span> output_file_name<span class="op">);</span></span></code></pre></div>
<p>Con ROOT si scriverebbe invece qualcosa del genere:</p>
<div class="sourceCode" id="cb6"><pre class="sourceCode cpp"><code class="sourceCode cpp"><span id="cb6-1"><a href="#cb6-1" aria-hidden="true" tabindex="-1"></a><span class="bu">std::</span>vector<span class="op"><</span><span class="dt">int</span><span class="op">></span> steps<span class="op">{</span><span class="dv">10</span><span class="op">,</span> <span class="dv">50</span><span class="op">,</span> <span class="dv">100</span><span class="op">,</span> <span class="dv">500</span><span class="op">,</span> <span class="dv">1000</span><span class="op">};</span></span>
<span id="cb6-2"><a href="#cb6-2" aria-hidden="true" tabindex="-1"></a>TGraph <span class="va">g_errore</span><span class="op">{};</span></span>
<span id="cb6-3"><a href="#cb6-3" aria-hidden="true" tabindex="-1"></a><span class="dt">double</span> true_value<span class="op">{</span><span class="dv">1</span><span class="op">};</span></span>
<span id="cb6-4"><a href="#cb6-4" aria-hidden="true" tabindex="-1"></a>fmt<span class="op">::</span>println<span class="op">(</span><span class="st">"Passi Errore"</span><span class="op">)</span></span>
<span id="cb6-5"><a href="#cb6-5" aria-hidden="true" tabindex="-1"></a><span class="cf">for</span> <span class="op">(</span><span class="dt">int</span> i<span class="op">{};</span> i <span class="op"><</span> size<span class="op">(</span>steps<span class="op">);</span> i<span class="op">++)</span> <span class="op">{</span></span>
<span id="cb6-6"><a href="#cb6-6" aria-hidden="true" tabindex="-1"></a> <span class="dt">double</span> estimated_value<span class="op">{</span>myInt<span class="op">.</span>integrate<span class="op">(</span>steps<span class="op">[</span>i<span class="op">],</span> f<span class="op">)};</span></span>
<span id="cb6-7"><a href="#cb6-7" aria-hidden="true" tabindex="-1"></a> <span class="dt">double</span> err<span class="op">{</span>fabs<span class="op">(</span>estimated_value <span class="op">-</span> true_value<span class="op">)};</span></span>
<span id="cb6-8"><a href="#cb6-8" aria-hidden="true" tabindex="-1"></a> fmt<span class="op">::</span>println<span class="op">(</span><span class="st">"</span><span class="sc">{:12d}</span><span class="st"> </span><span class="sc">{:20.8e}</span><span class="st">"</span><span class="op">,</span> steps<span class="op">[</span>i<span class="op">],</span> err<span class="op">);</span></span>
<span id="cb6-9"><a href="#cb6-9" aria-hidden="true" tabindex="-1"></a> <span class="va">g_errore</span><span class="op">.</span>SetPoint<span class="op">(</span>i<span class="op">,</span> myInt<span class="op">.</span>GetH<span class="op">(),</span> err<span class="op">);</span></span>
<span id="cb6-10"><a href="#cb6-10" aria-hidden="true" tabindex="-1"></a><span class="op">}</span></span></code></pre></div>
<h2 id="e-i-test">E i test?</h2>
<p>Da oggi non fornirò più direttamente il codice delle funzioni <code>test_??()</code>. Fornirò però una implementazione <em>completa</em> di ogni esercizio scritta nel linguaggio Julia. È un linguaggio molto semplice da leggere, ma il fatto che sia diverso dal C++ impedirà di “barare” facendo copia-e-incolla dei miei esempi.</p>
<p>Per ogni esercizio, il codice Julia esegue una serie di calcoli che vi forniscono i risultati attesi. Sarà semplice quindi per voi creare da soli le funzioni <code>test_??()</code>.</p>
<p>Il notebook Julia per la lezione corrente si trova all’indirizzo <a href="https://ziotom78.github.io/tnds-notebooks/lezione07/" class="uri">https://ziotom78.github.io/tnds-notebooks/lezione07/</a>.</p>
<h1 id="esercizio-7.1">Esercizio 7.1 - Integrazione alla Simpson (da consegnare)</h1>
<p>Implementare l’integrazione con il metodo di Simpson con un numero di passi fissato. Si può utilizzare lo stesso schema dell’esercizio precedente costruendo una classe derivata <code>Simpson</code>. Come nell’<a href="#esercizio-7.0">Esercizio 7.0</a>, stampare una tabella (o costruire un grafico) con la precisione raggiunta in funzione del numero di passi.</p>
<h2 id="il-metodo-simpson">Il metodo Simpson</h2>
<p>Nel metodo di integrazione alla Simpson, la funzione integranda è approssimata, negli intervalli <span class="math inline">[x_k, x_{k + 2}]</span> (dove <span class="math inline">k</span> è un intero pari e <span class="math inline">x_k = a + k h</span>) con un polinomio di secondo grado i cui nodi sono nei punti <span class="math inline">\bigl(x_k, f(x_k)\bigr)</span>, <span class="math inline">\bigl(x_{k + 1}, f(x_{k + 1})\bigr)</span>, <span class="math inline">\bigl(x_{k + 2}, f(x_{k + 2})\bigr)</span>. Esso fornisce una valutazione dell’integrale con una precisione pari a <span class="math inline">h^4</span>. La sua applicazione richiede che il numero di passi sia pari e la formula che approssima l’integrale è</p>
<p><span class="math display">
\int_a^b f(x)\,\mathrm{d}x = h \cdot \left(
\frac13 f(x_o) + \frac43 f(x_1) + \frac23 f(x_2) + \frac43 f(x_3) + \ldots
+ \frac23 f(x_{N - 2}) + \frac43 f(x_{N - 1}) + \frac13 f(x_N)\bigl)
\right),
</span></p>
<p>dove vale che</p>
<p><span class="math display">
h = \frac{b - a}N, \quad x_k = a + k h \qquad k = 0, 1, \ldots N.
</span></p>
<h2 id="test">Test</h2>
<p>Come scritto per l’esercizio 7.0, il notebook Julia per la lezione corrente all’indirizzo <a href="https://ziotom78.github.io/tnds-notebooks/lezione07/" class="uri">https://ziotom78.github.io/tnds-notebooks/lezione07/</a> contiene una serie di risultati che potete usare per scrivere i vostri test.</p>
<h1 id="esercizio-7.2">Esercizio 7.2 - Integrazione con la formula dei trapezi con precisione fissata (da consegnare)</h1>
<p>Concludiamo l’esercitazione implementando l’integrazione della funzione <span class="math inline">x \sin x</span> su <span class="math inline">[0, \pi]</span> con il metodo dei trapezi. In quest’ultimo esercizio implementiamo un algoritmo di integrazione numerica a precisione fissata invece che a numero di passi fissato. Negli esercizi precedenti il calcolo a numero di passi fissato non ci dà alcuna indicazione sulla qualità del risultato: integrare con 10 passi è sufficiente? L’idea che vogliamo sviluppare è che l’utente fornisca una precisione desiderata e l’algoritmo sia in grado di aumentare automaticamente il numero di passi fino a raggiungere la precisione richiesta sul valore dell’integrale. L’algoritmo dovrà accettare in input il valore della precisione e raddoppiare il numero di passi finché l’errore (stimato runtime, si veda sotto) non diventa inferiore alla precisione impostata.</p>
<ol type="1">
<li><p>Come nei casi precedenti si può costruire una classe dedicata per l’implementazione del metodo dei trapezi.</p></li>
<li><p>Possiamo implementare come nei casi precedenti un metodo privato</p>
<div class="sourceCode" id="cb7"><pre class="sourceCode cpp"><code class="sourceCode cpp"><span id="cb7-1"><a href="#cb7-1" aria-hidden="true" tabindex="-1"></a><span class="dt">double</span> calculate<span class="op">(</span><span class="dt">int</span> nstep<span class="op">,</span> <span class="at">const</span> FunzioneBase <span class="op">&</span> f<span class="op">)</span> <span class="kw">override</span> <span class="op">{</span> <span class="op">...</span> <span class="op">}</span></span></code></pre></div></li>
<li><p>Il punto che ci interessa maggiormente è realizzare un metodo dei trapezi che lavori a precisione fissata. Noi implementeremo il calcolo a precisione fissata solo per il metodo dei trapezi, quindi non toccheremo la classe <code>Integral</code>; di conseguenza, implementeremo solo un metodo <code>integrate_prec</code> in <code>Trapezi</code>:</p>
<div class="sourceCode" id="cb8"><pre class="sourceCode cpp"><code class="sourceCode cpp"><span id="cb8-1"><a href="#cb8-1" aria-hidden="true" tabindex="-1"></a><span class="dt">double</span> integrate_prec<span class="op">(</span><span class="dt">double</span> a<span class="op">,</span> <span class="dt">double</span> b<span class="op">,</span> <span class="dt">double</span> prec<span class="op">,</span> <span class="at">const</span> FunzioneBase <span class="op">&</span> f<span class="op">)</span> <span class="op">{</span> <span class="op">...</span></span>
<span id="cb8-2"><a href="#cb8-2" aria-hidden="true" tabindex="-1"></a><span class="op">}</span></span></code></pre></div>
<p>Ovviamente, è possibile definire una coppia di funzioni <code>integrate_prec</code> (pubblica) e <code>calculate_prec</code> (privata) già nella classe <code>Integral</code> e poi implementare <code>calculate_prec</code> in ogni classe derivata. Bisogna però fare attenzione, perché il calcolo dell’errore dipende dall’ordine dell’algoritmo, quindi sarebbe necessario definire un nuovo metodo astratto <code>calculate_error()</code> che accetti l’integrale calcolato con passo <span class="math inline">h</span> e con passo <span class="math inline">h/2</span> per determinare l’errore corretto, e invocarlo poi in <code>integrate_prec</code> (il metodo pubblico). Noi, come già detto, non faremo nulla di ciò: ci basta implementare l’algoritmo a precisione fissata per il solo metodo dei trapezi. La classe <code>Trapezoids</code> avrà quindi un’implementazione simile:</p>
<div class="sourceCode" id="cb9"><pre class="sourceCode cpp"><code class="sourceCode cpp"><span id="cb9-1"><a href="#cb9-1" aria-hidden="true" tabindex="-1"></a><span class="kw">class</span> Trapezoids <span class="op">:</span> <span class="kw">public</span> Integral <span class="op">{</span></span>
<span id="cb9-2"><a href="#cb9-2" aria-hidden="true" tabindex="-1"></a><span class="kw">private</span><span class="op">:</span></span>
<span id="cb9-3"><a href="#cb9-3" aria-hidden="true" tabindex="-1"></a> <span class="co">// Metodo dei trapezi con NUMERO DI PASSI fissato (come negli esercizi 7.0 e 7.1)</span></span>
<span id="cb9-4"><a href="#cb9-4" aria-hidden="true" tabindex="-1"></a> <span class="dt">double</span> calculate<span class="op">(</span><span class="dt">int</span> nstep<span class="op">,</span> <span class="at">const</span> FunzioneBase <span class="op">&</span> f<span class="op">)</span> <span class="kw">override</span> <span class="op">{</span> <span class="op">...</span> <span class="op">}</span></span>
<span id="cb9-5"><a href="#cb9-5" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb9-6"><a href="#cb9-6" aria-hidden="true" tabindex="-1"></a><span class="kw">public</span><span class="op">:</span></span>
<span id="cb9-7"><a href="#cb9-7" aria-hidden="true" tabindex="-1"></a> <span class="co">// Metodo dei trapezi con PRECISIONE FISSATA (nuovo!). Questo metodo si</span></span>
<span id="cb9-8"><a href="#cb9-8" aria-hidden="true" tabindex="-1"></a> <span class="co">// implementa direttamente come pubblico</span></span>
<span id="cb9-9"><a href="#cb9-9" aria-hidden="true" tabindex="-1"></a> <span class="dt">double</span> integrate_prec<span class="op">(</span><span class="dt">double</span> a<span class="op">,</span> <span class="dt">double</span> b<span class="op">,</span> <span class="dt">double</span> prec<span class="op">,</span> <span class="at">const</span> FunzioneBase <span class="op">&</span> f<span class="op">)</span> <span class="op">{</span></span>
<span id="cb9-10"><a href="#cb9-10" aria-hidden="true" tabindex="-1"></a> <span class="op">...</span></span>
<span id="cb9-11"><a href="#cb9-11" aria-hidden="true" tabindex="-1"></a> <span class="op">}</span></span>
<span id="cb9-12"><a href="#cb9-12" aria-hidden="true" tabindex="-1"></a><span class="op">};</span></span></code></pre></div></li>
</ol>
<p>Un algoritmo a precisione fissata si può implementare anche per il metodo midpoint e per il metodo di Simpson ma nel caso dei trapezoidi si presta ad una implementazione particolarmente efficiente (si veda la sezione sotto).</p>
<h2 id="il-metodo-dei-trapezi">Il metodo dei trapezi</h2>
<p>Il metodo dei trapezi rappresenta la prima formula chiusa di Newton-Cotes e si basa su una interpolazione polinomiale di ordine 1 sui due punti che delimitano l’intervallo. La formula che approssima l’integrale è</p>
<p><span class="math display">
\int_a^b f(x)\,\mathrm{d}x \approx \left(\frac12 f(x_0) + f(x_1) + \dots + f(x_{N - 1}) + \frac12 f(x_N)\right)h,
</span></p>
<p>dove</p>
<p><span class="math display">
h = \frac{b - a}N
</span></p>
<p>e</p>
<p><span class="math display">
x_k = a + k h, \qquad k = 0, 1, \ldots N.
</span></p>
<h2 id="stima-runtime-dellerrore">Stima runtime dell’errore</h2>
<p>Come possiamo stimare l’errore che stiamo commettendo nel calcolo di un integrale se non conosciamo il valore vero <span class="math inline">I</span> dell’integrale? Partendo dalla conoscenza dell’andamento teorico dell’errore in funzione del passo <span class="math inline">h</span>, possiamo trovare un modo semplice per la stima dell’errore che stiamo commettendo. Nel caso della regola dei trapezi, l’adamento dell’errore è <span class="math inline">\epsilon = k h^2</span>. Calcolando l’integrale <span class="math inline">I_N</span> utilizzando un passo <span class="math inline">h</span> e successivamente l’intgrale <span class="math inline">I_{2N}</span> con un passo <span class="math inline">h/2</span>, possiamo scrivere il seguente sistema di equazioni: <span class="math display">
\begin{cases}
I - I_N = k h^2,\\
I - I_{2N} = k\left(\frac{h}2\right)^2.
\end{cases}
</span></p>
<p>Sottraendo per esempio la prima equazione alla seconda, non è difficile ricavare che una stima dell’errore pari a <span class="math display">
\epsilon(N) = \frac43 \left|I_{2N} - I_N\right|.
</span></p>
<h2 id="cenni-sullimplementazione-1">Cenni sull’implementazione</h2>
<p>La condizione di uscita dell’algoritmo è basata sul confronto del risultato di due passaggi successivi. Se la differenza della stima dell’integrale tra due passaggi consecutivi è più piccola della precisione richiesta allora l’algoritmo si ferma. Questo poiché l’integrale calcolato con una partizione più fine è sempre una stima migliore dell’approssimazione con la partizione originaria. L’algoritmo dei trapezoidi si presta ad una implementazione ottimizzata descritta qui di seguito. Nel costruttore calcoliamo la prima approssimazione:</p>
<p><span class="math display">
\begin{cases}
\Sigma_0 &= \frac{f(a) + f(b)}2,\\
I_0 &= \Sigma_0 \times (b - a).
\end{cases}
</span></p>
<p>Al primo passo dell’algoritmo suddividiamo l’intervallo in due:</p>
<p><span class="math display">
\begin{cases}
\Sigma_1 = \Sigma_0 + f\left(x_{11}\right),\\
I_1 = \Sigma_1 \times \frac{b - a}2.
\end{cases}
</span></p>
<p>Al secondo passaggio otteniamo</p>
<p><span class="math display">
\begin{cases}
\Sigma_2 = \Sigma_1 + f\left(x_{21}\right) + f\left(x_{22}\right),\\
I_2 = \Sigma_2 \times \frac{b - a}4.
\end{cases}
</span></p>
<p>Al terzo passaggio otteniamo</p>
<p><span class="math display">
\begin{cases}
\Sigma_3 = \Sigma_1 + f\left(x_{31}\right) + f\left(x_{32}\right) + f\left(x_{33}\right) + f\left(x_{34}\right),\\
I_3 = \Sigma_3 \times \frac{b - a}8.
\end{cases}
</span> e così via secondolo schema in figura, con <span class="math inline">L = b - a</span>:</p>
<p><img src="https://labtnds.docs.cern.ch/Lezione7/intervallo.png" /></p>
<p>I valori dell’ultima approssimazione dell’integrale e dell’ultima somma calcolata sono memorizzati all’interno dell’oggetto. In questo modo, se viene richiesto di ricalcolare l’integrale, non è necessario ricominciare da capo.</p>
<h2 id="test-1">Test</h2>
<p>Come scritto per l’esercizio 7.0, il notebook Julia per la lezione corrente all’indirizzo <a href="https://ziotom78.github.io/tnds-notebooks/lezione07/" class="uri">https://ziotom78.github.io/tnds-notebooks/lezione07/</a> contiene una serie di risultati che potete usare per scrivere i vostri test.</p>
<h1 id="esercizio-7.3">Esercizio 7.3 - Integrazione di una funzione Gaussiana (facoltativo)</h1>
<p>Come esercizio facoltativo vi proponiamo il calcolo di un integrale non risolubile analiticamente. Assumendo che una misura abbia un valore vero <span class="math inline">\mu</span> e una deviazione standard <span class="math inline">\sigma</span>, calcoliamo la probabilità che una misura cada entro <span class="math inline">\pm t\sigma</span> dal valore vero per <span class="math inline">t</span> che va da 0 a 5. La probabilità è ovviamente data dalla <a href="https://en.wikipedia.org/wiki/Normal_distribution">distribuzione Gaussiana</a>:</p>
<p><span class="math display">
f(x) = \frac1{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x - \mu)^2}{2\sigma^2}\right).
</span></p>
<ol>
<li><p>Costruite una classe <code>Gaussian</code> generica. I parametri <span class="math inline">\mu</span> e <span class="math inline">\sigma</span> della funzione <span class="math inline">f(x)</span> possono essere passati nel costruttore.</p></li>
<li><p>Integrate la funzione Gaussiana con il metodo dei trapezi a precisione fissata in un intervallo <span class="math inline">[\mu - t\sigma, \mu + t\sigma]</span> con <span class="math inline">t</span> variabile da 0 a 5.</p></li>
<li><p>Producete un grafico dell’integrale ottenuto in funzione del numero di <span class="math inline">\sigma</span>. Vi torna?</p></li>
</ol>
<p>Il risultato non dovrebbe sorprendere:</p>
<p><img src="https://labtnds.docs.cern.ch/Lezione7/probability.png" /></p>
<h1 id="errori-comuni">Errori comuni</h1>
<p>Come di consueto, elenco alcuni errori molto comuni che ho trovato negli anni passati correggendo gli esercizi che gli studenti hanno consegnato all’esame:</p>
<ul>
<li><p>Molte volte gli studenti usano la funzione <code>abs</code> anziché <code>fabs</code> nel determinare l’ampiezza dell’intervallo di integrazione. Attenzione! Se <code>fabs</code> restituisce sempre un numero floating-point, la funzione <code>abs</code> lo fa <strong>solo</strong> se si include <code><cmath></code>, altrimenti restituisce un intero: se quindi non includete <code><cmath></code> e l’intervallo di integrazione <code>b-a</code> è inferiore a 1, avrete che <code>abs(b - a) == 0</code>!</p></li>
<li><p>Alcuni studenti calcolano le somme dei termini degli integrali saltando l’ultimo punto a destra dell’intervallo (e sottostimando quindi l’integrale di <span class="math inline">h \cdot \pi/2 \cdot \sin \pi/2</span>).</p></li>
<li><p>Attenzione ai coefficienti nella formula di Simpson, perché a volte gli studenti scambiano di posto il 4 con il 2!</p></li>
<li><p>Non confondete il significato di “numero di passi” quando calcolate l’errore di un metodo di integrazione: se per calcolare l’errore dovete stimare l’integrale con passo <code>h</code> e con passo <code>h/2</code>, l’errore che ottenete si riferisce al passo <code>h</code>, non al passo <code>h/2</code>! (Va detto però che quest’errore è più frequente negli esami scritti che negli esercizi.)</p></li>
</ul>
<p>Il notebook all’indirizzo <a href="https://ziotom78.github.io/tnds-notebooks/lezione07/" class="uri">https://ziotom78.github.io/tnds-notebooks/lezione07/</a> fornisce una lunga serie di test: se li implementate tutti, avete il 99.9% di essere sicuri che la vostra implementazione sia corretta!</p>
</main>
<script>
;(function() {
// Non-essential if user has JavaScript off. Just makes checkboxes look nicer.
var selector = '.task-list > li > input[type="checkbox"]';
var checkboxes = document.querySelectorAll(selector);
Array.from(checkboxes).forEach((checkbox) => {
var wasChecked = checkbox.checked;
checkbox.disabled = false;
checkbox.addEventListener('click', (ev) => {ev.target.checked = wasChecked});
});
})();
</script>
</body>
</html>