CODEWITHSUNDEEP
pythonpython-3.xnumpycontinuous-integrationscipy
user3440244
user3440244

Reputation: 371

solve an equation containing integrals

I need to solve an integral equation embedded with another integral equation by python 3.2 in win7.

There are 2 integral equations.

The code is here: 

import numpy as np
from scipy.optimize.minpack import fsolve
from numpy import exp
from scipy.integrate.quadpack import quad
import matplotlib.pyplot as plt
import sympy as syp

lb = 0 
def integrand2(x, a):
    print("integrand2 called")
    return x**(a-1) * exp(-x)

def integrand1(x, b, c):
    print("integrand1 called")
    integral , err = quad(integrand2, lb/b, syp.oo , args=(1+c))
    return c/(b*integral)

def intergralFunc1(b,c):
    integral,err = quad(integrand1, 0, 10, args=(b,c)) 
    print("integral is ", integral, " and err is ", err)
    print("b is ", b, " and c is ", c)
    return 10 - integral

def findGuess():
    vfunc = np.vectorize(intergralFunc1)
    f1 = np.linspace(0.01, 10,10)
    f2 = np.linspace(0.01, 10,10)
    result = vfunc(f1, f2)
    plt.plot(f1, result)
    plt.xlabel('f1')
    plt.subplot(211)

    plt.plot(f2, result)
    plt.xlabel('f2')
    plt.subplot(212)
    plt.show()

def solveFunction():
    sol= fsolve(intergralFunc1, 5, 5, full_output=True)
    return sol

if __name__ == '__main__':
     findGuess()

    sol = solveFunction()
    print("sol is ", sol)
    print("verification: \n")
    print("f(b,c) is ", intergralFunc1(sol[0],5))

I got the results that make no sense. 

     integral is  nan  and err is  nan
     b is  [ 5.]  and c is  5
     f(b,c) is  nan

Any help would be appreciated !!!

Upvotes: 0

Views: 780

Answers (2)

pv.
pv.

Reputation: 35145

One problem is in using sympy.oo for infinity.

In [36]: quad(integrand2, 0, syp.oo, args=(1+5,))
(nan, nan)

In [37]: quad(integrand2, 0, np.inf, args=(1+5,))
(120.0, 2.2281810652e-07)

In [38]: quad(integrand2, 0, float(syp.oo), args=(1+5,))
(120.0, 2.2281810652e-07)

It's not a floating-point number:

In [39]: type(syp.oo)
sympy.core.numbers.Infinity

Not sure exactly why it doesn't work as expected, however.

Upvotes: 0

hivert
hivert

Reputation: 10667

For some reason integrand1 is returning a numpy.ndarray where it is expected to return a float.

Also on my machine numpy reported

 y:289: UserWarning: Extremely bad integrand behavior occurs at some points of the
  integration interval.

which means that your problem is numerically very instable. Therefore nonsensical results are to be expected.

Upvotes: 1

Related Questions