0

I wrote a code to compute the integrale of a function of one variable, using the midpoint rule.

import scipy
import numpy as np
from scipy import integrate
from numpy import linspace, sum

a=1
b=1e4
n=100

def f(x): 
    return 2*x

def midpoint(f,a,b,n):
     h=float(b-a)/n
     x=linspace(a+h/2, b-h/2,n)
     return h*sum(f(x))

It is working. Now i want to use the same approach for f(x,y), where I still just want to integrate over only the x variable, so a 1D integration.

precision: y is not a constant, for example y=linspace(1,10,10). I want to plot (y,integrale of f(x,y) over x).

Any idea ?

PS: thanks to @SUTerliakov a possible solution is

import scipy
import numpy as np
from scipy import integrate
from numpy import linspace, sum
from functools import partial

a=1 #minimum value taken by the variable x
b=1e4 #maximum value taken by the variable x
n=100 #number of points of integration 

def f(x,y): 
    return 2*x*y

def midpoint(g,a,b,n): #g(x) is an arbitrary function
    h=float(b-a)/n
    x=linspace(a+h/2, b-h/2,n)
    return h*sum(g(x))

y_arr=linspace(1,10,10)

for y in y_arr: 
    f_partial=partial(f,y) #the function f(x,y) for y=cte
    print(midpoint(f_partial,a,b,n))

The "for" loop is slow if my function is more complicated and if I want to integrate over a large scale. Is it possible to find a solution without using a "for" loop ?

Thank you.

Efalir
  • 1
  • 2
  • If you want to integrate only with respect to `x` then `y` becomes a constant, so, in that case, it doesn't change much, you can just add a `y` as a parameter of `midpoint` and give it a constant value. – Enric Grau-Luque May 18 '22 at 10:23
  • Hi @EnricGrau-Luque thank you for your answer, y is not a constant, i want to plot f(x,y) in function of y. I just added a comment about it. Sorry for that i did not precise it, thank you again :) – Efalir May 18 '22 at 11:06
  • You can generate your data for different `y` values with `[midpoint(partial(f, y=y), a, b, n) for y in y_arr]`, assuming `def f(x, y): return 2*x*y`, `y_arr = linspace(1,10,10)` and `from functools import partial`. `y` is constant for every integral - so you can fix it. No need to pass `y` to integration function, `partial` creates a function of `x` for given `y`. – STerliakov May 18 '22 at 12:05
  • Hi @SUTerliakov thank your for your help. I added the code to my message, it gave me the good result. However i would like to not use the for loop. Actually, here f(x,y) is just an example, my function is much more complicated and I want to increase the runtime. Have a nice day and thank you again :p – Efalir May 18 '22 at 12:46
  • Code you added is slightly wrong, `partial` should go to the last `for` loop (`y` is not defined in `midpoint` body). And if your function is not natively vectorizable, [for loop is good enough](https://stackoverflow.com/questions/35215161/most-efficient-way-to-map-function-over-numpy-array). – STerliakov May 18 '22 at 12:50
  • Oh, thank you. I made some changes, hope it is corrected now :) – Efalir May 18 '22 at 13:11

1 Answers1

0

PS: thanks to @SUTerliakov a possible solution is

import scipy
import numpy as np
from scipy import integrate
from numpy import linspace, sum
from functools import partial

a=1 #minimum value taken by the variable x
b=1e4 #maximum value taken by the variable x
n=100 #number of points of integration 

def f(x,y): 
    return 2*x*y

def midpoint(g,a,b,n): #g(x) is an arbitrary function
    h=float(b-a)/n
    x=linspace(a+h/2, b-h/2,n)
    return h*sum(g(x))

y_arr=linspace(1,10,10)

for y in y_arr: 
    f_partial=partial(f,y) #the function f(x,y) for y=cte
    print(midpoint(f_partial,a,b,n))
Efalir
  • 1
  • 2