Numerical integration#
Introduction#
In this chapter, for a given function \(f(x)\) we will compute the definite integral:
\[\int_a^b\,f(x)\,\textrm{d}x\]
where \(a\) and \(b\) are the limits of integration, and \(f(x)\) are the function values obtained either from a table of values or from an analytical function.
In numerical integration, the integral is estimated by \(I\), and the following holds:
\[\int_a^b\,f(x)\,\textrm{d}x= I + E,\]
where \(E\) is the error of the integral estimate.
We will compute the numerical integral based on a discrete sum:
\[I=\sum_{i=0}^{m}A_i\,f(x_i),\]
where \(A_i\) are the weights, \(x_i\) are the nodes on the interval \([a, b]\), and \(m+1\) is the number of nodes.
We will look at two different approaches to numerical integration:
the Newton-Cotes approach, which is based on equidistant nodes (a constant integration step), and
the Gaussian integration approach, where the nodes are placed so as to achieve exactness for polynomials.
Motivational example#
For numerical integration, we will use a concrete example to guide us:
\[\int_1^2 x\,\sin(x)\,\textrm{d}x\]
Let us prepare the nodes. Let the basic step be \(h=0.25\); in this case we have four subintervals and five nodes, at step \(2h\) there are three nodes, and at step \(4h\) only two (the endpoints):
import numpy as np
xg, hg = np.linspace(1, 2, 100, retstep=True) # dense points (for plotting)
x2v, h2v = np.linspace(1, 2, 2, retstep=True) # step h2v = 1 (2 nodes)
x3v, h3v = np.linspace(1, 2, 3, retstep=True) # step h3v = 0.5 (3 nodes)
x4v, h4v = np.linspace(1, 2, 4, retstep=True) # step h4v = 0.33.. (4 nodes)
x5v, h5v = np.linspace(1, 2, 5, retstep=True) # step h5v = 0.25 (5 nodes)
Let us also prepare the function values:
yg = xg * np.sin(xg)
y2v = x2v * np.sin(x2v)
y3v = x3v * np.sin(x3v)
y4v = x4v * np.sin(x4v)
y5v = x5v * np.sin(x5v)
Let us prepare a plot of the data:
import matplotlib.pyplot as plt
from matplotlib import rc # we import this so that the fonts in the figure are latex-compatible
#rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
rc('text', usetex=True)
%matplotlib inline
def fig_motivation():
plt.fill_between(xg, yg, alpha=0.25, facecolor='r')
plt.annotate(r'$\int_1^2\,x\,\sin(x)\,\textrm{d}x$', (1.3, 0.5), fontsize=22)
plt.plot(xg, yg, lw=3, alpha=0.5, label=r'$x\,\sin(x)$')
plt.plot(x2v, y2v, 's', alpha=0.5, label=f'$h={h2v}$', markersize=14)
plt.plot(x3v, y3v, 'o', alpha=0.5, label=f'$h={h3v}$', markersize=10)
plt.legend(loc=(1.01, 0))
plt.ylim(0, 2)
plt.show()
Let us display the data:
fig_motivation()
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:250, in TexManager._run_checked_subprocess(cls, command, tex, cwd)
249 try:
--> 250 report = subprocess.check_output(
251 command, cwd=cwd if cwd is not None else cls._texcache,
252 stderr=subprocess.STDOUT)
253 except FileNotFoundError as exc:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:466, in check_output(timeout, *popenargs, **kwargs)
464 kwargs['input'] = empty
--> 466 return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
467 **kwargs).stdout
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:548, in run(input, capture_output, timeout, check, *popenargs, **kwargs)
546 kwargs['stderr'] = PIPE
--> 548 with Popen(*popenargs, **kwargs) as process:
549 try:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:1026, in Popen.__init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, user, group, extra_groups, encoding, errors, text, umask, pipesize, process_group)
1023 self.stderr = io.TextIOWrapper(self.stderr,
1024 encoding=encoding, errors=errors)
-> 1026 self._execute_child(args, executable, preexec_fn, close_fds,
1027 pass_fds, cwd, env,
1028 startupinfo, creationflags, shell,
1029 p2cread, p2cwrite,
1030 c2pread, c2pwrite,
1031 errread, errwrite,
1032 restore_signals,
1033 gid, gids, uid, umask,
1034 start_new_session, process_group)
1035 except:
1036 # Cleanup if the child failed starting.
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:1955, in Popen._execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, gid, gids, uid, umask, start_new_session, process_group)
1954 if err_filename is not None:
-> 1955 raise child_exception_type(errno_num, err_msg, err_filename)
1956 else:
FileNotFoundError: [Errno 2] No such file or directory: 'latex'
The above exception was the direct cause of the following exception:
RuntimeError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backend_bases.py:2157, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
2154 # we do this instead of `self.figure.draw_without_rendering`
2155 # so that we can inject the orientation
2156 with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2157 self.figure.draw(renderer)
2158 if bbox_inches:
2159 if bbox_inches == "tight":
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:94, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
92 @wraps(draw)
93 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 94 result = draw(artist, renderer, *args, **kwargs)
95 if renderer._rasterizing:
96 renderer.stop_rasterizing()
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/figure.py:3264, in Figure.draw(self, renderer)
3261 # ValueError can occur when resizing a window.
3263 self.patch.draw(renderer)
-> 3264 mimage._draw_list_compositing_images(
3265 renderer, self, artists, self.suppressComposite)
3267 renderer.close_group('figure')
3268 finally:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axes/_base.py:3226, in _AxesBase.draw(self, renderer)
3223 if artists_rasterized:
3224 _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3226 mimage._draw_list_compositing_images(
3227 renderer, self, artists, self.get_figure(root=True).suppressComposite)
3229 renderer.close_group('axes')
3230 self.stale = False
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1405, in Axis.draw(self, renderer)
1402 renderer.open_group(__name__, gid=self.get_gid())
1404 ticks_to_draw = self._update_ticks()
-> 1405 tlb1, tlb2 = self._get_ticklabel_bboxes(ticks_to_draw, renderer)
1407 for tick in ticks_to_draw:
1408 tick.draw(renderer)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1332, in Axis._get_ticklabel_bboxes(self, ticks, renderer)
1330 if renderer is None:
1331 renderer = self.get_figure(root=True)._get_renderer()
-> 1332 return ([tick.label1.get_window_extent(renderer)
1333 for tick in ticks if tick.label1.get_visible()],
1334 [tick.label2.get_window_extent(renderer)
1335 for tick in ticks if tick.label2.get_visible()])
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1332, in <listcomp>(.0)
1330 if renderer is None:
1331 renderer = self.get_figure(root=True)._get_renderer()
-> 1332 return ([tick.label1.get_window_extent(renderer)
1333 for tick in ticks if tick.label1.get_visible()],
1334 [tick.label2.get_window_extent(renderer)
1335 for tick in ticks if tick.label2.get_visible()])
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:969, in Text.get_window_extent(self, renderer, dpi)
964 raise RuntimeError(
965 "Cannot get window extent of text w/o renderer. You likely "
966 "want to call 'figure.draw_without_rendering()' first.")
968 with cbook._setattr_cm(fig, dpi=dpi):
--> 969 bbox, info, descent = self._get_layout(self._renderer)
970 x, y = self.get_unitless_position()
971 x, y = self.get_transform().transform((x, y))
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:373, in Text._get_layout(self, renderer)
370 ys = []
372 # Full vertical extent of font, including ascenders and descenders:
--> 373 _, lp_h, lp_d = _get_text_metrics_with_cache(
374 renderer, "lp", self._fontproperties,
375 ismath="TeX" if self.get_usetex() else False,
376 dpi=self.get_figure(root=True).dpi)
377 min_dy = (lp_h - lp_d) * self._linespacing
379 for i, line in enumerate(lines):
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:69, in _get_text_metrics_with_cache(renderer, text, fontprop, ismath, dpi)
66 """Call ``renderer.get_text_width_height_descent``, caching the results."""
67 # Cached based on a copy of fontprop so that later in-place mutations of
68 # the passed-in argument do not mess up the cache.
---> 69 return _get_text_metrics_with_cache_impl(
70 weakref.ref(renderer), text, fontprop.copy(), ismath, dpi)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:77, in _get_text_metrics_with_cache_impl(renderer_ref, text, fontprop, ismath, dpi)
73 @functools.lru_cache(4096)
74 def _get_text_metrics_with_cache_impl(
75 renderer_ref, text, fontprop, ismath, dpi):
76 # dpi is unused, but participates in cache invalidation (via the renderer).
---> 77 return renderer_ref().get_text_width_height_descent(text, fontprop, ismath)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backends/backend_agg.py:211, in RendererAgg.get_text_width_height_descent(self, s, prop, ismath)
209 _api.check_in_list(["TeX", True, False], ismath=ismath)
210 if ismath == "TeX":
--> 211 return super().get_text_width_height_descent(s, prop, ismath)
213 if ismath:
214 ox, oy, width, height, descent, font_image = \
215 self.mathtext_parser.parse(s, self.dpi, prop)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backend_bases.py:566, in RendererBase.get_text_width_height_descent(self, s, prop, ismath)
562 fontsize = prop.get_size_in_points()
564 if ismath == 'TeX':
565 # todo: handle properties
--> 566 return self.get_texmanager().get_text_width_height_descent(
567 s, fontsize, renderer=self)
569 dpi = self.points_to_pixels(72)
570 if ismath:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:364, in TexManager.get_text_width_height_descent(cls, tex, fontsize, renderer)
362 if tex.strip() == '':
363 return 0, 0, 0
--> 364 dvifile = cls.make_dvi(tex, fontsize)
365 dpi_fraction = renderer.points_to_pixels(1.) if renderer else 1
366 with dviread.Dvi(dvifile, 72 * dpi_fraction) as dvi:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:293, in TexManager.make_dvi(cls, tex, fontsize)
290 with TemporaryDirectory(dir=dvifile.parent) as tmpdir:
291 Path(tmpdir, "file.tex").write_text(
292 cls._get_tex_source(tex, fontsize), encoding='utf-8')
--> 293 cls._run_checked_subprocess(
294 ["latex", "-interaction=nonstopmode", "-halt-on-error",
295 "-no-shell-escape", "file.tex"], tex, cwd=tmpdir)
296 Path(tmpdir, "file.dvi").replace(dvifile)
297 # Also move the tex source to the main cache directory, but
298 # only for backcompat.
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:254, in TexManager._run_checked_subprocess(cls, command, tex, cwd)
250 report = subprocess.check_output(
251 command, cwd=cwd if cwd is not None else cls._texcache,
252 stderr=subprocess.STDOUT)
253 except FileNotFoundError as exc:
--> 254 raise RuntimeError(
255 f'Failed to process string with tex because {command[0]} '
256 'could not be found') from exc
257 except subprocess.CalledProcessError as exc:
258 raise RuntimeError(
259 '{prog} was not able to process the following string:\n'
260 '{tex!r}\n\n'
(...) 267 exc=exc.output.decode('utf-8', 'backslashreplace'))
268 ) from None
RuntimeError: Failed to process string with tex because latex could not be found
<Figure size 640x480 with 1 Axes>
Let us compute the exact result analytically:
import sympy as sym
sym.init_printing()
x = sym.symbols('x')
I_exact = float(sym.integrate(x*sym.sin(x), (x, 1, 2)).evalf())
I_exact
\[\displaystyle 1.44042242098021\]
Newton-Cotes approach#
The figure shows a general case in which the distance between the nodes \(x_i\) is equal to \(h\) (this is an equidistant partition).
In this chapter we will first look at the trapezoidal and the composite trapezoidal rule, and later the Simpson’s and Romberg methods will follow.
Trapezoidal rule#
The trapezoidal rule interpolates the values of the integrand \(f(x)\) on the (sub)interval \([x_0, x_1]\) with a linear function. For two nodes, this means that we approximate the area under the graph of the function as:
\[I_{\textrm{trapezno}}=\sum_{i=0}^{n=1}A_i\,f(x_i)=\frac{h}{2}\cdot\left(f(x_0)+f(x_1)\right).\]
And the weights are:
\[A_0 = A_1 = \frac{1}{2}\,h.\]
Numerical implementation#
The numerical implementation is:
def trapezoidal(y, h):
"""
Trapezoidal rule of integration.
:param y: function values
:param h: integration step
"""
return (y[0] + y[-1])*h/2
Numerical example#
In our concrete case this means that we weight the first and last function value by \(h/2\). In our example \(h=1\):
I_trapezoidal = trapezoidal(y2v, h=h2v)
I_trapezoidal
\[\displaystyle 1.33003291922963\]
Let us prepare a figure:
def fig_trapezoidal():
plt.fill_between(x2v, y2v, alpha=0.25, facecolor='r')
plt.vlines(x2v, 0, y2v, color='r', linestyles='dashed', lw=1)
plt.annotate('$I_{\\textrm{trapezno}}$', (1.4, 0.5), fontsize=22)
plt.annotate('Error', fontsize=20, xy=(1.5, 1.4), xytext=(1.6, 1.8),
arrowprops=dict(facecolor='gray', shrink=0.05))
plt.plot(xg, yg, lw=3, alpha=0.5, label='$x\\,\\sin(x)$')
plt.plot(x2v, y2v, 'o', alpha=0.5, label=f'$h={h2v}$')
plt.legend(loc=(1.01, 0))
plt.ylim(0, 2)
plt.show()
We display:
fig_trapezoidal()
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:250, in TexManager._run_checked_subprocess(cls, command, tex, cwd)
249 try:
--> 250 report = subprocess.check_output(
251 command, cwd=cwd if cwd is not None else cls._texcache,
252 stderr=subprocess.STDOUT)
253 except FileNotFoundError as exc:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:466, in check_output(timeout, *popenargs, **kwargs)
464 kwargs['input'] = empty
--> 466 return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
467 **kwargs).stdout
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:548, in run(input, capture_output, timeout, check, *popenargs, **kwargs)
546 kwargs['stderr'] = PIPE
--> 548 with Popen(*popenargs, **kwargs) as process:
549 try:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:1026, in Popen.__init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, user, group, extra_groups, encoding, errors, text, umask, pipesize, process_group)
1023 self.stderr = io.TextIOWrapper(self.stderr,
1024 encoding=encoding, errors=errors)
-> 1026 self._execute_child(args, executable, preexec_fn, close_fds,
1027 pass_fds, cwd, env,
1028 startupinfo, creationflags, shell,
1029 p2cread, p2cwrite,
1030 c2pread, c2pwrite,
1031 errread, errwrite,
1032 restore_signals,
1033 gid, gids, uid, umask,
1034 start_new_session, process_group)
1035 except:
1036 # Cleanup if the child failed starting.
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:1955, in Popen._execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, gid, gids, uid, umask, start_new_session, process_group)
1954 if err_filename is not None:
-> 1955 raise child_exception_type(errno_num, err_msg, err_filename)
1956 else:
FileNotFoundError: [Errno 2] No such file or directory: 'latex'
The above exception was the direct cause of the following exception:
RuntimeError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backend_bases.py:2157, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
2154 # we do this instead of `self.figure.draw_without_rendering`
2155 # so that we can inject the orientation
2156 with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2157 self.figure.draw(renderer)
2158 if bbox_inches:
2159 if bbox_inches == "tight":
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:94, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
92 @wraps(draw)
93 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 94 result = draw(artist, renderer, *args, **kwargs)
95 if renderer._rasterizing:
96 renderer.stop_rasterizing()
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/figure.py:3264, in Figure.draw(self, renderer)
3261 # ValueError can occur when resizing a window.
3263 self.patch.draw(renderer)
-> 3264 mimage._draw_list_compositing_images(
3265 renderer, self, artists, self.suppressComposite)
3267 renderer.close_group('figure')
3268 finally:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axes/_base.py:3226, in _AxesBase.draw(self, renderer)
3223 if artists_rasterized:
3224 _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3226 mimage._draw_list_compositing_images(
3227 renderer, self, artists, self.get_figure(root=True).suppressComposite)
3229 renderer.close_group('axes')
3230 self.stale = False
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1405, in Axis.draw(self, renderer)
1402 renderer.open_group(__name__, gid=self.get_gid())
1404 ticks_to_draw = self._update_ticks()
-> 1405 tlb1, tlb2 = self._get_ticklabel_bboxes(ticks_to_draw, renderer)
1407 for tick in ticks_to_draw:
1408 tick.draw(renderer)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1332, in Axis._get_ticklabel_bboxes(self, ticks, renderer)
1330 if renderer is None:
1331 renderer = self.get_figure(root=True)._get_renderer()
-> 1332 return ([tick.label1.get_window_extent(renderer)
1333 for tick in ticks if tick.label1.get_visible()],
1334 [tick.label2.get_window_extent(renderer)
1335 for tick in ticks if tick.label2.get_visible()])
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1332, in <listcomp>(.0)
1330 if renderer is None:
1331 renderer = self.get_figure(root=True)._get_renderer()
-> 1332 return ([tick.label1.get_window_extent(renderer)
1333 for tick in ticks if tick.label1.get_visible()],
1334 [tick.label2.get_window_extent(renderer)
1335 for tick in ticks if tick.label2.get_visible()])
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:969, in Text.get_window_extent(self, renderer, dpi)
964 raise RuntimeError(
965 "Cannot get window extent of text w/o renderer. You likely "
966 "want to call 'figure.draw_without_rendering()' first.")
968 with cbook._setattr_cm(fig, dpi=dpi):
--> 969 bbox, info, descent = self._get_layout(self._renderer)
970 x, y = self.get_unitless_position()
971 x, y = self.get_transform().transform((x, y))
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:373, in Text._get_layout(self, renderer)
370 ys = []
372 # Full vertical extent of font, including ascenders and descenders:
--> 373 _, lp_h, lp_d = _get_text_metrics_with_cache(
374 renderer, "lp", self._fontproperties,
375 ismath="TeX" if self.get_usetex() else False,
376 dpi=self.get_figure(root=True).dpi)
377 min_dy = (lp_h - lp_d) * self._linespacing
379 for i, line in enumerate(lines):
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:69, in _get_text_metrics_with_cache(renderer, text, fontprop, ismath, dpi)
66 """Call ``renderer.get_text_width_height_descent``, caching the results."""
67 # Cached based on a copy of fontprop so that later in-place mutations of
68 # the passed-in argument do not mess up the cache.
---> 69 return _get_text_metrics_with_cache_impl(
70 weakref.ref(renderer), text, fontprop.copy(), ismath, dpi)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:77, in _get_text_metrics_with_cache_impl(renderer_ref, text, fontprop, ismath, dpi)
73 @functools.lru_cache(4096)
74 def _get_text_metrics_with_cache_impl(
75 renderer_ref, text, fontprop, ismath, dpi):
76 # dpi is unused, but participates in cache invalidation (via the renderer).
---> 77 return renderer_ref().get_text_width_height_descent(text, fontprop, ismath)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backends/backend_agg.py:211, in RendererAgg.get_text_width_height_descent(self, s, prop, ismath)
209 _api.check_in_list(["TeX", True, False], ismath=ismath)
210 if ismath == "TeX":
--> 211 return super().get_text_width_height_descent(s, prop, ismath)
213 if ismath:
214 ox, oy, width, height, descent, font_image = \
215 self.mathtext_parser.parse(s, self.dpi, prop)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backend_bases.py:566, in RendererBase.get_text_width_height_descent(self, s, prop, ismath)
562 fontsize = prop.get_size_in_points()
564 if ismath == 'TeX':
565 # todo: handle properties
--> 566 return self.get_texmanager().get_text_width_height_descent(
567 s, fontsize, renderer=self)
569 dpi = self.points_to_pixels(72)
570 if ismath:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:364, in TexManager.get_text_width_height_descent(cls, tex, fontsize, renderer)
362 if tex.strip() == '':
363 return 0, 0, 0
--> 364 dvifile = cls.make_dvi(tex, fontsize)
365 dpi_fraction = renderer.points_to_pixels(1.) if renderer else 1
366 with dviread.Dvi(dvifile, 72 * dpi_fraction) as dvi:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:293, in TexManager.make_dvi(cls, tex, fontsize)
290 with TemporaryDirectory(dir=dvifile.parent) as tmpdir:
291 Path(tmpdir, "file.tex").write_text(
292 cls._get_tex_source(tex, fontsize), encoding='utf-8')
--> 293 cls._run_checked_subprocess(
294 ["latex", "-interaction=nonstopmode", "-halt-on-error",
295 "-no-shell-escape", "file.tex"], tex, cwd=tmpdir)
296 Path(tmpdir, "file.dvi").replace(dvifile)
297 # Also move the tex source to the main cache directory, but
298 # only for backcompat.
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:254, in TexManager._run_checked_subprocess(cls, command, tex, cwd)
250 report = subprocess.check_output(
251 command, cwd=cwd if cwd is not None else cls._texcache,
252 stderr=subprocess.STDOUT)
253 except FileNotFoundError as exc:
--> 254 raise RuntimeError(
255 f'Failed to process string with tex because {command[0]} '
256 'could not be found') from exc
257 except subprocess.CalledProcessError as exc:
258 raise RuntimeError(
259 '{prog} was not able to process the following string:\n'
260 '{tex!r}\n\n'
(...) 267 exc=exc.output.decode('utf-8', 'backslashreplace'))
268 ) from None
RuntimeError: Failed to process string with tex because latex could not be found
<Figure size 640x480 with 1 Axes>
Error of the trapezoidal rule#
The difference between the analytical value of the integral and the numerical approximation \(I\) is the error of the method:
\[E = \int_a^bf(x)\,d x-I,\]
If the function \(f(x)\) is at least twice differentiable, one can derive (see, e.g., the reference: Burden, Faires, Burden: Numerical Analysis 10th Ed) an error estimate that applies only to the trapezoidal approximation over the entire integration interval:
\[E_{\textrm{trapezno}}=-\frac{h^3}{12}f''(\xi),\]
where \(h=b-a\) and \(\xi\) is an unknown value on the interval \([a,b]\).
Composite trapezoidal rule#
If we divide the interval \([a, b]\) into \(n\) subintervals and apply the trapezoidal rule of integration on each of them, we speak of the composite trapezoidal rule.
In this case, for each subinterval \(i\) we apply the trapezoidal rule, and thus for the bounds of the subinterval \(x_i\) and \(x_{i+1}\) we use the weights \(A_i=A_{i+i}=h/2\). Since the interior nodes are doubled, it follows that:
\[A_0=A_{n}=\frac{h}{2}\quad\textrm{and for the other nodes:}\quad A_i=h.\]
Here we have assumed subintervals of equal width:
\[h=\frac{x_{n}-x_0}{n}\]
It therefore follows that:
\[I_{\textrm{trapezno sest}}=\sum_{i=0}^{n}A_i\,f(x_i)=\left(\frac{y_0}{2} + y_1+y_2+\cdots+y_{n-1}+\frac{y_{n}}{2}\right)\,h.\]
Numerical implementation#
The numerical implementation is:
def trapezoidal_comp(y, h):
"""
Composite trapezoidal rule of integration.
:param y: function values
:param h: integration step
"""
return (np.sum(y) - 0.5*y[0] - 0.5*y[-1])*h
Numerical example#
Above we already prepared the data for two subintervals (three nodes):
x3v
array([1. , 1.5, 2. ])
h3v
\[\displaystyle 0.5\]
Let us compute the estimate of the integral with the composite trapezoidal rule:
I_trapezoidal_comp = trapezoidal_comp(y3v, h=h3v)
I_trapezoidal_comp
\[\displaystyle 1.41313769956786\]
Let us prepare a figure:
def fig_trapezoidal_comp():
plt.fill_between(x3v, y3v, alpha=0.25, facecolor='r')
plt.vlines(x3v, 0, y3v, color='r', linestyles='dashed', lw=1)
plt.annotate('$I_{\\textrm{trapezno sestavljeno}}$', (1.2, 0.5), fontsize=22)
plt.annotate('Error', fontsize=20, xy=(1.75, 1.68), xytext=(1.4, 1.8),
arrowprops=dict(facecolor='gray', shrink=0.05))
plt.annotate('Error', fontsize=20, xy=(1.2, 1.1), xytext=(1., 1.4),
arrowprops=dict(facecolor='gray', shrink=0.05))
plt.plot(xg, yg, lw=3, alpha=0.5, label='$x\\,\\sin(x)$')
plt.plot(x3v, y3v, 'o', alpha=0.5, label=f'$h={h3v}$')
plt.legend(loc=(1.01, 0))
plt.ylim(0, 2)
plt.show()
We display:
fig_trapezoidal_comp()
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:250, in TexManager._run_checked_subprocess(cls, command, tex, cwd)
249 try:
--> 250 report = subprocess.check_output(
251 command, cwd=cwd if cwd is not None else cls._texcache,
252 stderr=subprocess.STDOUT)
253 except FileNotFoundError as exc:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:466, in check_output(timeout, *popenargs, **kwargs)
464 kwargs['input'] = empty
--> 466 return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
467 **kwargs).stdout
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:548, in run(input, capture_output, timeout, check, *popenargs, **kwargs)
546 kwargs['stderr'] = PIPE
--> 548 with Popen(*popenargs, **kwargs) as process:
549 try:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:1026, in Popen.__init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, user, group, extra_groups, encoding, errors, text, umask, pipesize, process_group)
1023 self.stderr = io.TextIOWrapper(self.stderr,
1024 encoding=encoding, errors=errors)
-> 1026 self._execute_child(args, executable, preexec_fn, close_fds,
1027 pass_fds, cwd, env,
1028 startupinfo, creationflags, shell,
1029 p2cread, p2cwrite,
1030 c2pread, c2pwrite,
1031 errread, errwrite,
1032 restore_signals,
1033 gid, gids, uid, umask,
1034 start_new_session, process_group)
1035 except:
1036 # Cleanup if the child failed starting.
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:1955, in Popen._execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, gid, gids, uid, umask, start_new_session, process_group)
1954 if err_filename is not None:
-> 1955 raise child_exception_type(errno_num, err_msg, err_filename)
1956 else:
FileNotFoundError: [Errno 2] No such file or directory: 'latex'
The above exception was the direct cause of the following exception:
RuntimeError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backend_bases.py:2157, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
2154 # we do this instead of `self.figure.draw_without_rendering`
2155 # so that we can inject the orientation
2156 with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2157 self.figure.draw(renderer)
2158 if bbox_inches:
2159 if bbox_inches == "tight":
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:94, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
92 @wraps(draw)
93 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 94 result = draw(artist, renderer, *args, **kwargs)
95 if renderer._rasterizing:
96 renderer.stop_rasterizing()
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/figure.py:3264, in Figure.draw(self, renderer)
3261 # ValueError can occur when resizing a window.
3263 self.patch.draw(renderer)
-> 3264 mimage._draw_list_compositing_images(
3265 renderer, self, artists, self.suppressComposite)
3267 renderer.close_group('figure')
3268 finally:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axes/_base.py:3226, in _AxesBase.draw(self, renderer)
3223 if artists_rasterized:
3224 _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3226 mimage._draw_list_compositing_images(
3227 renderer, self, artists, self.get_figure(root=True).suppressComposite)
3229 renderer.close_group('axes')
3230 self.stale = False
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1405, in Axis.draw(self, renderer)
1402 renderer.open_group(__name__, gid=self.get_gid())
1404 ticks_to_draw = self._update_ticks()
-> 1405 tlb1, tlb2 = self._get_ticklabel_bboxes(ticks_to_draw, renderer)
1407 for tick in ticks_to_draw:
1408 tick.draw(renderer)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1332, in Axis._get_ticklabel_bboxes(self, ticks, renderer)
1330 if renderer is None:
1331 renderer = self.get_figure(root=True)._get_renderer()
-> 1332 return ([tick.label1.get_window_extent(renderer)
1333 for tick in ticks if tick.label1.get_visible()],
1334 [tick.label2.get_window_extent(renderer)
1335 for tick in ticks if tick.label2.get_visible()])
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1332, in <listcomp>(.0)
1330 if renderer is None:
1331 renderer = self.get_figure(root=True)._get_renderer()
-> 1332 return ([tick.label1.get_window_extent(renderer)
1333 for tick in ticks if tick.label1.get_visible()],
1334 [tick.label2.get_window_extent(renderer)
1335 for tick in ticks if tick.label2.get_visible()])
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:969, in Text.get_window_extent(self, renderer, dpi)
964 raise RuntimeError(
965 "Cannot get window extent of text w/o renderer. You likely "
966 "want to call 'figure.draw_without_rendering()' first.")
968 with cbook._setattr_cm(fig, dpi=dpi):
--> 969 bbox, info, descent = self._get_layout(self._renderer)
970 x, y = self.get_unitless_position()
971 x, y = self.get_transform().transform((x, y))
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:373, in Text._get_layout(self, renderer)
370 ys = []
372 # Full vertical extent of font, including ascenders and descenders:
--> 373 _, lp_h, lp_d = _get_text_metrics_with_cache(
374 renderer, "lp", self._fontproperties,
375 ismath="TeX" if self.get_usetex() else False,
376 dpi=self.get_figure(root=True).dpi)
377 min_dy = (lp_h - lp_d) * self._linespacing
379 for i, line in enumerate(lines):
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:69, in _get_text_metrics_with_cache(renderer, text, fontprop, ismath, dpi)
66 """Call ``renderer.get_text_width_height_descent``, caching the results."""
67 # Cached based on a copy of fontprop so that later in-place mutations of
68 # the passed-in argument do not mess up the cache.
---> 69 return _get_text_metrics_with_cache_impl(
70 weakref.ref(renderer), text, fontprop.copy(), ismath, dpi)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:77, in _get_text_metrics_with_cache_impl(renderer_ref, text, fontprop, ismath, dpi)
73 @functools.lru_cache(4096)
74 def _get_text_metrics_with_cache_impl(
75 renderer_ref, text, fontprop, ismath, dpi):
76 # dpi is unused, but participates in cache invalidation (via the renderer).
---> 77 return renderer_ref().get_text_width_height_descent(text, fontprop, ismath)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backends/backend_agg.py:211, in RendererAgg.get_text_width_height_descent(self, s, prop, ismath)
209 _api.check_in_list(["TeX", True, False], ismath=ismath)
210 if ismath == "TeX":
--> 211 return super().get_text_width_height_descent(s, prop, ismath)
213 if ismath:
214 ox, oy, width, height, descent, font_image = \
215 self.mathtext_parser.parse(s, self.dpi, prop)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backend_bases.py:566, in RendererBase.get_text_width_height_descent(self, s, prop, ismath)
562 fontsize = prop.get_size_in_points()
564 if ismath == 'TeX':
565 # todo: handle properties
--> 566 return self.get_texmanager().get_text_width_height_descent(
567 s, fontsize, renderer=self)
569 dpi = self.points_to_pixels(72)
570 if ismath:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:364, in TexManager.get_text_width_height_descent(cls, tex, fontsize, renderer)
362 if tex.strip() == '':
363 return 0, 0, 0
--> 364 dvifile = cls.make_dvi(tex, fontsize)
365 dpi_fraction = renderer.points_to_pixels(1.) if renderer else 1
366 with dviread.Dvi(dvifile, 72 * dpi_fraction) as dvi:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:293, in TexManager.make_dvi(cls, tex, fontsize)
290 with TemporaryDirectory(dir=dvifile.parent) as tmpdir:
291 Path(tmpdir, "file.tex").write_text(
292 cls._get_tex_source(tex, fontsize), encoding='utf-8')
--> 293 cls._run_checked_subprocess(
294 ["latex", "-interaction=nonstopmode", "-halt-on-error",
295 "-no-shell-escape", "file.tex"], tex, cwd=tmpdir)
296 Path(tmpdir, "file.dvi").replace(dvifile)
297 # Also move the tex source to the main cache directory, but
298 # only for backcompat.
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:254, in TexManager._run_checked_subprocess(cls, command, tex, cwd)
250 report = subprocess.check_output(
251 command, cwd=cwd if cwd is not None else cls._texcache,
252 stderr=subprocess.STDOUT)
253 except FileNotFoundError as exc:
--> 254 raise RuntimeError(
255 f'Failed to process string with tex because {command[0]} '
256 'could not be found') from exc
257 except subprocess.CalledProcessError as exc:
258 raise RuntimeError(
259 '{prog} was not able to process the following string:\n'
260 '{tex!r}\n\n'
(...) 267 exc=exc.output.decode('utf-8', 'backslashreplace'))
268 ) from None
RuntimeError: Failed to process string with tex because latex could not be found
<Figure size 640x480 with 1 Axes>
Error of the composite trapezoidal rule#
The error of the composite trapezoidal rule stems from the error of the trapezoidal rule; here we make such an error \(n\) times.
Since \(h\cdot n=b-a\) holds, we derive the error of the composite trapezoidal rule as:
\[E_{\textrm{trapezno sest}}=-\frac{h^2(b-a)}{12}f''(\eta),\]
\(\eta\) is a value on the interval \([a,b]\). The error is of second order \(\mathcal{O}(h^2)\).
A better approximation of the integral#
Next we will look at the so-called Richardson extrapolation, in which we compute a better approximation from the integral estimates with step \(h\) and \(2h\).
If we compute the integral \(I\) numerically at two different steps \(h\) and \(2\,h\), the following holds:
\[\int_a^b f(x)\,\textrm{d}x = I_h + E_h = I_{2h} + E_{2h},\]
where \(I_h\) and \(I_{2h}\) are the approximations of the integral with step \(h\) and \(2h\), and \(E_h\) and \(E_{2h}\) are the error estimates at step \(h\) and \(2h\). We derive \(I_{2h}-I_{h} = E_{h}-E_{2h}\)
Next we write:
\[E_h=-\frac{h^2(b-a)}{12}f''(\eta)=h^2\,K.\]
Under the assumption that \(f''\left (\eta \right )\) is the same at step \(h\) and \(2h\) (\(\eta\) is actually different at step \(h\) and \(2h\)), we write:
\[E_{2h}=-\frac{(2h)^2(b-a)}{12}f''(\eta)=4\,h^2\,K\]
It follows that:
\[I_{2h}-I_h=-3K\,h^2.\]
Now we can estimate the error at step \(h\):
\[E_h=h^2\,K=\frac{I_h-I_{2h}}{3}.\]
Based on the error estimate, we can compute a better numerical approximation \(I_{h}^*\):
\[I_h^* = I_h + \frac{1}{3}\,(I_{h}-I_{2h})\]
or
\[I_h^* = \frac{4}{3}\,I_h - \frac{1}{3}\,I_{2h}\]
Numerical example#
We previously computed the integral with the trapezoidal rule at step \(h=1\) and at step \(h=0.5\) already; the results were:
[I_trapezoidal, I_trapezoidal_comp]
\[\displaystyle \left[ 1.33003291922963, \ 1.41313769956786\right]\]
Using the formula above, we compute a better approximation:
I_trapezoidal_better = 4/3*I_trapezoidal_comp - 1/3*I_trapezoidal
print(f'Exact result: {I_exact}\nBetter approximation: {I_trapezoidal_better}')
Exact result: 1.4404224209802097
Better approximation: 1.4408392930139313
The composite trapezoidal rule is also implemented in the numpy package, via the numpy.trapz function:
trapz(y, x=None, dx=1.0, axis=-1)
yrepresents the table of function values,xis an optional parameter and defines the nodes; if the parameter is not defined, equidistant nodes spaced bydxare assumed,dxdefines the (constant) integration step, with a default value of 1,axisdefines the axis along which to integrate (in caseyis a multidimensional numeric array).
The function returns the computed integral according to the composite trapezoidal rule. More information can be found in the documentation.
numpy implementation of the composite trapezoidal rule#
Let us look at an example:
#%%timeit
I_trapezoidal_np = np.trapezoid(y3v, dx=h3v)
I_trapezoidal_np
\[\displaystyle 1.41313769956786\]
Simpson’s and other methods#
Above we looked at the trapezoidal rule, which is based on a linear interpolating function on an individual subinterval. With higher-order interpolation we can derive further integration methods.
We need to compute:
\[\int_{a}^b f(x)\,dx.\]
We tabulate the integrand \(f(x)\) and interpolate the table with the Lagrange interpolating polynomial \(P_n(x)\) of degree \(n\):
\[P_n(x)=\sum_{i=0}^{n}\,f(x_i)\,l_i(x),\]
where \(y_i=f(x_i)\) are the function values at the nodes \(x_i\) and the Lagrange polynomial \(l_i\) is defined as:
\[l_i(x)=\prod_{j=0, j\ne i}^n \frac{x-x_j}{x_i-x_j}.\]
For the numerical computation of the integral \(\int_{a}^b f(x)\,dx\) (the limits are: \(a=x_0\), \(b=x_n\)), instead of the function \(f(x)\) we substitute the interpolating polynomial \(P_n(x)\) into the integral:
\[I=\int_{x_0}^{x_{n}} P_n(x)\,dx=\int_{x_0}^{x_{n}} \sum_{i=0}^{n}\,f(x_i)\,l_i(x)\,dx.\]
Since integration is a linear operation, we can swap the integration and the sum:
\[I=\sum_{i=0}^{n}\,f(x_i)\,\underbrace{\int_{x_0}^{x_{n}} l_i(x)\,dx}_{A_i}.\]
We integrate the Lagrange polynomial and obtain the weights \(A_i\):
\[A_i = \int_{x_0}^{x_{n}} l_i(x)\,dx\]
Derivation of the trapezoidal rule using Lagrange polynomials#
Let us look at how we can computationally derive the weights \(A_i\) for the case of the trapezoidal rule using a first-degree Lagrange interpolating polynomial.
First, let us define the variable x, the nodes x0 and x1, and the step h in symbolic form:
x, x0, x1, h = sym.symbols('x x0, x1, h')
Let us prepare a Python function that returns, in symbolic form, a list of \(n\) coefficients of the Lagrange polynomials \([l_0(x), l_1(x),\dots, l_{n-1}(x)]\) of degree \(n-1\):
def lagrange(n, x, nodes_prefix='x'):
if isinstance(nodes_prefix, str):
nodes = sym.symbols(f'{nodes_prefix}:{n}')
coeffs = []
for i in range(0, n):
numer = []
denom = []
for j in range(0, n):
if i == j:
continue
numer.append(x - nodes[j])
denom.append(nodes[i] - nodes[j])
numer = sym.Mul(*numer)
denom = sym.Mul(*denom)
coeffs.append(numer/denom)
return coeffs
First, let us look at the Lagrange polynomials for linear interpolation (\(n=2\)):
lag = lagrange(2, x)
lag
\[\displaystyle \left[ \frac{x - x_{1}}{x_{0} - x_{1}}, \ \frac{x - x_{0}}{- x_{0} + x_{1}}\right]\]
Now we integrate the Lagrange polynomial \(l_0(x)\) over the entire interval:
int0 = sym.integrate(lag[0], (x, x0, x1))
int0
\[\displaystyle - \frac{x_{0}^{2}}{2 x_{0} - 2 x_{1}} + \frac{x_{0} x_{1}}{x_{0} - x_{1}} + \frac{x_{1}^{2}}{2 x_{0} - 2 x_{1}} - \frac{x_{1}^{2}}{x_{0} - x_{1}}\]
We simplify the expression and obtain:
int1 = int0.factor()
int1
\[\displaystyle - \frac{x_{0} - x_{1}}{2}\]
Since the width of the subinterval is constant, \(x_1=h_0+h\), we perform the substitution:
substitutions = {x1: x0+h}
int1.subs(substitutions)
\[\displaystyle \frac{h}{2}\]
The steps above for the Lagrange polynomial \(l_0(x)\) can be generalized to a list of Lagrange polynomials:
x, x0, x1, h = sym.symbols('x, x0, x1, h')
substitutions = {x1: x0+h}
A_trapezoid = [sym.integrate(li, (x, x0, x1)).factor().subs(substitutions)
for li in lagrange(2, x)] # for each Lagrange polynomial `li` in the list lagrange(2,x)
A_trapezoid
\[\displaystyle \left[ \frac{h}{2}, \ \frac{h}{2}\right]\]
We derived the weights that we used in the trapezoidal method:
\[A_0=h/2\qquad A_1=h/2.\]
The trapezoidal rule is:
\[I_{\textrm{trapezno}}=\frac{h}{2}\left(y_0+y_1\right)\]
The error estimate is (reference: Burden, Faires, Burden: Numerical Analysis 10th Ed):
\[E_{\textrm{trapezno}}= -\frac{h^3}{12}f''(\xi).\]
Computing the weights for Simpson’s rule#
After showing the concise computation for the trapezoidal rule above, we can proceed similarly for quadratic interpolation over three points (\(n=3\)).
The computation of the weights is analogous to the one above:
x, x0, x1, x2, h = sym.symbols('x, x0, x1, x2, h')
substitutions = {x1: x0+h, x2: x0+2*h}
A_Simpson1_3 = [sym.integrate(li, (x, x0, x2)).factor().subs(substitutions).factor()
for li in lagrange(3, x)]
A_Simpson1_3
\[\displaystyle \left[ \frac{h}{3}, \ \frac{4 h}{3}, \ \frac{h}{3}\right]\]
Simpson’s rule (this rule is also called Simpson’s 1/3 rule) is:
\[I_{\textrm{Simpsonovo}}=\frac{h}{3}\left(y_0+4\,y_1+y_2\right)\]
The error estimate is (reference: Burden, Faires, Burden: Numerical Analysis 10th Ed):
\[E_{\textrm{Simpsonovo}}= -\frac{h^5}{90}f^{(4)}(\xi).\]
It should be pointed out that the error is locally of fifth order \(\mathcal{O}(h^5)\) and is defined by the unknown value of the fourth derivative \(f^{(4)}(\xi)\); consequently, this rule is exact for polynomials of degree 3 or less.
Example of use:
I_Simps = h3v/3 * np.sum(y3v * [1, 4, 1])
I_Simps
\[\displaystyle 1.44083929301393\]
Let us prepare a figure. Since Simpson’s rule is based on quadratic interpolation, we must first prepare the interpolating polynomial (we use scipy.interpolate):
from scipy import interpolate
def fig_simpson():
y_interpolate = interpolate.lagrange(x3v, y3v)
plt.fill_between(xg, y_interpolate(xg), alpha=0.25, facecolor='r')
plt.vlines(x3v, 0, y3v, color='r', linestyles='dashed', lw=1)
plt.annotate('$I_{\\textrm{Simpsonovo}}$', (1.2, 0.5), fontsize=22)
plt.annotate('Error', fontsize=20, xy=(1.75, 1.7), xytext=(1.4, 1.8),
arrowprops=dict(facecolor='gray', shrink=0.05))
plt.annotate('Error', fontsize=20, xy=(1.2, 1.1), xytext=(1., 1.4),
arrowprops=dict(facecolor='gray', shrink=0.05))
plt.plot(xg, yg, lw=3, alpha=0.5, label='$x\\,\\sin(x)$')
plt.plot(x3v, y3v, 'o', alpha=0.5, label=f'$h={h3v}$')
plt.legend(loc=(1.01, 0))
plt.ylim(0, 2)
plt.show()
We display:
fig_simpson()
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:250, in TexManager._run_checked_subprocess(cls, command, tex, cwd)
249 try:
--> 250 report = subprocess.check_output(
251 command, cwd=cwd if cwd is not None else cls._texcache,
252 stderr=subprocess.STDOUT)
253 except FileNotFoundError as exc:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:466, in check_output(timeout, *popenargs, **kwargs)
464 kwargs['input'] = empty
--> 466 return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
467 **kwargs).stdout
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:548, in run(input, capture_output, timeout, check, *popenargs, **kwargs)
546 kwargs['stderr'] = PIPE
--> 548 with Popen(*popenargs, **kwargs) as process:
549 try:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:1026, in Popen.__init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, user, group, extra_groups, encoding, errors, text, umask, pipesize, process_group)
1023 self.stderr = io.TextIOWrapper(self.stderr,
1024 encoding=encoding, errors=errors)
-> 1026 self._execute_child(args, executable, preexec_fn, close_fds,
1027 pass_fds, cwd, env,
1028 startupinfo, creationflags, shell,
1029 p2cread, p2cwrite,
1030 c2pread, c2pwrite,
1031 errread, errwrite,
1032 restore_signals,
1033 gid, gids, uid, umask,
1034 start_new_session, process_group)
1035 except:
1036 # Cleanup if the child failed starting.
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:1955, in Popen._execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, gid, gids, uid, umask, start_new_session, process_group)
1954 if err_filename is not None:
-> 1955 raise child_exception_type(errno_num, err_msg, err_filename)
1956 else:
FileNotFoundError: [Errno 2] No such file or directory: 'latex'
The above exception was the direct cause of the following exception:
RuntimeError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backend_bases.py:2157, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
2154 # we do this instead of `self.figure.draw_without_rendering`
2155 # so that we can inject the orientation
2156 with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2157 self.figure.draw(renderer)
2158 if bbox_inches:
2159 if bbox_inches == "tight":
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:94, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
92 @wraps(draw)
93 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 94 result = draw(artist, renderer, *args, **kwargs)
95 if renderer._rasterizing:
96 renderer.stop_rasterizing()
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/figure.py:3264, in Figure.draw(self, renderer)
3261 # ValueError can occur when resizing a window.
3263 self.patch.draw(renderer)
-> 3264 mimage._draw_list_compositing_images(
3265 renderer, self, artists, self.suppressComposite)
3267 renderer.close_group('figure')
3268 finally:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axes/_base.py:3226, in _AxesBase.draw(self, renderer)
3223 if artists_rasterized:
3224 _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3226 mimage._draw_list_compositing_images(
3227 renderer, self, artists, self.get_figure(root=True).suppressComposite)
3229 renderer.close_group('axes')
3230 self.stale = False
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1405, in Axis.draw(self, renderer)
1402 renderer.open_group(__name__, gid=self.get_gid())
1404 ticks_to_draw = self._update_ticks()
-> 1405 tlb1, tlb2 = self._get_ticklabel_bboxes(ticks_to_draw, renderer)
1407 for tick in ticks_to_draw:
1408 tick.draw(renderer)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1332, in Axis._get_ticklabel_bboxes(self, ticks, renderer)
1330 if renderer is None:
1331 renderer = self.get_figure(root=True)._get_renderer()
-> 1332 return ([tick.label1.get_window_extent(renderer)
1333 for tick in ticks if tick.label1.get_visible()],
1334 [tick.label2.get_window_extent(renderer)
1335 for tick in ticks if tick.label2.get_visible()])
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1332, in <listcomp>(.0)
1330 if renderer is None:
1331 renderer = self.get_figure(root=True)._get_renderer()
-> 1332 return ([tick.label1.get_window_extent(renderer)
1333 for tick in ticks if tick.label1.get_visible()],
1334 [tick.label2.get_window_extent(renderer)
1335 for tick in ticks if tick.label2.get_visible()])
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:969, in Text.get_window_extent(self, renderer, dpi)
964 raise RuntimeError(
965 "Cannot get window extent of text w/o renderer. You likely "
966 "want to call 'figure.draw_without_rendering()' first.")
968 with cbook._setattr_cm(fig, dpi=dpi):
--> 969 bbox, info, descent = self._get_layout(self._renderer)
970 x, y = self.get_unitless_position()
971 x, y = self.get_transform().transform((x, y))
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:373, in Text._get_layout(self, renderer)
370 ys = []
372 # Full vertical extent of font, including ascenders and descenders:
--> 373 _, lp_h, lp_d = _get_text_metrics_with_cache(
374 renderer, "lp", self._fontproperties,
375 ismath="TeX" if self.get_usetex() else False,
376 dpi=self.get_figure(root=True).dpi)
377 min_dy = (lp_h - lp_d) * self._linespacing
379 for i, line in enumerate(lines):
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:69, in _get_text_metrics_with_cache(renderer, text, fontprop, ismath, dpi)
66 """Call ``renderer.get_text_width_height_descent``, caching the results."""
67 # Cached based on a copy of fontprop so that later in-place mutations of
68 # the passed-in argument do not mess up the cache.
---> 69 return _get_text_metrics_with_cache_impl(
70 weakref.ref(renderer), text, fontprop.copy(), ismath, dpi)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:77, in _get_text_metrics_with_cache_impl(renderer_ref, text, fontprop, ismath, dpi)
73 @functools.lru_cache(4096)
74 def _get_text_metrics_with_cache_impl(
75 renderer_ref, text, fontprop, ismath, dpi):
76 # dpi is unused, but participates in cache invalidation (via the renderer).
---> 77 return renderer_ref().get_text_width_height_descent(text, fontprop, ismath)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backends/backend_agg.py:211, in RendererAgg.get_text_width_height_descent(self, s, prop, ismath)
209 _api.check_in_list(["TeX", True, False], ismath=ismath)
210 if ismath == "TeX":
--> 211 return super().get_text_width_height_descent(s, prop, ismath)
213 if ismath:
214 ox, oy, width, height, descent, font_image = \
215 self.mathtext_parser.parse(s, self.dpi, prop)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backend_bases.py:566, in RendererBase.get_text_width_height_descent(self, s, prop, ismath)
562 fontsize = prop.get_size_in_points()
564 if ismath == 'TeX':
565 # todo: handle properties
--> 566 return self.get_texmanager().get_text_width_height_descent(
567 s, fontsize, renderer=self)
569 dpi = self.points_to_pixels(72)
570 if ismath:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:364, in TexManager.get_text_width_height_descent(cls, tex, fontsize, renderer)
362 if tex.strip() == '':
363 return 0, 0, 0
--> 364 dvifile = cls.make_dvi(tex, fontsize)
365 dpi_fraction = renderer.points_to_pixels(1.) if renderer else 1
366 with dviread.Dvi(dvifile, 72 * dpi_fraction) as dvi:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:293, in TexManager.make_dvi(cls, tex, fontsize)
290 with TemporaryDirectory(dir=dvifile.parent) as tmpdir:
291 Path(tmpdir, "file.tex").write_text(
292 cls._get_tex_source(tex, fontsize), encoding='utf-8')
--> 293 cls._run_checked_subprocess(
294 ["latex", "-interaction=nonstopmode", "-halt-on-error",
295 "-no-shell-escape", "file.tex"], tex, cwd=tmpdir)
296 Path(tmpdir, "file.dvi").replace(dvifile)
297 # Also move the tex source to the main cache directory, but
298 # only for backcompat.
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:254, in TexManager._run_checked_subprocess(cls, command, tex, cwd)
250 report = subprocess.check_output(
251 command, cwd=cwd if cwd is not None else cls._texcache,
252 stderr=subprocess.STDOUT)
253 except FileNotFoundError as exc:
--> 254 raise RuntimeError(
255 f'Failed to process string with tex because {command[0]} '
256 'could not be found') from exc
257 except subprocess.CalledProcessError as exc:
258 raise RuntimeError(
259 '{prog} was not able to process the following string:\n'
260 '{tex!r}\n\n'
(...) 267 exc=exc.output.decode('utf-8', 'backslashreplace'))
268 ) from None
RuntimeError: Failed to process string with tex because latex could not be found
<Figure size 640x480 with 1 Axes>
scipy.integrate.newton_cotes#
The coefficients of the Newton-Cotes integration approach can also be obtained using scipy.integrate.newton_cotes() (documentation):
newton_cotes(rn, equal=0)
where the parameters are:
rn, which defines the number of subintervals (a non-constant step is also possible, see the documentation),equal, which defines whether equally wide subintervals are required.
The function returns a tuple, in which the first element represents the numeric array of weights and the second element the error estimate.
Let us look at an example:
from scipy import integrate
integrate.newton_cotes(2)
(array([0.33333333, 1.33333333, 0.33333333]), -0.011111111111111112)
Computing the weights for Simpson’s 3/8 rule#
We can continue with cubic interpolation over four points (\(n=4\)):
x, x0, x1, x2, x3, h = sym.symbols('x, x0, x1, x2, x3, h')
substitutions = {x1: x0+h, x2: x0+2*h, x3: x0+3*h}
A_Simpson3_8 = [sym.integrate(li, (x, x0, x3)).factor().subs(substitutions).factor()
for li in lagrange(4, x)]
A_Simpson3_8
\[\displaystyle \left[ \frac{3 h}{8}, \ \frac{9 h}{8}, \ \frac{9 h}{8}, \ \frac{3 h}{8}\right]\]
Simpson’s 3/8 rule is:
\[I_{\textrm{Simpsonovo 3/8}}=\frac{3h}{8}\left(y_0+3\,y_1+3\,y_2+y_3\right)\]
The error estimate is (reference: Burden, Faires, Burden: Numerical Analysis 10th Ed):
\[E_{\textrm{Simpsonovo 3/8}}= -\frac{3\,h^5}{80}f^{(4)}(\xi).\]
Let us look at an example of use. We use the previously prepared table of function values at four points:
y4v
array([0.84147098, 1.2959172 , 1.65901326, 1.81859485])
I_Simps38 = 3*h4v/8 * np.sum(y4v * [1, 3, 3, 1])
I_Simps38
\[\displaystyle 1.44060715408392\]
Let us also prepare a plot:
def fig_simpson38():
y_interpolate = interpolate.lagrange(x4v, y4v)
plt.fill_between(xg, y_interpolate(xg), alpha=0.25, facecolor='r')
plt.vlines(x4v, 0, y4v, color='r', linestyles='dashed', lw=1)
plt.annotate('$I_{\\textrm{Simpsonovo 3/8}}$', (1.2, 0.5), fontsize=22)
plt.annotate('Error', fontsize=20, xy=(1.75, 1.7), xytext=(1.4, 1.8),
arrowprops=dict(facecolor='gray', shrink=0.05))
plt.annotate('Error', fontsize=20, xy=(1.5, 1.47), xytext=(1.1, 1.6),
arrowprops=dict(facecolor='gray', shrink=0.05))
plt.annotate('Error', fontsize=20, xy=(1.2, 1.1), xytext=(1., 1.4),
arrowprops=dict(facecolor='gray', shrink=0.05))
plt.plot(xg, yg, lw=3, alpha=0.5, label='$x\\,\\sin(x)$')
plt.plot(x4v, y4v, 'o', alpha=0.5, label=f'$h={h4v:.6f}$')
plt.legend(loc=(1.01, 0))
plt.ylim(0, 2)
plt.show()
We display:
fig_simpson38()
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:250, in TexManager._run_checked_subprocess(cls, command, tex, cwd)
249 try:
--> 250 report = subprocess.check_output(
251 command, cwd=cwd if cwd is not None else cls._texcache,
252 stderr=subprocess.STDOUT)
253 except FileNotFoundError as exc:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:466, in check_output(timeout, *popenargs, **kwargs)
464 kwargs['input'] = empty
--> 466 return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
467 **kwargs).stdout
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:548, in run(input, capture_output, timeout, check, *popenargs, **kwargs)
546 kwargs['stderr'] = PIPE
--> 548 with Popen(*popenargs, **kwargs) as process:
549 try:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:1026, in Popen.__init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, user, group, extra_groups, encoding, errors, text, umask, pipesize, process_group)
1023 self.stderr = io.TextIOWrapper(self.stderr,
1024 encoding=encoding, errors=errors)
-> 1026 self._execute_child(args, executable, preexec_fn, close_fds,
1027 pass_fds, cwd, env,
1028 startupinfo, creationflags, shell,
1029 p2cread, p2cwrite,
1030 c2pread, c2pwrite,
1031 errread, errwrite,
1032 restore_signals,
1033 gid, gids, uid, umask,
1034 start_new_session, process_group)
1035 except:
1036 # Cleanup if the child failed starting.
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:1955, in Popen._execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, gid, gids, uid, umask, start_new_session, process_group)
1954 if err_filename is not None:
-> 1955 raise child_exception_type(errno_num, err_msg, err_filename)
1956 else:
FileNotFoundError: [Errno 2] No such file or directory: 'latex'
The above exception was the direct cause of the following exception:
RuntimeError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backend_bases.py:2157, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
2154 # we do this instead of `self.figure.draw_without_rendering`
2155 # so that we can inject the orientation
2156 with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2157 self.figure.draw(renderer)
2158 if bbox_inches:
2159 if bbox_inches == "tight":
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:94, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
92 @wraps(draw)
93 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 94 result = draw(artist, renderer, *args, **kwargs)
95 if renderer._rasterizing:
96 renderer.stop_rasterizing()
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/figure.py:3264, in Figure.draw(self, renderer)
3261 # ValueError can occur when resizing a window.
3263 self.patch.draw(renderer)
-> 3264 mimage._draw_list_compositing_images(
3265 renderer, self, artists, self.suppressComposite)
3267 renderer.close_group('figure')
3268 finally:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axes/_base.py:3226, in _AxesBase.draw(self, renderer)
3223 if artists_rasterized:
3224 _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3226 mimage._draw_list_compositing_images(
3227 renderer, self, artists, self.get_figure(root=True).suppressComposite)
3229 renderer.close_group('axes')
3230 self.stale = False
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1405, in Axis.draw(self, renderer)
1402 renderer.open_group(__name__, gid=self.get_gid())
1404 ticks_to_draw = self._update_ticks()
-> 1405 tlb1, tlb2 = self._get_ticklabel_bboxes(ticks_to_draw, renderer)
1407 for tick in ticks_to_draw:
1408 tick.draw(renderer)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1332, in Axis._get_ticklabel_bboxes(self, ticks, renderer)
1330 if renderer is None:
1331 renderer = self.get_figure(root=True)._get_renderer()
-> 1332 return ([tick.label1.get_window_extent(renderer)
1333 for tick in ticks if tick.label1.get_visible()],
1334 [tick.label2.get_window_extent(renderer)
1335 for tick in ticks if tick.label2.get_visible()])
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1332, in <listcomp>(.0)
1330 if renderer is None:
1331 renderer = self.get_figure(root=True)._get_renderer()
-> 1332 return ([tick.label1.get_window_extent(renderer)
1333 for tick in ticks if tick.label1.get_visible()],
1334 [tick.label2.get_window_extent(renderer)
1335 for tick in ticks if tick.label2.get_visible()])
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:969, in Text.get_window_extent(self, renderer, dpi)
964 raise RuntimeError(
965 "Cannot get window extent of text w/o renderer. You likely "
966 "want to call 'figure.draw_without_rendering()' first.")
968 with cbook._setattr_cm(fig, dpi=dpi):
--> 969 bbox, info, descent = self._get_layout(self._renderer)
970 x, y = self.get_unitless_position()
971 x, y = self.get_transform().transform((x, y))
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:373, in Text._get_layout(self, renderer)
370 ys = []
372 # Full vertical extent of font, including ascenders and descenders:
--> 373 _, lp_h, lp_d = _get_text_metrics_with_cache(
374 renderer, "lp", self._fontproperties,
375 ismath="TeX" if self.get_usetex() else False,
376 dpi=self.get_figure(root=True).dpi)
377 min_dy = (lp_h - lp_d) * self._linespacing
379 for i, line in enumerate(lines):
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:69, in _get_text_metrics_with_cache(renderer, text, fontprop, ismath, dpi)
66 """Call ``renderer.get_text_width_height_descent``, caching the results."""
67 # Cached based on a copy of fontprop so that later in-place mutations of
68 # the passed-in argument do not mess up the cache.
---> 69 return _get_text_metrics_with_cache_impl(
70 weakref.ref(renderer), text, fontprop.copy(), ismath, dpi)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:77, in _get_text_metrics_with_cache_impl(renderer_ref, text, fontprop, ismath, dpi)
73 @functools.lru_cache(4096)
74 def _get_text_metrics_with_cache_impl(
75 renderer_ref, text, fontprop, ismath, dpi):
76 # dpi is unused, but participates in cache invalidation (via the renderer).
---> 77 return renderer_ref().get_text_width_height_descent(text, fontprop, ismath)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backends/backend_agg.py:211, in RendererAgg.get_text_width_height_descent(self, s, prop, ismath)
209 _api.check_in_list(["TeX", True, False], ismath=ismath)
210 if ismath == "TeX":
--> 211 return super().get_text_width_height_descent(s, prop, ismath)
213 if ismath:
214 ox, oy, width, height, descent, font_image = \
215 self.mathtext_parser.parse(s, self.dpi, prop)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backend_bases.py:566, in RendererBase.get_text_width_height_descent(self, s, prop, ismath)
562 fontsize = prop.get_size_in_points()
564 if ismath == 'TeX':
565 # todo: handle properties
--> 566 return self.get_texmanager().get_text_width_height_descent(
567 s, fontsize, renderer=self)
569 dpi = self.points_to_pixels(72)
570 if ismath:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:364, in TexManager.get_text_width_height_descent(cls, tex, fontsize, renderer)
362 if tex.strip() == '':
363 return 0, 0, 0
--> 364 dvifile = cls.make_dvi(tex, fontsize)
365 dpi_fraction = renderer.points_to_pixels(1.) if renderer else 1
366 with dviread.Dvi(dvifile, 72 * dpi_fraction) as dvi:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:293, in TexManager.make_dvi(cls, tex, fontsize)
290 with TemporaryDirectory(dir=dvifile.parent) as tmpdir:
291 Path(tmpdir, "file.tex").write_text(
292 cls._get_tex_source(tex, fontsize), encoding='utf-8')
--> 293 cls._run_checked_subprocess(
294 ["latex", "-interaction=nonstopmode", "-halt-on-error",
295 "-no-shell-escape", "file.tex"], tex, cwd=tmpdir)
296 Path(tmpdir, "file.dvi").replace(dvifile)
297 # Also move the tex source to the main cache directory, but
298 # only for backcompat.
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:254, in TexManager._run_checked_subprocess(cls, command, tex, cwd)
250 report = subprocess.check_output(
251 command, cwd=cwd if cwd is not None else cls._texcache,
252 stderr=subprocess.STDOUT)
253 except FileNotFoundError as exc:
--> 254 raise RuntimeError(
255 f'Failed to process string with tex because {command[0]} '
256 'could not be found') from exc
257 except subprocess.CalledProcessError as exc:
258 raise RuntimeError(
259 '{prog} was not able to process the following string:\n'
260 '{tex!r}\n\n'
(...) 267 exc=exc.output.decode('utf-8', 'backslashreplace'))
268 ) from None
RuntimeError: Failed to process string with tex because latex could not be found
<Figure size 640x480 with 1 Axes>
Composite Simpson’s rule#
We divide the interval \([a, b]\) into an even number \(n\) of subintervals of equal width \(h=(b - a)/n\), which defines the nodes \(x_i=a+i\,h\) for \(i=0,1,\dots,n\). In this case we write the composite Simpson’s rule:
\[\int_{a}^{b}f(x)\,\textrm{d}x=
\frac{h}{3}\left(
f(x_0)
+4\sum_{i=1}^{n/2}f(x_{2i-1})
+2\sum_{i=1}^{n/2-1}f(x_{2i})
+f(x_n)
\right)
\underbrace{
-\frac{b-a}{180}h^4\,f^{(4)}(\eta)
}_{E_{\textrm{Sestavljeno Simpsonovo 1/3}}}
,\]
where \(\eta\) is an unknown value on the interval \([a, b]\). The error is of fourth order \(\mathcal{O}(h^4)\).
Numerical implementation:
def simpson_comp(y, h):
"""
Composite Simpson's rule of integration.
:param y: function values
:param h: integration step
"""
return h/3 * (y[0] + 4*np.sum(y[1:-1:2]) + 2*np.sum(y[2:-1:2]) + y[-1])
I_Simps_comp = simpson_comp(y5v, h=h5v)
I_Simps_comp
\[\displaystyle 1.44044796026391\]
Let us prepare a figure:
def fig_simpson_comp():
y_interpolate = interpolate.lagrange(x5v, y5v)
plt.fill_between(xg, y_interpolate(xg), alpha=0.25, facecolor='r')
plt.vlines(x5v, 0, y5v, color='r', linestyles='dashed', lw=1)
plt.annotate('$I_{\\textrm{Simpsonovo sestavljeno}}$', (1.2, 0.5), fontsize=22)
plt.annotate('Error', fontsize=20, xy=(1.70, 1.68), xytext=(1.4, 1.8),
arrowprops=dict(facecolor='gray', shrink=0.05))
plt.plot(xg, yg, lw=3, alpha=0.5, label='$x\\,\\sin(x)$')
plt.plot(x5v, y5v, 'o', alpha=0.5, label=f'$h={h5v}$')
plt.legend(loc=(1.01, 0))
plt.ylim(0, 2)
plt.show()
We display:
fig_simpson_comp()
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:250, in TexManager._run_checked_subprocess(cls, command, tex, cwd)
249 try:
--> 250 report = subprocess.check_output(
251 command, cwd=cwd if cwd is not None else cls._texcache,
252 stderr=subprocess.STDOUT)
253 except FileNotFoundError as exc:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:466, in check_output(timeout, *popenargs, **kwargs)
464 kwargs['input'] = empty
--> 466 return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
467 **kwargs).stdout
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:548, in run(input, capture_output, timeout, check, *popenargs, **kwargs)
546 kwargs['stderr'] = PIPE
--> 548 with Popen(*popenargs, **kwargs) as process:
549 try:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:1026, in Popen.__init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, user, group, extra_groups, encoding, errors, text, umask, pipesize, process_group)
1023 self.stderr = io.TextIOWrapper(self.stderr,
1024 encoding=encoding, errors=errors)
-> 1026 self._execute_child(args, executable, preexec_fn, close_fds,
1027 pass_fds, cwd, env,
1028 startupinfo, creationflags, shell,
1029 p2cread, p2cwrite,
1030 c2pread, c2pwrite,
1031 errread, errwrite,
1032 restore_signals,
1033 gid, gids, uid, umask,
1034 start_new_session, process_group)
1035 except:
1036 # Cleanup if the child failed starting.
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:1955, in Popen._execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, gid, gids, uid, umask, start_new_session, process_group)
1954 if err_filename is not None:
-> 1955 raise child_exception_type(errno_num, err_msg, err_filename)
1956 else:
FileNotFoundError: [Errno 2] No such file or directory: 'latex'
The above exception was the direct cause of the following exception:
RuntimeError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backend_bases.py:2157, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
2154 # we do this instead of `self.figure.draw_without_rendering`
2155 # so that we can inject the orientation
2156 with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2157 self.figure.draw(renderer)
2158 if bbox_inches:
2159 if bbox_inches == "tight":
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:94, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
92 @wraps(draw)
93 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 94 result = draw(artist, renderer, *args, **kwargs)
95 if renderer._rasterizing:
96 renderer.stop_rasterizing()
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/figure.py:3264, in Figure.draw(self, renderer)
3261 # ValueError can occur when resizing a window.
3263 self.patch.draw(renderer)
-> 3264 mimage._draw_list_compositing_images(
3265 renderer, self, artists, self.suppressComposite)
3267 renderer.close_group('figure')
3268 finally:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axes/_base.py:3226, in _AxesBase.draw(self, renderer)
3223 if artists_rasterized:
3224 _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3226 mimage._draw_list_compositing_images(
3227 renderer, self, artists, self.get_figure(root=True).suppressComposite)
3229 renderer.close_group('axes')
3230 self.stale = False
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1405, in Axis.draw(self, renderer)
1402 renderer.open_group(__name__, gid=self.get_gid())
1404 ticks_to_draw = self._update_ticks()
-> 1405 tlb1, tlb2 = self._get_ticklabel_bboxes(ticks_to_draw, renderer)
1407 for tick in ticks_to_draw:
1408 tick.draw(renderer)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1332, in Axis._get_ticklabel_bboxes(self, ticks, renderer)
1330 if renderer is None:
1331 renderer = self.get_figure(root=True)._get_renderer()
-> 1332 return ([tick.label1.get_window_extent(renderer)
1333 for tick in ticks if tick.label1.get_visible()],
1334 [tick.label2.get_window_extent(renderer)
1335 for tick in ticks if tick.label2.get_visible()])
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1332, in <listcomp>(.0)
1330 if renderer is None:
1331 renderer = self.get_figure(root=True)._get_renderer()
-> 1332 return ([tick.label1.get_window_extent(renderer)
1333 for tick in ticks if tick.label1.get_visible()],
1334 [tick.label2.get_window_extent(renderer)
1335 for tick in ticks if tick.label2.get_visible()])
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:969, in Text.get_window_extent(self, renderer, dpi)
964 raise RuntimeError(
965 "Cannot get window extent of text w/o renderer. You likely "
966 "want to call 'figure.draw_without_rendering()' first.")
968 with cbook._setattr_cm(fig, dpi=dpi):
--> 969 bbox, info, descent = self._get_layout(self._renderer)
970 x, y = self.get_unitless_position()
971 x, y = self.get_transform().transform((x, y))
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:373, in Text._get_layout(self, renderer)
370 ys = []
372 # Full vertical extent of font, including ascenders and descenders:
--> 373 _, lp_h, lp_d = _get_text_metrics_with_cache(
374 renderer, "lp", self._fontproperties,
375 ismath="TeX" if self.get_usetex() else False,
376 dpi=self.get_figure(root=True).dpi)
377 min_dy = (lp_h - lp_d) * self._linespacing
379 for i, line in enumerate(lines):
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:69, in _get_text_metrics_with_cache(renderer, text, fontprop, ismath, dpi)
66 """Call ``renderer.get_text_width_height_descent``, caching the results."""
67 # Cached based on a copy of fontprop so that later in-place mutations of
68 # the passed-in argument do not mess up the cache.
---> 69 return _get_text_metrics_with_cache_impl(
70 weakref.ref(renderer), text, fontprop.copy(), ismath, dpi)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:77, in _get_text_metrics_with_cache_impl(renderer_ref, text, fontprop, ismath, dpi)
73 @functools.lru_cache(4096)
74 def _get_text_metrics_with_cache_impl(
75 renderer_ref, text, fontprop, ismath, dpi):
76 # dpi is unused, but participates in cache invalidation (via the renderer).
---> 77 return renderer_ref().get_text_width_height_descent(text, fontprop, ismath)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backends/backend_agg.py:211, in RendererAgg.get_text_width_height_descent(self, s, prop, ismath)
209 _api.check_in_list(["TeX", True, False], ismath=ismath)
210 if ismath == "TeX":
--> 211 return super().get_text_width_height_descent(s, prop, ismath)
213 if ismath:
214 ox, oy, width, height, descent, font_image = \
215 self.mathtext_parser.parse(s, self.dpi, prop)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backend_bases.py:566, in RendererBase.get_text_width_height_descent(self, s, prop, ismath)
562 fontsize = prop.get_size_in_points()
564 if ismath == 'TeX':
565 # todo: handle properties
--> 566 return self.get_texmanager().get_text_width_height_descent(
567 s, fontsize, renderer=self)
569 dpi = self.points_to_pixels(72)
570 if ismath:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:364, in TexManager.get_text_width_height_descent(cls, tex, fontsize, renderer)
362 if tex.strip() == '':
363 return 0, 0, 0
--> 364 dvifile = cls.make_dvi(tex, fontsize)
365 dpi_fraction = renderer.points_to_pixels(1.) if renderer else 1
366 with dviread.Dvi(dvifile, 72 * dpi_fraction) as dvi:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:293, in TexManager.make_dvi(cls, tex, fontsize)
290 with TemporaryDirectory(dir=dvifile.parent) as tmpdir:
291 Path(tmpdir, "file.tex").write_text(
292 cls._get_tex_source(tex, fontsize), encoding='utf-8')
--> 293 cls._run_checked_subprocess(
294 ["latex", "-interaction=nonstopmode", "-halt-on-error",
295 "-no-shell-escape", "file.tex"], tex, cwd=tmpdir)
296 Path(tmpdir, "file.dvi").replace(dvifile)
297 # Also move the tex source to the main cache directory, but
298 # only for backcompat.
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:254, in TexManager._run_checked_subprocess(cls, command, tex, cwd)
250 report = subprocess.check_output(
251 command, cwd=cwd if cwd is not None else cls._texcache,
252 stderr=subprocess.STDOUT)
253 except FileNotFoundError as exc:
--> 254 raise RuntimeError(
255 f'Failed to process string with tex because {command[0]} '
256 'could not be found') from exc
257 except subprocess.CalledProcessError as exc:
258 raise RuntimeError(
259 '{prog} was not able to process the following string:\n'
260 '{tex!r}\n\n'
(...) 267 exc=exc.output.decode('utf-8', 'backslashreplace'))
268 ) from None
RuntimeError: Failed to process string with tex because latex could not be found
<Figure size 640x480 with 1 Axes>
A better estimate of the integral#
The error of the composite Simpson’s method is defined by:
\[E_{\textrm{Sestavljeno Simpsonovo 1/3}}=
-\frac{b-a}{180}h^4\,f^{(4)}(\eta),\]
where \(\eta\) is an unknown value on the interval \([a, b]\).
We determine an improved estimate of the integral in a manner similar to the composite trapezoidal method; we estimate the integral \(I\) at two different steps \(h\) and \(2\,h\), and the following holds exactly:
\[\int_a^b f(x)\,\textrm{d}x = I_h + E_h = I_{2h} + E_{2h},\]
where \(I_h\) is the approximation of the integral with step \(h\) and \(E_h\) is the error estimate at step \(h\); the same holds analogously for \(I_{2h}\) and \(E_{2h}\).
If we assume that \(f^{(4)}\left (\eta \right )\) is the same in both cases, we can determine the difference \(I_{2h}-I_{h} = E_{h}-E_{2h}\).
Next we write:
\[E_h=-\frac{b-a}{180}h^4\,f^{(4)}(\eta)=h^4\,K.\]
Under the assumption that \(f^{(4)}\left (\eta \right )\) is the same at step \(h\) and \(2h\) (\(\eta\) is actually different at step \(h\) and \(2h\)), we write:
\[E_{2h}=-\frac{(b-a)}{180}(2h)^4\,f^{(4)}(\eta)=16\,h^4\,K\]
It follows that:
\[I_{2h}-I_h=-15K\,h^4.\]
Now we can estimate the error at step \(h\):
\[E_h=h^4\,K=\frac{I_h-I_{2h}}{15}.\]
Based on the error estimate, we can compute a better approximation \(I_{h}^*\):
\[I_h^* = I_h + \frac{1}{15}\,(I_{h}-I_{2h})\]
or
\[I_h^* = \frac{16}{15}\,I_h - \frac{1}{15}\,I_{2h}\]
Numerical example#
We previously computed the integral with Simpson’s rule at step \(h=0.5\) and at step \(h=0.25\) already; the results were:
[I_Simps, I_Simps_comp]
\[\displaystyle \left[ 1.44083929301393, \ 1.44044796026391\right]\]
Using the formula above, we compute a better approximation:
I_Simps_better = 16/15*I_Simps_comp - 1/15*I_Simps
print(f'Exact result: {I_exact}\nBetter approximation: {I_Simps_better}')
Exact result: 1.4404224209802097
Better approximation: 1.4404218714139077
We obtain a better numerical approximation, but we lose the error estimate!
Simpson’s method in scipy.integrate.simps#
The composite Simpson’s rule is implemented in scipy in scipy.integrate.simpson() (documentation):
simpson(y, x=None, dx=1, axis=-1, even='avg')
where the parameters are:
ythe table of function values that we integrate,xthe nodes, an optional parameter; ifx=None, equidistant subintervals of widthdxare assumed,dxthe width of the equidistant subintervals, i.e., the integration step,axisthe integration axis (important in the case of a multidimensional numeric array)evencan be {'avg','first','last'} and defines the integration approach in the case of an odd number of subintervals.
Let us look at an example; first we import the simpson function:
from scipy.integrate import simpson
#%%timeit
simpson(y5v, dx=h5v)
\[\displaystyle 1.44044796026391\]
Romberg method#
The Romberg method is based on Richardson extrapolation. Let us assume that we integrate the integral \(\int_a^b f(x)\textrm{d}x\) on the interval \([a,b]\), which we divide into \(2^{n-1}\) subintervals (\(n=1,2,4,8, \dots\)).
We denote the result of the trapezoidal rule at \(n=1\) by \(R_{\underbrace{1}_{n},\underbrace{1}_{j}}\), where \(j\) denotes the accuracy of the obtained result \(\mathcal{O}(h^{2j})\).
If we apply the trapezoidal integration rule at \(n=1,2,4 \dots\) subintervals, we compute:
\(R_{1,1}\), step \(h_1=h\), order of accuracy \(\mathcal{O}(h_1^2)\),
\(R_{2,1}\), step \(h_2=h/2\), order of accuracy \(\mathcal{O}(h_2^2)\),
\(R_{3,1}\), step \(h_3=h/4\), order of accuracy \(\mathcal{O}(h_3^2)\),
\(R_{4,1}\), step \(h_4=h/8\), order of accuracy \(\mathcal{O}(h_4^2)\),
\(\dots\)
\(R_{n,1}\), step \(h_4=h/(2^{n-1})\), order of accuracy \(\mathcal{O}(h_n^2)\).
Based on Richardson extrapolation we compute a better approximation
\(R_{2,2} = R_{2,1} + \frac{1}{3}\left(R_{2,1}-R_{1,1}\right)\), step \(h_2=h/2\), order of accuracy \(\mathcal{O}(h_2^4)\),
\(R_{3,2} = R_{3,1} + \frac{1}{3}\left(R_{3,1}-R_{2,1}\right)\), step \(h_3=h/4\), order of accuracy \(\mathcal{O}(h_3^4)\),
\(\dots\)
\(R_{n,2} = R_{n,1} + \frac{1}{3}\left(R_{n,1}-R_{n-1,1}\right)\), step \(h_n=h/(2^{n-1})\), order of accuracy \(\mathcal{O}(h_n^4)\).
Then we continue with Richardson extrapolation:
\(R_{3,3} = R_{3,2} + \frac{1}{15}\left(R_{3,2}-R_{2,2}\right)\), step \(h_3=h/4\), order of accuracy \(\mathcal{O}(h_3^6)\),
\(\dots\)
\(R_{n,3} = R_{n,2} + \frac{1}{15}\left(R_{n,2}-R_{n-1,2}\right)\), step \(h_n=h/(2^{n-1})\), order of accuracy \(\mathcal{O}(h_n^6)\).
Richardson extrapolation can be generalized:
\(R_{n,j} = R_{n,j-1} + \frac{1}{4^{j-1}-1}\left(R_{n,j-1}-R_{n-1,j-1}\right)\), step \(h_n=h/(2^{n-1})\), order of accuracy \(\mathcal{O}(h_n^{2j})\)
Here, the better approximation \(R_{2,2}\) at step \(h/2\) is equal to the result obtained by Simpson’s method at step \(h/2\). Similarly, the better approximation \(R_{3,2}\) at step \(h/4\) is equal to the numerical integral of Simpson’s method at step \(h/4\). And further, \(R_{3,3}\) is equal to the corrected approximation of Simpson’s method at step \(h/4\).
Let us prepare the numerical data:
x9v, h9v = np.linspace(1, 2, 9, retstep=True) # step h9v = 0.125 (9 nodes)
y9v = x9v * np.sin(x9v)
Let us look at the example from above. First, we compute the integral with the composite trapezoidal method at different steps (second-order error):
## O(h^2)
R1_1 = trapezoidal_comp(y9v[::8], 8*h9v) # h=1.0
R2_1 = trapezoidal_comp(y9v[::4], 4*h9v) # h=0.5
R3_1 = trapezoidal_comp(y9v[::2], 2*h9v) # h=0.25
R4_1 = trapezoidal_comp(y9v, h9v) # h=0.125
[R1_1, R2_1, R3_1, R4_1]
\[\displaystyle \left[ 1.33003291922963, \ 1.41313769956786, \ 1.4336203950899, \ 1.43872310575291\right]\]
Then we compute better approximations (we obtain a fourth-order error):
## O(h^4)
R2_2 = R2_1 + 1/3 * (R2_1 - R1_1)
R3_2 = R3_1 + 1/3 * (R3_1 - R2_1)
R4_2 = R4_1 + 1/3 * (R4_1 - R3_1)
[R2_2, R3_2, R4_2]
\[\displaystyle \left[ 1.44083929301393, \ 1.44044796026391, \ 1.44042400930725\right]\]
The results represent the result of Simpson’s method at step \(h=0.5\), \(h=0.25\) and \(h=0.125\):
[simpson_comp(y9v[::4], 4*h9v), simpson_comp(y9v[::2], 2*h9v), simpson_comp(y9v, h9v)]
\[\displaystyle \left[ 1.44083929301393, \ 1.44044796026391, \ 1.44042400930725\right]\]
We again compute better approximations (we obtain a sixth-order error):
## O(h6)
R3_3 = R3_2 + 1/15 * (R3_2 - R2_2)
R4_3 = R4_2 + 1/15 * (R4_2 - R3_2)
[R3_3, R4_3]
\[\displaystyle \left[ 1.44042187141391, \ 1.44042241257681\right]\]
The result represents the better result of Simpson’s method at step \(h=0.25\) and \(h=0.125\):
a4h, a2h, ah = [simpson_comp(y9v[::4], 4*h9v), simpson_comp(y9v[::2], 2*h9v), simpson_comp(y9v, h9v)]
[16/15*a2h-1/15*a4h, 16/15*ah-1/15*a2h]
\[\displaystyle \left[ 1.44042187141391, \ 1.44042241257681\right]\]
We again compute better approximations (we obtain an eighth-order error):
## O(h8)
R4_4 = R4_3 + 1/63 * (R4_3 - R3_3)
R4_4
\[\displaystyle 1.4404224211667\]
The Romberg method thus offers high result accuracy at a relatively low numerical cost. The error is estimated by:
\[E = \left|R_{n,n}-R_{n-1,n-1}\right|.\]
Romberg method in scipy.integrate.romb#
The Romberg method is implemented in scipy in scipy.integrate.romb() (documentation):
romb(y, dx=1.0, axis=-1, show=False)
where the parameters are:
ythe table of function values that we integrate,dxthe width of the equidistant subintervals, i.e., the integration step,axisthe integration axis (important in the case of a multidimensional numeric array),showifTrue, displays the elements of the Richardson extrapolation.
Let us look at the example from above:
from scipy.integrate import romb
y9v
array([0.84147098, 1.01505104, 1.18623077, 1.34872795, 1.49624248,
1.62261343, 1.72197541, 1.78891084, 1.81859485])
romb(y9v, dx=h9v, show=True)
Richardson Extrapolation Table for Romberg Integration
======================================================
1.33003
1.41314 1.44084
1.43362 1.44045 1.44042
1.43872 1.44042 1.44042 1.44042
======================================================
\[\displaystyle 1.4404224211667\]
Gaussian integration approach#
The Newton-Cotes approach is based on a polynomial of degree \(n\), and the error is of degree \(n+1\). This means that the integration yields an exact result if the integrated function is a polynomial of degree \(n\) or less.
The Gaussian integration approach is fundamentally different: the goal is to replace the integral of the function \(f(x)\) with a weighted sum of function values at discrete values \(f(x_i)\):
\[\int_a^bf(x)\,\textrm{d}x\approx \sum_{i=0}^{n-1} A_i\, f(x_i).\]
Here both the weight \(A_i\) and the location of the node \(x_i\) are unknown. For polynomials of degree \(n\) we will also need more points \((x_i, f(x_i))\).
In what follows, we will learn that we can compute a numerically exact integral very efficiently. An additional advantage of Gaussian integration is that it can compute the integral of functions with singularities (e.g.: \(\int_0^1\sin(x)/\sqrt{(x)}\,dx\)).
Gaussian quadrature with one node#
Let us assume that we want to integrate a polynomial of degree \(n=1\) (a linear function):
\[f(x)=P_1(x)=a_0+a_1\,x.\]
Let us compute:
\[\int_{x_L}^{x_D}P_1(x)\,\textrm{d}x=\left(a_0\,x+a_1\,\frac{x^2}{2}\right)_{x_L}^{x_D}=-a_0\,x_L+a_0\,x_D-\frac{a_1\,x_L^2}{2}+\frac{a_1\,x_D^2}{2}.\]
On the other hand, we want to compute the integral from the appropriately weighted \(A_0\) value of the function \(f(x_0)\) at an unknown node \(x_0\) (only one node!):
\[\int_{x_L}^{x_D}P_1(x)\,\textrm{d}x = A_0\,P_1(x_0)= A_0\,a_0+A_0\,a_1\,x_0.\]
By equating the expressions above, we derive:
\[-a_0\,x_L+a_0\,x_D-\frac{a_1\,x_L^2}{2}+\frac{a_1\,x_D^2}{2}=A_0\,a_0+A_0\,a_1\,x_0.\]
\(a_0\) and \(a_1\) are the coefficients of the linear function, which can take arbitrary values, so it holds that:
\[a_0\,\left(-x_L+x_D-A_0\right)=0\qquad\textrm{and}\qquad a_1\left(-\frac{x_L^2}{2}+\frac{x_D^2}{2}-A_0\,x_0\right)=0.\]
This is a system of linear equations with the solution:
\[A_0= x_D-x_L, \qquad x_0=\frac{x_L+x_D}{2}.\]
If the integrated function is linear, we will compute the true value based only on the value at a single point!
So that the Gaussian integration approach is independent of the integration limits \(x_L\), \(x_D\), we introduce standard limits.
Standard limits: \(x_L=-1\) and \(x_D=1\)#
For generality, we transform the limits \(x\in[x_L, x_D]\) into \(\xi\in[-1, +1]\) using:
\[x=\frac{x_D+x_L}{2}+\frac{x_D-x_L}{2}\xi\]
and
\[\textrm{d}x=\frac{x_D-x_L}{2}\textrm{d}\xi.\]
The following holds:
\[\int_{x_L}^{x_D}f(x)\,\textrm{d}x=\int_{-1}^1 g\left(\xi\right)\,\textrm{d}\xi,\]
where:
\[g(\xi)=\frac{x_D-x_L}{2}\,f\left(\frac{x_L+x_d}{2}+\frac{x_D-x_L}{2}\xi\right).\]
In the case of standard limits, at one Gauss point the weight is \(A_0=2\) and \(x_0=0\) is the value at which we must evaluate the function \(f(x_0)\).
Computational derivation of the weights and the Gauss point#
Let us define the symbols and set up the equation:
a_0, a_1, x, x_L, x_D, A_0, x_0 = sym.symbols('a_0, a_1, x, x_L, x_D, A_0, x_0') # symbols
P1 = a_0 + a_1*x # linear polynomial
eq = sym.Eq(P1.integrate((x, x_L, x_D)).expand(), A_0*P1.subs(x, x_0)) # theoretical integral = estimate with weights
eq
\[\displaystyle a_{0} x_{D} - a_{0} x_{L} + \frac{a_{1} x_{D}^{2}}{2} - \frac{a_{1} x_{L}^{2}}{2} = A_{0} \left(a_{0} + a_{1} x_{0}\right)\]
Let us prepare two equations (for the first we assume \(a_0=0, a_1=1\), for the second we assume \(a_0=1, a_1=0\)) and solve the system for A_0 and x_0:
sym.solve([eq.subs(a_0, 0).subs(a_1, 1), eq.subs(a_1, 0).subs(a_0, 1)], [A_0, x_0])
\[\displaystyle \left[ \left( x_{D} - x_{L}, \ \frac{x_{D} + x_{L}}{2}\right)\right]\]
For additional explanation, I recommend this video.
Gaussian integration method with multiple nodes#
We want to derive a Gaussian integration method that will take into account \(v\) nodes on the interval \([a,b]\) and will exactly compute the integral of polynomials up to degree \(n=2v-1\). The following must hold:
\[\int_a^bP_{2v-1}(x)\,\textrm{d}x = \sum_{i=0}^{v-1} A_i\,P_{2v-1}(x_i).\]
In the derivation we will restrict ourselves to the standard limits \(x_L=-1\), \(x_D=1\),
where the polynomial of degree \(n=2v-1\) is defined as:
\[P_{2v-1}(x)=\sum_{i=0}^{2v-1} a_i\,x^i.\]
With two Gauss points/nodes we exactly compute the integral of a polynomial up to third order, and with three Gauss nodes we exactly compute the integral of a polynomial up to fifth order!
Computational derivation#
Let us first prepare the symbolic notation of the polynomial and the corresponding variables:
def P_etc(n=1, a='a', x='x'): # n is the degree of the polynomial
ai = sym.symbols(f'{a}:{n+1}') # list of a_i
x = sym.symbols(x) # variable x
return ai, x, sum([ai[i]*x**i for i in range(n+1)])
Now let us find the weights \(A_i\) and the nodes \(x_i\) for the case of two Gauss nodes; the polynomial is therefore of third degree.
v = 2 # number of nodes
ai, x, P = P_etc(n=2*v-1)
xi = sym.symbols(f'x:{v}') # list of x_i
Ai = sym.symbols(f'A:{v}') # list of w_i
print(f'Nodes: {xi}\nWeights: {Ai}')
Nodes: (x0, x1)
Weights: (A0, A1)
The polynomial:
P
\[\displaystyle a_{0} + a_{1} x + a_{2} x^{2} + a_{3} x^{3}\]
Similarly to the one node above, here let us define the equations:
eqs = [sym.Eq(P.integrate((x, -1, 1)).coeff(a_),\
sum([Ai[i]*P.subs(x, xi[i]) \
for i in range(v)]).expand().coeff(a_)) \
for a_ in ai]
eqs
\[\displaystyle \left[ 2 = A_{0} + A_{1}, \ 0 = A_{0} x_{0} + A_{1} x_{1}, \ \frac{2}{3} = A_{0} x_{0}^{2} + A_{1} x_{1}^{2}, \ 0 = A_{0} x_{0}^{3} + A_{1} x_{1}^{3}\right]\]
We solve them for the unknown \(x_i\) and \(w_i\):
sol = sym.solve(eqs, sym.flatten((xi, Ai)))
sol
\[\displaystyle \left[ \left( - \frac{\sqrt{3}}{3}, \ \frac{\sqrt{3}}{3}, \ 1, \ 1\right), \ \left( \frac{\sqrt{3}}{3}, \ - \frac{\sqrt{3}}{3}, \ 1, \ 1\right)\right]\]
We determined a list of two (equal) solutions.
First, the nodes are defined: \(x_0=-\sqrt{3}/3\) and \(x_1=\sqrt{3}/3\), with the corresponding weights \(A_0=A_1=1\).
The code above is derived in general - you can increase the number of nodes v and compute the nodes and the corresponding weights.
Below is a table of nodes and weights for one, two, and three nodes (for the limits \(a=-1\), \(b=1\)):
Number of points 1:
\(i\) |
Node \(x_i\) |
Weight \(A_i\) |
|---|---|---|
0 |
0 |
2 |
Number of points 2:
\(i\) |
Node \(x_i\) |
Weight \(A_i\) |
|---|---|---|
0 |
\(-\frac{\sqrt{3}}{3}\) |
1 |
1 |
\(+\frac{\sqrt{3}}{3}\) |
1 |
Number of points 3:
\(i\) |
Node \(x_i\) |
Weight \(A_i\) |
|---|---|---|
0 |
\(-\frac{\sqrt{15}}{5}\) |
\(\frac{5}{9}\) |
1 |
\(0\) |
\(\frac{8}{9}\) |
2 |
\(\frac{\sqrt{15}}{5}\) |
\(\frac{5}{9}\) |
For more nodes and also the error estimate, see Mathworld Legendre-Gauss Quadrature.
For the case of three Gauss points, we compute the numerical integral (standard limits):
\[I_{\textrm{Gauss3}}= \frac{5}{9}\,f\left(-\frac{\sqrt{15}}{5}\right) + \frac{8}{9}\,f\left(0\right)+\frac{5}{9}\,f\left(\frac{\sqrt{15}}{5}\right).\]
Numerical implementation#
The numerical implementation (including the transformation of the limits) for one, two, or three nodes:
def gaussian(fun, a, b, nodes=1):
"""
Gaussian integration method.
:param fun: object of the function to integrate
:param a: lower limit
:param b: upper limit
:param nodes: number of nodes (1, 2, or 3)
"""
def g(xi): # function for transforming the limits
return (b-a)/2 * fun((b+a +xi * (b-a))/2)
if nodes == 1:
return 2*g(0.)
elif nodes == 2:
return 1. * g(-np.sqrt(3)/3) + 1. * g(np.sqrt(3)/3)
elif nodes == 3:
return 5/9 * g(-np.sqrt(15)/5) +8/9 * g(0.) + 5/9 * g(np.sqrt(15)/5)
Let us look at an example. First, let us define the function we want to integrate:
def f(x):
return x*np.sin(x)
Now we pass the function (in this concrete case f) and not the value (e.g., f(0.)) to the gaussian function. First for one node, then two and three:
gaussian(fun=f, a=1., b=2., nodes=1)
\[\displaystyle 1.49624247990608\]
gaussian(fun=f, a=1., b=2., nodes=2)
\[\displaystyle 1.44014401845179\]
gaussian(fun=f, a=1., b=2., nodes=3)
\[\displaystyle 1.44042294912151\]
scipy.integrate.fixed_quad#
Within scipy, the Gaussian integration method is implemented in the scipy.integrate.fixed_quad() function (documentation):
fixed_quad(func, a, b, args=(), n=5)
where the parameters are:
funcis the name of the function we call,ais the lower limit,bis the upper limit,argsis a tuple of any additional arguments of the functionfunc,nis the number of nodes of the Gaussian integration, defaultn=5.
The function returns a tuple with the integration result val and the value None: (val, None)
Let us look at the example from above:
from scipy.integrate import fixed_quad
fixed_quad(f, a=1, b=2, n=3)[0]
\[\displaystyle 1.44042294912151\]
The result is equal to the previous one:
gaussian(fun=f, a=1., b=2., nodes=3)
\[\displaystyle 1.44042294912151\]
scipy.integrate#
scipy.integrate is a powerful tool for numerical integration (see the documentation). In what follows, we will look at some of the functions.
Integration functions that require a defining function:#
quad(func, a, b[, args, full_output, ...])computes the definite integral offunc(x)over the limits [a,b],dblquad(func, a, b, gfun, hfun[, args, ...])computes the definite integral offunc(x,y),tplquad(func, a, b, gfun, hfun, qfun, rfun)computes the definite integral offunc(x,y,z),nquad(func, ranges[, args, opts, full_output])computes the definite integral of the \(n\)-dimensional functionfunc( ...),romberg(function, a, b[, args, tol, rtol, ...])Romberg integration for the functionfunction,
scipy.integrate.quad#
Let us look at a very useful function for integration, scipy.integrate.quad() (documentation):
quad(func, a, b, args=(), full_output=0, epsabs=1.49e-08, epsrel=1.49e-08, limit=50,
points=None, weight=None, wvar=None, wopts=None, maxp1=50, limlst=50)
The selected parameters are:
functhe Python function we integrate,athe lower limit of integration (np.infcan be used for an infinite limit),bthe upper limit of integration (np.infcan be used for an infinite limit),full_outputto display all results, default 0,epsabsthe allowed absolute error,epsrelthe allowed relative error.
Let us look at the example from above:
from scipy import integrate
integrate.quad(f, a=1, b=2)
\[\displaystyle \left( 1.44042242098021, \ 1.59919013695854 \cdot 10^{-14}\right)\]
scipy.integrate.dblquad#
This is a function similar to quad, but for integration over two variables (documentation).
The function scipy.integrate.dblquad() computes the double integral of func(y, x) over the limits from x = [a, b] and y = [gfun(x), hfun(x)].
dblquad(func, a, b, gfun, hfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)
The selected parameters are:
funcis the Python function we integrate,ais the lower limit of integration ofx(np.infcan be used for an infinite limit),bis the upper limit of integration ofx(np.infcan be used for an infinite limit),gfunis the Python function that defines the lower limit ofyas a function ofx,hfunis the Python function that defines the upper limit ofyas a function ofx.
Let us look at an example of computing the area of a half-circle with radius 1:
\[\int_{-1}^{1}\left(\int_0^{\sqrt{1-x^2}}1\,\textrm{d}y\right)\,\textrm{d}x=\frac{\pi}{2}\]
Let us define the appropriate Python functions and compute the result:
def func(y, x): # the integration function is simple, a constant = 1!
return 1.
def gfun(x): # lower limit = 0
return 0.
def hfun(x): # upper limit
return np.sqrt(1-x**2)
integrate.dblquad(func=func, a=-1, b=1, gfun=gfun, hfun=hfun)
\[\displaystyle \left( 1.5707963267949, \ 1.00023545002159 \cdot 10^{-9}\right)\]
Integration functions that require a table of values:#
trapz(y[, x, dx, axis])composite trapezoidal rule,cumulative_trapezoid(y[, x, dx, axis, initial])cumulative integral of the integrand (returns the result at each node),simps(y[, x, dx, axis, even])Simpson’s method,romb(y[, dx, axis, show])Romberg method.
Here we will look at the cumulative_trapezoid function on an example. As an example, consider a mass \(m=1\) (we omit the units) acted upon by a force \(F=\sin(t)\). From Newton’s 2nd law it follows that: \(F=m\,\ddot x\). We integrate the acceleration to compute the velocity; then we integrate once more for the displacement. Let us first define the function for the acceleration:
def acceleration(t):
F = np.sin(t)
m = 1
return F/m
We are interested in what happens over a time of 3 seconds:
t, h = np.linspace(0, 3, 100, retstep=True)
Let us compute the table of accelerations and then integrate for the velocity (here it is important to define the initial value):
a = acceleration(t)
v = integrate.cumulative_trapezoid(y=a, dx=h, initial=0)
Now let us integrate the velocity once more to compute the path:
s = integrate.cumulative_trapezoid(y=v, dx=h, initial=0)
Let us display the result:
plt.plot(t, a, label='Acceleration')
plt.plot(t, v, label='Velocity')
plt.plot(t, s, label='Path')
plt.xlabel('$t$ [s]')
plt.legend();
Error in callback <function _draw_all_if_interactive at 0x7fae1b1a9c60> (for post_execute), with arguments args (),kwargs {}:
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:250, in TexManager._run_checked_subprocess(cls, command, tex, cwd)
249 try:
--> 250 report = subprocess.check_output(
251 command, cwd=cwd if cwd is not None else cls._texcache,
252 stderr=subprocess.STDOUT)
253 except FileNotFoundError as exc:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:466, in check_output(timeout, *popenargs, **kwargs)
464 kwargs['input'] = empty
--> 466 return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
467 **kwargs).stdout
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:548, in run(input, capture_output, timeout, check, *popenargs, **kwargs)
546 kwargs['stderr'] = PIPE
--> 548 with Popen(*popenargs, **kwargs) as process:
549 try:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:1026, in Popen.__init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, user, group, extra_groups, encoding, errors, text, umask, pipesize, process_group)
1023 self.stderr = io.TextIOWrapper(self.stderr,
1024 encoding=encoding, errors=errors)
-> 1026 self._execute_child(args, executable, preexec_fn, close_fds,
1027 pass_fds, cwd, env,
1028 startupinfo, creationflags, shell,
1029 p2cread, p2cwrite,
1030 c2pread, c2pwrite,
1031 errread, errwrite,
1032 restore_signals,
1033 gid, gids, uid, umask,
1034 start_new_session, process_group)
1035 except:
1036 # Cleanup if the child failed starting.
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:1955, in Popen._execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, gid, gids, uid, umask, start_new_session, process_group)
1954 if err_filename is not None:
-> 1955 raise child_exception_type(errno_num, err_msg, err_filename)
1956 else:
FileNotFoundError: [Errno 2] No such file or directory: 'latex'
The above exception was the direct cause of the following exception:
RuntimeError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/pyplot.py:278, in _draw_all_if_interactive()
276 def _draw_all_if_interactive() -> None:
277 if matplotlib.is_interactive():
--> 278 draw_all()
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/_pylab_helpers.py:131, in Gcf.draw_all(cls, force)
129 for manager in cls.get_all_fig_managers():
130 if force or manager.canvas.figure.stale:
--> 131 manager.canvas.draw_idle()
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backend_bases.py:1893, in FigureCanvasBase.draw_idle(self, *args, **kwargs)
1891 if not self._is_idle_drawing:
1892 with self._idle_draw_cntx():
-> 1893 self.draw(*args, **kwargs)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backends/backend_agg.py:382, in FigureCanvasAgg.draw(self)
379 # Acquire a lock on the shared font cache.
380 with (self.toolbar._wait_cursor_for_draw_cm() if self.toolbar
381 else nullcontext()):
--> 382 self.figure.draw(self.renderer)
383 # A GUI class may be need to update a window using this draw, so
384 # don't forget to call the superclass.
385 super().draw()
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:94, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
92 @wraps(draw)
93 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 94 result = draw(artist, renderer, *args, **kwargs)
95 if renderer._rasterizing:
96 renderer.stop_rasterizing()
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/figure.py:3264, in Figure.draw(self, renderer)
3261 # ValueError can occur when resizing a window.
3263 self.patch.draw(renderer)
-> 3264 mimage._draw_list_compositing_images(
3265 renderer, self, artists, self.suppressComposite)
3267 renderer.close_group('figure')
3268 finally:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axes/_base.py:3226, in _AxesBase.draw(self, renderer)
3223 if artists_rasterized:
3224 _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3226 mimage._draw_list_compositing_images(
3227 renderer, self, artists, self.get_figure(root=True).suppressComposite)
3229 renderer.close_group('axes')
3230 self.stale = False
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1405, in Axis.draw(self, renderer)
1402 renderer.open_group(__name__, gid=self.get_gid())
1404 ticks_to_draw = self._update_ticks()
-> 1405 tlb1, tlb2 = self._get_ticklabel_bboxes(ticks_to_draw, renderer)
1407 for tick in ticks_to_draw:
1408 tick.draw(renderer)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1332, in Axis._get_ticklabel_bboxes(self, ticks, renderer)
1330 if renderer is None:
1331 renderer = self.get_figure(root=True)._get_renderer()
-> 1332 return ([tick.label1.get_window_extent(renderer)
1333 for tick in ticks if tick.label1.get_visible()],
1334 [tick.label2.get_window_extent(renderer)
1335 for tick in ticks if tick.label2.get_visible()])
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1332, in <listcomp>(.0)
1330 if renderer is None:
1331 renderer = self.get_figure(root=True)._get_renderer()
-> 1332 return ([tick.label1.get_window_extent(renderer)
1333 for tick in ticks if tick.label1.get_visible()],
1334 [tick.label2.get_window_extent(renderer)
1335 for tick in ticks if tick.label2.get_visible()])
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:969, in Text.get_window_extent(self, renderer, dpi)
964 raise RuntimeError(
965 "Cannot get window extent of text w/o renderer. You likely "
966 "want to call 'figure.draw_without_rendering()' first.")
968 with cbook._setattr_cm(fig, dpi=dpi):
--> 969 bbox, info, descent = self._get_layout(self._renderer)
970 x, y = self.get_unitless_position()
971 x, y = self.get_transform().transform((x, y))
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:373, in Text._get_layout(self, renderer)
370 ys = []
372 # Full vertical extent of font, including ascenders and descenders:
--> 373 _, lp_h, lp_d = _get_text_metrics_with_cache(
374 renderer, "lp", self._fontproperties,
375 ismath="TeX" if self.get_usetex() else False,
376 dpi=self.get_figure(root=True).dpi)
377 min_dy = (lp_h - lp_d) * self._linespacing
379 for i, line in enumerate(lines):
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:69, in _get_text_metrics_with_cache(renderer, text, fontprop, ismath, dpi)
66 """Call ``renderer.get_text_width_height_descent``, caching the results."""
67 # Cached based on a copy of fontprop so that later in-place mutations of
68 # the passed-in argument do not mess up the cache.
---> 69 return _get_text_metrics_with_cache_impl(
70 weakref.ref(renderer), text, fontprop.copy(), ismath, dpi)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:77, in _get_text_metrics_with_cache_impl(renderer_ref, text, fontprop, ismath, dpi)
73 @functools.lru_cache(4096)
74 def _get_text_metrics_with_cache_impl(
75 renderer_ref, text, fontprop, ismath, dpi):
76 # dpi is unused, but participates in cache invalidation (via the renderer).
---> 77 return renderer_ref().get_text_width_height_descent(text, fontprop, ismath)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backends/backend_agg.py:211, in RendererAgg.get_text_width_height_descent(self, s, prop, ismath)
209 _api.check_in_list(["TeX", True, False], ismath=ismath)
210 if ismath == "TeX":
--> 211 return super().get_text_width_height_descent(s, prop, ismath)
213 if ismath:
214 ox, oy, width, height, descent, font_image = \
215 self.mathtext_parser.parse(s, self.dpi, prop)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backend_bases.py:566, in RendererBase.get_text_width_height_descent(self, s, prop, ismath)
562 fontsize = prop.get_size_in_points()
564 if ismath == 'TeX':
565 # todo: handle properties
--> 566 return self.get_texmanager().get_text_width_height_descent(
567 s, fontsize, renderer=self)
569 dpi = self.points_to_pixels(72)
570 if ismath:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:364, in TexManager.get_text_width_height_descent(cls, tex, fontsize, renderer)
362 if tex.strip() == '':
363 return 0, 0, 0
--> 364 dvifile = cls.make_dvi(tex, fontsize)
365 dpi_fraction = renderer.points_to_pixels(1.) if renderer else 1
366 with dviread.Dvi(dvifile, 72 * dpi_fraction) as dvi:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:293, in TexManager.make_dvi(cls, tex, fontsize)
290 with TemporaryDirectory(dir=dvifile.parent) as tmpdir:
291 Path(tmpdir, "file.tex").write_text(
292 cls._get_tex_source(tex, fontsize), encoding='utf-8')
--> 293 cls._run_checked_subprocess(
294 ["latex", "-interaction=nonstopmode", "-halt-on-error",
295 "-no-shell-escape", "file.tex"], tex, cwd=tmpdir)
296 Path(tmpdir, "file.dvi").replace(dvifile)
297 # Also move the tex source to the main cache directory, but
298 # only for backcompat.
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:254, in TexManager._run_checked_subprocess(cls, command, tex, cwd)
250 report = subprocess.check_output(
251 command, cwd=cwd if cwd is not None else cls._texcache,
252 stderr=subprocess.STDOUT)
253 except FileNotFoundError as exc:
--> 254 raise RuntimeError(
255 f'Failed to process string with tex because {command[0]} '
256 'could not be found') from exc
257 except subprocess.CalledProcessError as exc:
258 raise RuntimeError(
259 '{prog} was not able to process the following string:\n'
260 '{tex!r}\n\n'
(...) 267 exc=exc.output.decode('utf-8', 'backslashreplace'))
268 ) from None
RuntimeError: Failed to process string with tex because latex could not be found
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:250, in TexManager._run_checked_subprocess(cls, command, tex, cwd)
249 try:
--> 250 report = subprocess.check_output(
251 command, cwd=cwd if cwd is not None else cls._texcache,
252 stderr=subprocess.STDOUT)
253 except FileNotFoundError as exc:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:466, in check_output(timeout, *popenargs, **kwargs)
464 kwargs['input'] = empty
--> 466 return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
467 **kwargs).stdout
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:548, in run(input, capture_output, timeout, check, *popenargs, **kwargs)
546 kwargs['stderr'] = PIPE
--> 548 with Popen(*popenargs, **kwargs) as process:
549 try:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:1026, in Popen.__init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, user, group, extra_groups, encoding, errors, text, umask, pipesize, process_group)
1023 self.stderr = io.TextIOWrapper(self.stderr,
1024 encoding=encoding, errors=errors)
-> 1026 self._execute_child(args, executable, preexec_fn, close_fds,
1027 pass_fds, cwd, env,
1028 startupinfo, creationflags, shell,
1029 p2cread, p2cwrite,
1030 c2pread, c2pwrite,
1031 errread, errwrite,
1032 restore_signals,
1033 gid, gids, uid, umask,
1034 start_new_session, process_group)
1035 except:
1036 # Cleanup if the child failed starting.
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/subprocess.py:1955, in Popen._execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, gid, gids, uid, umask, start_new_session, process_group)
1954 if err_filename is not None:
-> 1955 raise child_exception_type(errno_num, err_msg, err_filename)
1956 else:
FileNotFoundError: [Errno 2] No such file or directory: 'latex'
The above exception was the direct cause of the following exception:
RuntimeError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backend_bases.py:2157, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
2154 # we do this instead of `self.figure.draw_without_rendering`
2155 # so that we can inject the orientation
2156 with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2157 self.figure.draw(renderer)
2158 if bbox_inches:
2159 if bbox_inches == "tight":
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:94, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
92 @wraps(draw)
93 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 94 result = draw(artist, renderer, *args, **kwargs)
95 if renderer._rasterizing:
96 renderer.stop_rasterizing()
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/figure.py:3264, in Figure.draw(self, renderer)
3261 # ValueError can occur when resizing a window.
3263 self.patch.draw(renderer)
-> 3264 mimage._draw_list_compositing_images(
3265 renderer, self, artists, self.suppressComposite)
3267 renderer.close_group('figure')
3268 finally:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axes/_base.py:3226, in _AxesBase.draw(self, renderer)
3223 if artists_rasterized:
3224 _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3226 mimage._draw_list_compositing_images(
3227 renderer, self, artists, self.get_figure(root=True).suppressComposite)
3229 renderer.close_group('axes')
3230 self.stale = False
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1405, in Axis.draw(self, renderer)
1402 renderer.open_group(__name__, gid=self.get_gid())
1404 ticks_to_draw = self._update_ticks()
-> 1405 tlb1, tlb2 = self._get_ticklabel_bboxes(ticks_to_draw, renderer)
1407 for tick in ticks_to_draw:
1408 tick.draw(renderer)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1332, in Axis._get_ticklabel_bboxes(self, ticks, renderer)
1330 if renderer is None:
1331 renderer = self.get_figure(root=True)._get_renderer()
-> 1332 return ([tick.label1.get_window_extent(renderer)
1333 for tick in ticks if tick.label1.get_visible()],
1334 [tick.label2.get_window_extent(renderer)
1335 for tick in ticks if tick.label2.get_visible()])
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/axis.py:1332, in <listcomp>(.0)
1330 if renderer is None:
1331 renderer = self.get_figure(root=True)._get_renderer()
-> 1332 return ([tick.label1.get_window_extent(renderer)
1333 for tick in ticks if tick.label1.get_visible()],
1334 [tick.label2.get_window_extent(renderer)
1335 for tick in ticks if tick.label2.get_visible()])
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:969, in Text.get_window_extent(self, renderer, dpi)
964 raise RuntimeError(
965 "Cannot get window extent of text w/o renderer. You likely "
966 "want to call 'figure.draw_without_rendering()' first.")
968 with cbook._setattr_cm(fig, dpi=dpi):
--> 969 bbox, info, descent = self._get_layout(self._renderer)
970 x, y = self.get_unitless_position()
971 x, y = self.get_transform().transform((x, y))
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:373, in Text._get_layout(self, renderer)
370 ys = []
372 # Full vertical extent of font, including ascenders and descenders:
--> 373 _, lp_h, lp_d = _get_text_metrics_with_cache(
374 renderer, "lp", self._fontproperties,
375 ismath="TeX" if self.get_usetex() else False,
376 dpi=self.get_figure(root=True).dpi)
377 min_dy = (lp_h - lp_d) * self._linespacing
379 for i, line in enumerate(lines):
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:69, in _get_text_metrics_with_cache(renderer, text, fontprop, ismath, dpi)
66 """Call ``renderer.get_text_width_height_descent``, caching the results."""
67 # Cached based on a copy of fontprop so that later in-place mutations of
68 # the passed-in argument do not mess up the cache.
---> 69 return _get_text_metrics_with_cache_impl(
70 weakref.ref(renderer), text, fontprop.copy(), ismath, dpi)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/text.py:77, in _get_text_metrics_with_cache_impl(renderer_ref, text, fontprop, ismath, dpi)
73 @functools.lru_cache(4096)
74 def _get_text_metrics_with_cache_impl(
75 renderer_ref, text, fontprop, ismath, dpi):
76 # dpi is unused, but participates in cache invalidation (via the renderer).
---> 77 return renderer_ref().get_text_width_height_descent(text, fontprop, ismath)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backends/backend_agg.py:211, in RendererAgg.get_text_width_height_descent(self, s, prop, ismath)
209 _api.check_in_list(["TeX", True, False], ismath=ismath)
210 if ismath == "TeX":
--> 211 return super().get_text_width_height_descent(s, prop, ismath)
213 if ismath:
214 ox, oy, width, height, descent, font_image = \
215 self.mathtext_parser.parse(s, self.dpi, prop)
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/backend_bases.py:566, in RendererBase.get_text_width_height_descent(self, s, prop, ismath)
562 fontsize = prop.get_size_in_points()
564 if ismath == 'TeX':
565 # todo: handle properties
--> 566 return self.get_texmanager().get_text_width_height_descent(
567 s, fontsize, renderer=self)
569 dpi = self.points_to_pixels(72)
570 if ismath:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:364, in TexManager.get_text_width_height_descent(cls, tex, fontsize, renderer)
362 if tex.strip() == '':
363 return 0, 0, 0
--> 364 dvifile = cls.make_dvi(tex, fontsize)
365 dpi_fraction = renderer.points_to_pixels(1.) if renderer else 1
366 with dviread.Dvi(dvifile, 72 * dpi_fraction) as dvi:
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:293, in TexManager.make_dvi(cls, tex, fontsize)
290 with TemporaryDirectory(dir=dvifile.parent) as tmpdir:
291 Path(tmpdir, "file.tex").write_text(
292 cls._get_tex_source(tex, fontsize), encoding='utf-8')
--> 293 cls._run_checked_subprocess(
294 ["latex", "-interaction=nonstopmode", "-halt-on-error",
295 "-no-shell-escape", "file.tex"], tex, cwd=tmpdir)
296 Path(tmpdir, "file.dvi").replace(dvifile)
297 # Also move the tex source to the main cache directory, but
298 # only for backcompat.
File /opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/matplotlib/texmanager.py:254, in TexManager._run_checked_subprocess(cls, command, tex, cwd)
250 report = subprocess.check_output(
251 command, cwd=cwd if cwd is not None else cls._texcache,
252 stderr=subprocess.STDOUT)
253 except FileNotFoundError as exc:
--> 254 raise RuntimeError(
255 f'Failed to process string with tex because {command[0]} '
256 'could not be found') from exc
257 except subprocess.CalledProcessError as exc:
258 raise RuntimeError(
259 '{prog} was not able to process the following string:\n'
260 '{tex!r}\n\n'
(...) 267 exc=exc.output.decode('utf-8', 'backslashreplace'))
268 ) from None
RuntimeError: Failed to process string with tex because latex could not be found
<Figure size 640x480 with 1 Axes>
Practice questions#
Question 1: Using the tools of the sympy package, symbolically determine the value of the first moment of area of a cross-section in the shape of a quarter circle shown in the figure:
\[S_{xx} = \int_{A}y~\text{d}A=\int_0^r \int_0^{\pi/2}y(\varphi)~r~\text{d}\varphi~\text{d}r = \int_{0}^{\pi/2}\frac{r^3}{3}\, \sin (\varphi)\, \text{d}\varphi\]
Also prepare a numerical function for the dependence of the integrand \(f(\varphi) = \Big(\frac{r^3}{3}\, \sin (\varphi) \Big)\) for a given value of the radius \(R\).
Data:
\(R = 7.5\) mm
Question 2: Use the numerical function prepared in the previous problem and compute and graphically display the values of the integrand from the first problem \(\Big(f(\varphi) = \frac{r^3}{3}\sin (\varphi) \Big)\) at 100 discrete points on the interval \(\varphi \in [0, \pi/2]\).
Also compute the value of the definite integral \(S\) numerically, using the trapezoidal and Simpson’s 1/3 methods from the numpy or scipy packages.
\[ S = \int_{0}^{\pi/2}\frac{r^3}{3}\sin (\varphi)\text{d}\varphi\]
Trapezoidal rule#
Question 3: Using your own implementation of the trapezoidal rule (the basic one, not the composite one), compute the definite integral
\[I = \int_{0}^{2} -0.5x^3 - x^2 + 8 ~\text{d}x\]
Compare the obtained value with the result of the numpy.trapz function (or scipy.integrate.trapezoid), where you divide the observed interval into 10 equidistant segments.
Improved approximation with Richardson extrapolation:#
\[\bar{I} = \frac{4I_h - I_{2h}}{3}\]
Question 4 : Use Richardson extrapolation and determine a more accurate value of the integral from the previous problem by first dividing the observed interval \(x \in [0, 2]\) into 2, then into 4 equidistant segments and using the composite trapezoidal method from the numpy package.
Simpson’s 1/3 rule#
Question 5: Compute the definite integral from the previous problem with your own code of Simpson’s 1/3 method (the basic one, not the composite one).
\[I = \int_{0}^{2} -0.5x^3 - x^2 + 8 ~\text{d}x\]
Compare the obtained value with the result of the scipy.integrate.simpson function, where you divide the observed interval into 7 equidistant points.
Question 6: The measured values of a car’s acceleration in the direction of travel over time are given. Using the cumulative trapezoidal method (scipy.integrate.cumulative_trapezoid), determine the values of the velocity and the distance traveled by the car over the observed time interval.
At time \(t = 0\) the car is at rest (hint: the initial argument).
#Data
t = np.linspace(0, 7, 10) # s
a = np.array([ 3.47, 4.03, 4.36, 4.48, 4.27, 3.78, 3.19, 2.73, 2.54, 2.53]) # m/s^2
Gaussian integration method (Gaussian quadrature)#
Question 7: Using the composite trapezoidal method from scipy and the Gaussian integration method (scipy.integrate.quad), compute the values of the integrals:
\[I_1 = \int_{0}^{1}\frac{\sin(x)}{x} \text{d}x\]
\[I_2 = \int_{1}^{2}\frac{\sin(x)}{x} \text{d}x\]
\[I_3 = \int_{-0.5}^{0.5}\frac{\sin(x)}{x} \text{d}x\]
When using the trapezoidal method, divide the observed interval into 10 points.
Question 8: Compute the integral:
\[I = \int_{-0.5}^{0.5}\frac{\sin(x)}{x} \text{d}x\]
using the scipy.integrate.fixed_quad function as well, to which you can also pass the parameter for the number of nodes in the integration, n. First use integration with one node, then with two, and finally with three nodes.