Python Forum
How to run this triple integration in python
Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
How to run this triple integration in python
#1
When I run the following triple integral in python, it gives me an error:

 IntegrationWarning: The occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved
Here is the code, it is to plot a triple integral f(t1, t2, m, k) against a variable k.


import numpy as np
import scipy.special
from scipy import integrate
from scipy.special import kn
import matplotlib.pyplot as plt

import math
import time, sys

H = 4
ti = -H
a = 0.1

f = lambda t1, t2, m, k: (k**3 * m**3) * (1/(H * t1)) * (1/(H * t2)) * (2/k)**2

X = np.arange(0, 50, 0.1)
g = float('inf')

def F(x):
    res = np.zeros_like(x)
    for i, val in enumerate(x):
        y, err = integrate.tplquad(f, ti, a, lambda x: ti, lambda x: a,  lambda x, y: 0, lambda x, y: g, args=(val,))
        res[i] = y
    return res

plt.plot(X, F(X))
plt.title("P(k), H=2.5")
plt.xlabel("k")
plt.ylabel("P")
plt.savefig('P(k).png')
plt.show()
How can I run the code without IntegrationWarning?
Reply
#2
There are more errors

Output:
/tmp/tmp.py:16: RuntimeWarning: divide by zero encountered in scalar divide f = lambda t1, t2, m, k: (k**3 * m**3) * (1/(H * t1)) * (1/(H * t2)) * (2/k)**2 /tmp/tmp.py:16: RuntimeWarning: invalid value encountered in scalar multiply f = lambda t1, t2, m, k: (k**3 * m**3) * (1/(H * t1)) * (1/(H * t2)) * (2/k)**2 /usr/lib/python3/dist-packages/scipy/integrate/_quadpack_py.py:1225: IntegrationWarning: The occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. The error may be underestimated. quad_r = quad(f, low, high, args=args, full_output=self.full_output,
Reply


Possibly Related Threads…
Thread Author Replies Views Last Post
  SyntaxError: EOF while scanning triple-quoted string literal Caroline 5 8,703 Feb-07-2018, 07:47 AM
Last Post: micseydel
  Integration of Python and C using .dll files zbiegler 16 15,006 Mar-17-2017, 01:40 PM
Last Post: sparkz_alot

Forum Jump:

User Panel Messages

Announcements
Announcement #1 8/1/2020
Announcement #2 8/2/2020
Announcement #3 8/6/2020