I've made a simple program for numerically aproximating double integral, which accepts that the bounds of the inner integral are functions:
import numpy as np import time def double_integral(func, limits, res=1000): t = time.clock() t1 = time.clock() t2 = time.clock() s = 0 a, b = limits, limits outer_values = np.linspace(a, b, res) c_is_func = callable(limits) d_is_func = callable(limits) for y in outer_values: if c_is_func: c = limits(y) else: c = limits if d_is_func: d = limits(y) else: d = limits dA = ((b - a) / res) * ((d - c) / res) inner_values = np.linspace(c, d, res) for x in inner_values: t2 = time.clock() - t2 s += func(x, y) * dA t1 = time.clock() - t1 t = time.clock() - t return s, t, t1 / res, t2 / res**2
This is, however, terribly slow. When res=1000, such that the integral is a sum of a million parts, it takes about 5 seconds to run, but the answer is only correct to about the 3rd decimal in my experience. Is there any way to speed this up?
The code i am running to check the integral is
def f(x, y): if (4 - y**2 - x**2) < 0: return 0 #This is to avoid taking the root of negarive #'s return np.sqrt(4 - y**2 - x**2) def c(y): return np.sqrt(2 * y - y**2) def d(y): return np.sqrt(4 - y**2) # b d # S S f(x,y) dx dy # a c a, b, = 0, 2 print(double_integral(f, [a, b, c, d]))
The integral is eaqual to 16/9