Week 7: Riemann Integrals in 1D and 2D#
Demo by Karl Johan Funch Måstrup Kristensen, Magnus Troen, and Jakob Lemvig
from sympy import *
from dtumathtools import *
init_printing()
You are now about to do integration. And in many circumstances, SymPy will be of good help. SymPy’s \(\verb|integrate()|\) function takes care of many integrals, but you will often experience moments when \(\verb|integrate()|\) has difficulties with an integral, in which case you must use your own knowledge from the course to make the expression more “edible” for SymPy. Sometimes SymPy needs it cut into smaller pieces to not suffocate.
Integration in One Dimension#
Very simply, anti-derivatives are found using SymPy with \(\verb|integrate(f, var)|\). \(\verb|var|\) can be ledt out in cases where the function only has one variable, but as soon as multiple variables appear (or symbols defined with \(\verb|symbols()|\)) one must define which variable to integrate to.
Consider for example a function \(f:\mathbb{R} \to \mathbb{R}\) given by
We find an anti-derivative
x = symbols("x", real=True)
f = sqrt(1 - x**2)
f
![../_images/54b4a481594fb25b01401b043e003d74d98d069988c257cbc8b8e718c39c0535.png](../_images/54b4a481594fb25b01401b043e003d74d98d069988c257cbc8b8e718c39c0535.png)
F = integrate(f, x) # Technically the x is implicit here so integrate(f) would also work
plot(f, F, legend=True)
![../_images/efaa2720f390842ae44b28c2034e29cb627eddbabb462da036c1172f02cda930.png](../_images/efaa2720f390842ae44b28c2034e29cb627eddbabb462da036c1172f02cda930.png)
<sympy.plotting.backends.matplotlibbackend.matplotlib.MatplotlibBackend at 0x7fbc0c418400>
Notice, though, that SymPy does not add constant terms after it has integrated,
F
![../_images/a9ae242b4491c16c35e082f558af0a0664591934f605be7f768267546158acba.png](../_images/a9ae242b4491c16c35e082f558af0a0664591934f605be7f768267546158acba.png)
so a constant must be added afterwards manually if one in an exercise is asked to provide the indefinite integral,
C = symbols('C')
F = integrate(f, x) + C
F
![../_images/45847b67c1d76e113e67c92103a48fac1a1656decc28949f5506d2382751b5d3.png](../_images/45847b67c1d76e113e67c92103a48fac1a1656decc28949f5506d2382751b5d3.png)
To find the definite integral
one can make use of the anti-derivative as shown above and/or let SymPy do it with \(\verb|integrate(f, (var, a, b))|\). For example, if we for the Riemann integral are given:
# Simulated calculation by hand using the fundamental theorem of calculus
display(F.subs(x, 1) - F.subs(x, -1))
# SymPy's integrate function
integrate(f, (x, -1, 1))
![../_images/7f570e7a473d02cd6a6e85cef91377fa2886d900ce4e6ebe4899d7f55a2e0817.png](../_images/7f570e7a473d02cd6a6e85cef91377fa2886d900ce4e6ebe4899d7f55a2e0817.png)
![../_images/7f570e7a473d02cd6a6e85cef91377fa2886d900ce4e6ebe4899d7f55a2e0817.png](../_images/7f570e7a473d02cd6a6e85cef91377fa2886d900ce4e6ebe4899d7f55a2e0817.png)
Integration in Multiple Dimensions#
For functions of multiple variables the definite integral
can be determined using SymPy with
integrate(f, x1, x2)
and for definite integrals, the limits are added in tuples
just as for integration in one dimension:
integrate(f, (x1, a1, b1), (x2, a2, b2))
For example, have a look at the function
x1, x2 = symbols('x1 x2')
f2 = (x1*x2)**3 - x1**2 + sin(3 * x2) + x1*x2**2 + 4
f2
![../_images/0d2022cc42acfed7e5cf8f3cb8c33fd289d53d4a480c7f0b6d4de4500d953ead.png](../_images/0d2022cc42acfed7e5cf8f3cb8c33fd289d53d4a480c7f0b6d4de4500d953ead.png)
p = dtuplot.plot3d(f2, (x1, -1.2, 1.2), (x2, -1.2, 1.2), zlim=(-2, 6), use_cm=True, colorbar=False,
camera = {'azim': 25, 'elev': 20}, wireframe=True, show=False)
p.show()
![../_images/8e8702053420561abd2e62d05109ee747692cd7683546c4df99de655bbb22a00.png](../_images/8e8702053420561abd2e62d05109ee747692cd7683546c4df99de655bbb22a00.png)
An anti-derivative can be found with
integrate(f2, x1, x2) + C
![../_images/45f7e9968ef6a9e261864e7ccbae7535a8b791d5761bc1ffb366eee0aa479219.png](../_images/45f7e9968ef6a9e261864e7ccbae7535a8b791d5761bc1ffb366eee0aa479219.png)
and a definite volume under the graph can be found using the definite integral
integrate(f2, (x1, -1.2, 1.2), (x2, -1.2, 1.2))
![../_images/74c908dca15bb50c4c636807a3e240e3b46daa61dd09892f9053499b9fc92d93.png](../_images/74c908dca15bb50c4c636807a3e240e3b46daa61dd09892f9053499b9fc92d93.png)
Tips and Tricks for making Sympy Integrate#
When SymPy
runs into problems with an integral, here are some tricks that often will help.
Define Variables with Constraints#
State information about the variables you define. Is it for example a real variable? Then tell that to SymPy. The most frequent of such definitions are
real = True
positive = True
negative = True
nonnegative = True
nonzero = True
Some examples could be:
x1 = symbols('x1', real = True, nonnegative = True) # x1 is real and non-negative
x2 = symbols('x2', nonzero = True) # x2 is different from zero
Simplify before Integration#
Call simplify
on the integrand before integration is carried out. This can in many cases reduce the time required for SymPy
to integrate.
The work flow could look like this:
f = 1 / (x1 + x2)
integrand = simplify(f)
F = integrate(integrand, x1, x2)
F
![../_images/405a3751ff0a03f0e2f3a7b3c14a7f38f4644b8f72b1882b4ac32d3f22f29201.png](../_images/405a3751ff0a03f0e2f3a7b3c14a7f38f4644b8f72b1882b4ac32d3f22f29201.png)
Utilize the Linearity of Integrals#
Utilize the linearity of integrals. If you are to integrate the expression
then it will often be to your advantage to work with
\end{equation*} \alpha \int_{c}^d\int_{a}^{b} f_1(x_1,x_2) ,\mathrm dx,\mathrm dy + \beta \int_{c}^{d}\int_{a}^{b} f_2(x_1,x_2),\mathrm dx,\mathrm dy. \end{equation*}
An example could be
a = 3
b = 10
c = 2
d = 5
f = S(27)/3 * cos(x1) - 4 * atan(x2 -x1)
integrand1 = simplify(cos(x1)) # Here it is not necessary to use simplify
integrand2 = simplify(atan(x2 - x1))
F = S(27)/3 * integrate(integrand1, (x1, a,b), (x2, c,d)) - 4 * integrate(integrand2, (x1, a,b), (x2, c,d))
F
![../_images/0f0673d0f5624dfc1425aa6187478aa48dcec5416b62e910ee0705663e565fe3.png](../_images/0f0673d0f5624dfc1425aa6187478aa48dcec5416b62e910ee0705663e565fe3.png)
Find an Anti-Derivative first and Insert Limits next#
Sometimes SymPy
can find an anti-derivative but not evaluate the definite integral.
In such cases you can instead do the following:
f = tan(ln(x1**2))**3 / x1
integrand = simplify(tan(ln(x1**2))**3)
F = integrate(integrand, (x1,a,b))
F
![../_images/c26a8137e8df310e20a25b224cadfdb447b1267a0c36c560f4052094bfd0d1ae.png](../_images/c26a8137e8df310e20a25b224cadfdb447b1267a0c36c560f4052094bfd0d1ae.png)
We can instead find an anti-derivative and then insert the limits ourselves:
F = integrate(f, x1)
F.subs(x1, b) - F.subs(x1, a)
![../_images/2b0ba14bed22e5f384deb6efadb5bc25c86f55b1a6719292392f14e9d91c6c75.png](../_images/2b0ba14bed22e5f384deb6efadb5bc25c86f55b1a6719292392f14e9d91c6c75.png)
.doit()#
If SymPy
gives an output where the integration is not being carried out, one can in some cases force SymPy
to evaluate the integral using
This risks giving an error, though, and then one must use one of the other methods.
Numerical Methods#
You can in some cases risk that SymPy
either uses incredibly long time or cannot find a result for an integral at all, even after you have helped it out as much as you can. In such cases you might have no other choice than to solve the integral numerically. Methods for that require the numerical libraries numpy
and scipy
. You have already used those in Math1 as well as in your physics and programming courses, but if you have not installed them yet, then you can do that with:
or if one uses pip3
# quad is a numerical integrator for functions of one variable
# nquad is a numerical integrator for functions of multiple variables
from scipy.integrate import quad, nquad
import numpy as np
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[16], line 3
1 # quad is a numerical integrator for functions of one variable
2 # nquad is a numerical integrator for functions of multiple variables
----> 3 from scipy.integrate import quad, nquad
4 import numpy as np
ModuleNotFoundError: No module named 'scipy'
Below follows, respectively, an example of a numerical solution in one and in two dimensions.
One Dimension#
Let us have a look at the integral
f1 = tan(ln(x**2))**3 / x
# Note: SymPy won't determine the limited integral
integrate(f1, (x,pi,10))
This is one of those cases where SymPy
is able to find an analytical solution if we just help out a bit. But let us begin by solving the integral numerically.
Firstly we will convert \(f_1\) to a python/lambda function. This is done with the function lambdify
and corresponds, roughly speaking, to defining \(f_1\) by def f1(x):
f1_num = lambdify(x, f1, 'numpy')
# Now we can give f1_num a numerical value and get a numerical value back
f1_num(3)
Now, quad
can be used to determine the integral
F1_num, error1 = quad(f1_num, np.pi, 10)
F1_num, error1
The function quad
provides two outputs. The numerical approximation of the integral as well as an estimated error (deviation) from the exact solution.
We can try to compare this with the analytical solution. First we can force SymPy to evaluate the integral with \(\verb|.doit()|\), but this jsut results in an error:
# Remove # to see the error yourself. Maple
# integrate(f1, (x, pi, 10)).doit()
Instead we can use SymPy
to determine an anti-derivative, and then ourselves compute the definite integral:
# ad: anti-derivative
F1_ad = integrate(f1, x)
F1_analytic = trigsimp(F1_ad.subs(x,10) - F1_ad.subs(x,pi))
F1_ad, F1_analytic
Let us compare the two results:
F1_analytic.evalf() - F1_num
So, the numerical solution gives a fairly good approximation, but the two results are not entirely equal. Note that the true error that quad
incurs is smaller than the estimated error of \(3.4 \cdot 10^{−11}\). And that is not bad.
Two Dimensions#
Let us have a look at the expression
Here SymPy
has no problems finding a solution, so we can easily compare the analytical to the numerical method.
x,y = symbols('x y', real = True)
f2 = sin(x)*y**2
F2_analytic = integrate(f2, (y,-10,20), (x,0,pi))
F2_analytic
To evaluate the integral numerically we must again convert \(f_2\), using lambdify
.
f2_num = lambdify((x,y), f2, 'numpy') # lambdify can take multiple variables
# Now f2_num is a numerical function of two variables
f2_num(3,4)
Now the integral can be determined using nquad
. For a function of \(n\) variables, nquad
is called by:
F2_num, error2 = nquad(f2_num, [[0, np.pi], [-10,20]])
F2_num, error2
The two methods can now be compared:
F2_num - F2_analytic