1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
|
import openturns as ot
from matplotlib import pyplot as plt
from openturns.viewer import View
import math as m
oldPrecision = ot.PlatformInfo.GetNumericalPrecision()
ot.PlatformInfo.SetNumericalPrecision(16)
a = -m.pi
b = m.pi
f = ot.SymbolicFunction(["x", "y"], ["1+cos(x)*sin(y)"])
ll = [ot.SymbolicFunction(["x"], [" 2+cos(x)"])]
u = [ot.SymbolicFunction(["x"], ["-2-cos(x)"])]
# Draw the graph of the integrand and the bounds:
g = ot.Graph("IteratedQuadrature example", "x", "y", True, "upper right")
g.add(f.draw([a, a], [b, b]))
curve = ll[0].draw(a, b).getDrawable(0)
curve.setLineWidth(2)
curve.setColor("red")
g.add(curve)
curve = u[0].draw(a, b).getDrawable(0)
curve.setLineWidth(2)
curve.setColor("red")
g.add(curve)
# Evaluate the integral with high precision:
Iref = ot.IteratedQuadrature(
ot.GaussKronrod(100000, 1e-13, ot.GaussKronrodRule(ot.GaussKronrodRule.G11K23))
).integrate(f, a, b, ll, u)
# Evaluate the integral with the default GaussKronrod algorithm:
f = ot.MemoizeFunction(f)
I1 = ot.IteratedQuadrature(ot.GaussKronrod()).integrate(f, a, b, ll, u)
sample1 = f.getInputHistory()
print(
"I1=",
I1,
"#evals=",
sample1.getSize(),
"err=",
abs(100.0 * (1.0 - I1[0] / Iref[0])),
"%",
)
cloud = ot.Cloud(sample1)
cloud.setPointStyle("fcircle")
cloud.setColor("green")
g.add(cloud)
f.clearHistory()
# Evaluate the integral with the default IteratedQuadrature algorithm:
I2 = ot.IteratedQuadrature().integrate(f, a, b, ll, u)
sample2 = f.getInputHistory()
# print('I2=', I2, '#evals=', sample2.getSize(), \
# 'err=', abs(100.0*(1.0-I2[0]/Iref[0])), '%')
cloud = ot.Cloud(sample2)
cloud.setPointStyle("fcircle")
cloud.setColor("gold")
g.add(cloud)
fig = plt.figure(figsize=(4, 4))
axis = fig.add_subplot(111)
axis.set_xlim(auto=True)
View(g, figure=fig, axes=[axis], add_legend=False)
ot.PlatformInfo.SetNumericalPrecision(oldPrecision)
|
