CODEWITHSUNDEEP
pythonintegrationnumerical-methods
dustin
dustin

Reputation: 4416

Python: double integration on an infinite domain

I am trying to integrate my function over u and xx and then store the values in a matrix so I can plot them with imshow or pcolormesh. The bounds on the integration are 0 < u < inf and -inf < xx < inf. At the moment, I am only taking the bounds to be 10 until I can figure this out.

import numpy as np
import pylab as pl
from scipy.integrate import dblquad

b = 50.0

x = np.linspace(-10, 10, 1000)
y = np.linspace(0, 10, 1000)

T = pl.zeros([len(x), len(y)])


for xi in enumerate(x):
    for yi in enumerate(y):
        def f(xi, yi, u, xx):
            return ((np.exp(u * (b - yi)) - np.exp(-u * (b - yi))) /
                    (np.exp(u * b) - np.exp(-u * b)) * np.cos(u * (xx - xi)))


def fint(u, xx):
    return T + dblquad(f, -10, 10, 0.1, 10, args = (u, xx))[0]

This is the code I have so far but I know it isn't working properly; unfortunately, I don't know what the problem is. Maybe I can't have the two for loops in my definition of f or my my fint is wrong.

Upvotes: 0

Views: 962

Answers (2)

hunse
hunse

Reputation: 3255

It's not completely clear from your question what you're trying to do. But here's how I interpreted it: You have a double integral over two variables u and xx, which also takes two parameters xi and yi. You want to evaluate the integral over xx and u at many different values of xi and yj, and store these values in T. Assuming this is what you want to do (and correct me if I'm wrong), here's how I would do it.

import numpy as np
from scipy.integrate import dblquad

b = 50.0

x = np.linspace(-10, 10, 1000)
y = np.linspace(0, 10, 1000)

def f(xx, u, xi, yj):
    return ((np.exp(u * (b - yj)) - np.exp(-u * (b - yj))) /
            (np.exp(u * b) - np.exp(-u * b)) * np.cos(u * (xx - xi)))

T = np.zeros([len(x), len(y)])
for i, xi in enumerate(x):
    for j, yj in enumerate(y):
        T[i, j] += dblquad(
            f, -10, 10, lambda x: 0.1, lambda x: 10, args=(xi, yj))[0]

Upvotes: 1

Jblasco
Jblasco

Reputation: 3967

fint is the only thing that calls f. You are not calling fint, which means that f is not being used at all, just defined about a million times. I would consider defining the function just once and calling it a million times instead.

Upvotes: 0

Related Questions