File: plot_taylor_approximation.py

package info (click to toggle)
openturns 1.26-4
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 67,708 kB
  • sloc: cpp: 261,605; python: 67,030; ansic: 4,378; javascript: 406; sh: 185; xml: 164; makefile: 101
file content (135 lines) | stat: -rw-r--r-- 3,445 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
"""
Taylor approximations
=====================
"""

# %%
# In this example we build a local approximation of a model using the
# Taylor decomposition with the :class:`~openturns.LinearTaylor` class.
#
# We consider the function :math:`\vect{h} : \Rset^2 \rightarrow \Rset^2` defined by:
#
# .. math::
#    \vect{h}(\vect{x}) = \Tr{\left( \cos(x_1 + x_2), (x_2 + 1) e^{x_1 - 2 x_2} \right)}.
#
# for any :math:`\vect{x} \in \Rset^2`.
# The metamodel :math:`\widehat{\vect{h}}` is an approximation of the model :math:`\vect{h}`:
#
# .. math::
#    \vect{y} \approx \widehat{\vect{h}}(\vect{x})
#
# for any :math:`\vect{x} \in \Rset^2`. We note  :math:`h_k`, :math:`k=1, 2` the marginal outputs
# In this example, we consider two different approximations:
#
# - the first order Taylor expansion,
# - the second order Taylor expansion.

# %%
import openturns as ot
import openturns.viewer as otv


# %%
# Define the model
# ~~~~~~~~~~~~~~~~

# %%
# Prepare some data.
formulas = ["cos(x1 + x2)", "(x2 + 1) * exp(x1 - 2 * x2)"]
model = ot.SymbolicFunction(["x1", "x2"], formulas)

# %%
# Center of the approximation.
x0 = [-0.4, -0.4]

# %%
# Drawing bounds.
a = -0.4
b = 0.0

# %%
# First order Taylor expansion
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# %%
# Let :math:`\vect{x}_0 \in \Rset^2` be a reference point where
# the linear approximation is evaluated.
# The first order Taylor expansion of the output marginal :math:`h_k` is:
#
# .. math::
#    \widehat{h_k}(\vect{x}) = h_k(\vect{x}_0) +
#      \sum_{i=1}^{2} \frac{\partial h_k}{\partial x_i}(\vect{x}_0) \left(x_i - x_{0,i} \right)
#
# for any :math:`\vect{x} \in \Rset^2`.


# %%
# Create a linear (first-order) Taylor approximation.
algo = ot.LinearTaylor(x0, model)
algo.run()
responseSurface = algo.getMetaModel()

# %%
# Plot the second output of our model with :math:`x_1=x_{0,1}`.

# %%
graph = ot.ParametricFunction(model, [0], [x0[1]]).getMarginal(1).draw(a, b)
graph.setLegends(["Model"])
curve = (
    ot.ParametricFunction(responseSurface, [0], [x0[1]])
    .getMarginal(1)
    .draw(a, b)
    .getDrawable(0)
)
curve.setLegend("Taylor")
curve.setLineStyle("dashed")
graph.add(curve)
graph.setLegendPosition("upper right")
view = otv.View(graph)

# %%
# Second order Taylor expansion
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# %%
# Let :math:`\vect{x}_0 \in \Rset^2` be a reference point where
# the quadratic approximation is evaluated.
# The second order Taylor expansion of the output marginal :math:`h_k` is:
#
# .. math::
#
#     \widehat{h_k}(\vect{x}) = h_k(\vect{x}_0) + \sum_{i=1}^{2}
#     \frac{\partial h_k}{\partial x_i}(\vect{x}_0) \left(x_i - x_{0,i} \right) +
#     \frac{1}{2} \sum_{i,j=1}^{2}
#     \frac{\partial^2 h_k}{\partial x_i \partial x_j}(\vect{x}_0) \left(x_i - x_{0,i} \right) \left(x_j - x_{0,j} \right)
#
# for any :math:`\vect{x} \in \Rset^2`.

# %%
# Create a quadratic (second-order) Taylor approximation.

# %%
algo = ot.QuadraticTaylor(x0, model)
algo.run()
responseSurface = algo.getMetaModel()

# %%
# Plot second output of our model with :math:`x_1=x_{0,1}`.

# %%
graph = ot.ParametricFunction(model, [0], [x0[1]]).getMarginal(1).draw(a, b)
graph.setLegends(["Model"])
curve = (
    ot.ParametricFunction(responseSurface, [0], [x0[1]])
    .getMarginal(1)
    .draw(a, b)
    .getDrawable(0)
)
curve.setLegend("Taylor")
curve.setLineStyle("dashed")
graph.add(curve)
graph.setLegendPosition("upper right")
view = otv.View(graph)

# %%
view.ShowAll()