CODEWITHSUNDEEP
pythonnumerical-methodsnumerical-integration
JD_PM
JD_PM

Reputation: 159

How to use numerical methods to approximate the solution of an integral

I have the following integrals (more details: https://math.stackexchange.com/questions/3193669/how-to-evaluate-the-line-integral-checking-stokes-theorem)

enter image description here

C_3 can be evaluated with trigonometric tricks. Then you can solve it by:

import sympy as sp
t = sp.symbols('t')
sp.Integral((1-sp.cos(t)-sp.sin(t))**2 * sp.exp(1-sp.cos(t)-sp.sin(t)) * (sp.sin(t)-sp.cos(t)), (t, 0, 2*sp.pi))

The problem are C_1 and C_2. These cannot be evaluated with tricks. Then, I have to use numerical methods.

What do you suggest? I have been trying with N() but got nothing.

Thanks.

Upvotes: 1

Views: 191

Answers (2)

Nico Schlömer
Nico Schlömer

Reputation: 58721

Alternative: Use quadpy (a project of mine):

import quadpy
from numpy import cos, sin, exp, pi

c1, err1 = quadpy.quad(
    lambda t: (1 + sin(t)) * exp(1 + cos(t)) * (-sin(t)), 0.0, 2 * pi,
)

c2, err2 = quadpy.quad(
    lambda t: ((1 + cos(t)) ** 2 + exp(1 + cos(t))) * cos(t), 0.0, 2 * pi,
)

print("C1 = ", c1, ", estimated error: ", err1)
print("C2 = ", c2, ", estimated error: ", err2)

C1 =  -9.652617076333142 , estimated error:  1.3725463615061705e-09
C2 =  15.9358023895608 , estimated error:  6.646678031309946e-11

Upvotes: 1

s.ouchene
s.ouchene

Reputation: 1869

You can use scipy.integrate.quad function:

from scipy.integrate import quad
from numpy import cos, sin, exp, pi

f1 = lambda t: (1 + sin(t))*exp(1+cos(t))*(-sin(t))
f2 = lambda t: ((1 + cos(t))**2 + exp(1+cos(t)))*cos(t)

C1, err1 = quad(f1, 0, 2*pi)
C2, err2 = quad(f2, 0, 2*pi)

print("C1 = ", C1, ", estimated error: ", err1)
print("C2 = ", C2, ", estimated error: ", err2)

Output:

C1 =  -9.652617083240306, estimated error:  2.549444932020608e-09
C2 =  15.93580239041989, estimated error:  3.4140955340600243e-10

EDIT: You can also specify the precision via the arguments: epsrel: relative error, epsabs: absolute error. But this is a little bit tricky (See this): we specify an absolute error target of zero. This condition cannot be satisfied, and so the relative error target will determine when the integration stops.

C1, err1 = quad(f1, 0, 2*pi, epsrel=1e-10, epsabs=0)
print("C1 = ", C1, ", estimated error: ", err1)

Output: 

C1 =  -9.652617083240308 , estimated error:  1.4186554373311127e-13

Upvotes: 1

Related Questions