[feature] simplify manual integral usage

[Mom Mam]

hello, i saw that i should post my feature request here, and only then open a ticket.

i think it would be nice to simplify the usage of "numerical_integral":

1.

instead of using labmda, pass the integration variable to the function just as in the "integral" function

instead of this:

numerical_integral(lambda x: x*y, 0, 100)

do this:

numerical_integral(x*y, x, 0, 100)

2.

make it possible to use variables inside the functio

instead of this:

x = SR.var('x')

func = x^2

numerical_integral(lambda  new_x: func.subs({x: new_x}),  0, 100)

do this:

x = SR.var('x')

func = x^2

numerical_integral(lambda  x: func,  0, 100)

# or this

numerical_integral(func, x,  0, 100)

3.

support integrals for vectors

instead of this:

x = SR.var('p')

vec = vector([p, p^2, 0]))

new_x = numerical_integral(vec[0], p,  0, 100)

new_y = numerical_integral(vec[1], p,  0, 100)

new_z = numerical_integral(vec[2], p,  0, 100)

result = vector([new_x, new_y, new_z]))

do this:

x = SR.var('p')

vec = vector([p, p^2, 0]))

result numerical_integral(vec, p,  0, 100)

4.

support double integrals

instead of this:

numerical_integral(numerical_integral(x*y^2, 0, 100)[0], 0, 100)

do this:

numerical_integral(x*y^2, [x,y], [0,0], [100,100])

end goal:

so there wont be messes like this:

    x = numerical_integral(lambda new_beta: numerical_integral(lambda new_alpha: func[0].subs({alpha: new_alpha, beta: new_beta}), 0, 100)[0], 0, 100)

    y = numerical_integral(lambda new_beta: numerical_integral(lambda new_alpha: func[1].subs({alpha: new_alpha, beta: new_beta}), 0, 100)[0], 0, 100)

    z = numerical_integral(lambda new_beta: numerical_integral(lambda new_alpha: func[2].subs({alpha: new_alpha, beta: new_beta}), 0, 100)[0], 0, 100)

and instead it would be simple as this:

    vec = numerical_integral(func, [alpha, beta], [0,0], [100,100])

sorry for the mess, i wish i could use markdown

[Nils Bruin]

On Monday, 13 June 2022 at 13:38:55 UTC-7 zfrh...@gmail.com wrote:

hello, i saw that i should post my feature request here, and only then open a ticket.

i think it would be nice to simplify the usage of "numerical_integral":

1.

instead of using labmda, pass the integration variable to the function just as in the "integral" function

instead of this:

numerical_integral(lambda x: x*y, 0, 100)

do this:

numerical_integral(x*y, x, 0, 100)

You don't have to use lambda. As long as the specified integrand is a (one-argument) callable, you're OK. However, because it's numerical integration, the callable must evaluate to an actual number. Variables cannot be in the output. So, unless "y" is bound to a specific value, the above will not work.

You can use

numerical_integral(x.function(x), 0, 100)

if you prefer.

2.

make it possible to use variables inside the functio

instead of this:

x = SR.var('x')

func = x^2

numerical_integral(lambda  new_x: func.subs({x: new_x}),  0, 100)

do this:

x = SR.var('x')

func = x^2

numerical_integral(lambda  x: func,  0, 100)

# or this

numerical_integral(func, x,  0, 100)

it works if you set func = (x^2).function(x).

You have to be specific about which variable your expression should be considered a function of. Since symbolic arguments have no special meaning, growing a different call signature (f,x,0,100) seems a little misplaced.

You can write things a little more glibly if you define:

Cx=CallableSymbolicExpressionRing([x])

then you can use

numerical_integral(Cx(sin(x)),0,100)

3.

support integrals for vectors

instead of this:

x = SR.var('p')

vec = vector([p, p^2, 0]))

new_x = numerical_integral(vec[0], p,  0, 100)

new_y = numerical_integral(vec[1], p,  0, 100)

new_z = numerical_integral(vec[2], p,  0, 100)

result = vector([new_x, new_y, new_z]))

do this:

x = SR.var('p')

vec = vector([p, p^2, 0]))

result numerical_integral(vec, p,  0, 100)

Check out:

import sage.numerical.gauss_legendre

 

4.

support double integrals

instead of this:

numerical_integral(numerical_integral(x*y^2, 0, 100)[0], 0, 100)

do this:

numerical_integral(x*y^2, [x,y], [0,0], [100,100])

The numerical_integral routine you're referring to is mainly just wrapping scipy. It has multiple numerical integration code as well:

https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.dblquad.html#scipy.integrate.dblquad

 this might fit your goals better.

Note that an interface for a numerical double integrator should accept integration domains that are a little more general than just rectangular regions. At least (like dblquad above), it should allow for iterated integrals, where the bounds of the inner integral are allowed to be functions of the outer variable.

end goal:

so there wont be messes like this:

    x = numerical_integral(lambda new_beta: numerical_integral(lambda new_alpha: func[0].subs({alpha: new_alpha, beta: new_beta}), 0, 100)[0], 0, 100)

    y = numerical_integral(lambda new_beta: numerical_integral(lambda new_alpha: func[1].subs({alpha: new_alpha, beta: new_beta}), 0, 100)[0], 0, 100)

    z = numerical_integral(lambda new_beta: numerical_integral(lambda new_alpha: func[2].subs({alpha: new_alpha, beta: new_beta}), 0, 100)[0], 0, 100)

You can clean up that code yourself already. To compute

int_a^b i(nt_c^d f(x,y) dy ) dx

def double_int(f,a,b,c,d):

    def g(x):

        return numerical_integral( lambda y: f(x,y), c,d)

    return numerical_integral(g,a,b)

Just make sure that you pass in a value for f that is a callable in two arguments, e.g.

f = func[0].function(beta,alpha)

x= double_int(f,0,100,0,100)

(but mind you: using dblquad may be significantly better)

[Mom Mam]

Check out:
import sage.numerical.gauss_legendre

looks interesting but i cant find where to put the limits

---

The numerical_integral routine you’re referring to is mainly just wrapping scipy. It has multiple numerical integration code as well:
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.dblquad.html#scipy.integrate.dblquad

it doesnt seems to support vectors, maybe it can be useful but i want very easy usage and generic,
i dont want separate function for each object of specific usecase.
like the integral function, it looks so great and easy, works with vectors, double integrals, you can specify the integral parameters, limits… basically perfect
the only problem with that integral is that i wait for hours and nothing happens, i think the integral is just too difficult, so i use the numerical integral instead.

---

You can clean up that code yourself already. To compute> int_a^b i(nt_c^d f(x,y) dy ) dx

its still difficult to understand :/

---

i have the next code:

## sagemath ##

def circles_force(circle_1, circle_2, constant, manual):
  alpha, beta = SR.var('alpha, beta')
  assume(0 < alpha, alpha < 2*pi, 0 < beta, beta < 2*pi)

  # finding dot
  if round(circle_1[1][0]/circle_1[1].norm(), 4) == 0:
    new_x_1 = circle_1[1].cross_product(vector([1,0,0])).normalized()
  else:
    new_x_1 = circle_1[1].cross_product(vector([0,1,0])).normalized()

  new_y_1 = circle_1[1].cross_product(new_x_1)
  new_x_1 = new_x_1*circle_1[1].norm()

  if round(circle_2[1][0]/circle_2[1].norm(), 4) == 0:
    new_x_2 = circle_2[1].cross_product(vector([1,0,0])).normalized()
  else:
    new_x_2 = circle_2[1].cross_product(vector([0,1,0])).normalized()

  new_y_2 = circle_2[1].cross_product(new_x_2)
  new_x_2 = new_x_2*circle_2[1].norm()

  relative_place_1 = sin(alpha)*new_x_1 + cos(alpha)*new_y_1
  relative_place_2 = sin(beta)*new_x_2 + cos(beta)*new_y_2

  v_1 = circle_1[1].cross_product(relative_place_1).normalized()
  v_2 = circle_2[1].cross_product(relative_place_2).normalized()

  absolute_place_1 = circle_1[0] + relative_place_1
  absolute_place_2 = circle_2[0] + relative_place_2

  R = absolute_place_2 - absolute_place_1

  # "their" force calculation
  f_1 = constant * v_2.cross_product(-R.normalized()).cross_product(v_1) / R.norm()^2
  f_2 = constant * v_1.cross_product(R.normalized()).cross_product(v_2) / R.norm()^2

  # rotating force
  f_1_rotating = f_1.cross_product(relative_place_1)
  f_2_rotating = f_2.cross_product(relative_place_2)

  # integrals

  def manual_integral(func, i, f):
    ## works very bad
    # x = monte_carlo_integral(lambda new_alpha,new_beta: func[0].subs({alpha: new_alpha, beta: new_beta}), i, f, accuracy)
    # y = monte_carlo_integral(lambda new_alpha,new_beta: func[1].subs({alpha: new_alpha, beta: new_beta}), i, f, accuracy)
    # z = monte_carlo_integral(lambda new_alpha,new_beta: func[2].subs({alpha: new_alpha, beta: new_beta}), i, f, accuracy)

    x = numerical_integral(lambda new_beta: numerical_integral(lambda new_alpha: func[0].subs({alpha: new_alpha, beta: new_beta}), i[0], f[0])[0], i[1], f[1])
    y = numerical_integral(lambda new_beta: numerical_integral(lambda new_alpha: func[1].subs({alpha: new_alpha, beta: new_beta}), i[0], f[0])[0], i[1], f[1])
    z = numerical_integral(lambda new_beta: numerical_integral(lambda new_alpha: func[2].subs({alpha: new_alpha, beta: new_beta}), i[0], f[0])[0], i[1], f[1])

    return [ vector([x[0],y[0],z[0]]), vector([x[1],y[1],z[1]]) ]

  # calculating the integral
  if manual:
    F_1_T, F_1_E = manual_integral(f_1, [0,0], [2*pi,2*pi])
    F_2_T, F_2_E = manual_integral(f_2, [0,0], [2*pi,2*pi])

    F_1_rotating_T, F_1_rotating_E = manual_integral(f_1_rotating, [0,0], [2*pi,2*pi])
    F_2_rotating_T, F_2_rotating_E = manual_integral(f_2_rotating, [0,0], [2*pi,2*pi])

  else:
    F_1_T = f_1.integral(alpha, 0, 2*pi).integral(beta, 0, 2*pi)
    F_1_E = vector([0, 0, 0])

    F_2_T = f_2.integral(alpha, 0, 2*pi).integral(beta, 0, 2*pi)
    F_2_E = vector([0, 0, 0])

    F_1_rotating_T = f_1_rotating.integral(alpha, 0, 2*pi).integral(beta, 0, 2*pi)
    F_1_rotating_E = vector([0, 0, 0])

    F_2_rotating_T = f_2_rotating.integral(alpha, 0, 2*pi).integral(beta, 0, 2*pi)
    F_2_rotating_E = vector([0, 0, 0])

  return [F_1_T, F_2_T, F_1_rotating_T, F_2_rotating_T, F_1_E, F_2_E, F_1_rotating_E, F_2_rotating_E]

## configurations

manual = true

# constants
constant = 1

# set points (center, direction)
circle_1 = (vector([0,0,0]),vector([1,0,0]))
circle_2 = (vector([5,0,0]),vector([0,1,0]))

res = circles_force(circle_1, circle_2, constant, manual)
F_1_T, F_2_T, F_1_rotating_T, F_2_rotating_T, F_1_E, F_2_E, F_1_rotating_E, F_2_rotating_E = map(lambda vec: vector([round(vec[0],3),round(vec[1],3),round(vec[2],3)]), res)

print(f"{F_1_T          = }")
print(f"{F_2_T          = }")
print(f"{F_1_rotating_T = }")
print(f"{F_2_rotating_T = }")
#print(f"{F_2_E          = }")
#print(f"{F_2_rotating_E = }")

#ratio = F_2_T.norm() / F_2_rotating_T.norm()
#print(f"{ratio          = }")

so all i want is this:

  def manual_integral(func, i, f):
    x = numerical_integral(lambda new_beta: numerical_integral(lambda new_alpha: func[0].subs({alpha: new_alpha, beta: new_beta}), i[0], f[0])[0], i[1], f[1])
    y = numerical_integral(lambda new_beta: numerical_integral(lambda new_alpha: func[1].subs({alpha: new_alpha, beta: new_beta}), i[0], f[0])[0], i[1], f[1])
    z = numerical_integral(lambda new_beta: numerical_integral(lambda new_alpha: func[2].subs({alpha: new_alpha, beta: new_beta}), i[0], f[0])[0], i[1], f[1])

    return [ vector([x[0],y[0],z[0]]), vector([x[1],y[1],z[1]]) ]

  F_1_T, F_1_E = manual_integral(f_1, [0,0], [2*pi,2*pi])
  F_2_T, F_2_E = manual_integral(f_2, [0,0], [2*pi,2*pi])

  F_1_rotating_T, F_1_rotating_E = manual_integral(f_1_rotating, [0,0], [2*pi,2*pi])
  F_2_rotating_T, F_2_rotating_E = manual_integral(f_2_rotating, [0,0], [2*pi,2*pi])

look a lot more simple like this:

    F_1_T, F_1_E = numerical_integral(f_1, [alphta, beta], [0,0], [2*pi,2*pi])
    F_2_T, F_2_E = numerical_integral(f_2, [alphta, beta], [0,0], [2*pi,2*pi])

    F_1_rotating_T, F_1_rotating_E = numerical_integral(f_1_rotating, [alphta, beta], [0,0], [2*pi,2*pi])
    F_2_rotating_T, F_2_rotating_E = numerical_integral(f_2_rotating, [alphta, beta], [0,0], [2*pi,2*pi])

easy to use, easy to understand

вт, 14 июн. 2022 г. в 04:02, Nils Bruin <nbr...@sfu.ca>:

--
You received this message because you are subscribed to a topic in the Google Groups "sage-devel" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sage-devel/baK6jz9z1og/unsubscribe.
To unsubscribe from this group and all its topics, send an email to sage-devel+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sage-devel/2fc49d91-6f52-4676-bb0f-4177a6badac4n%40googlegroups.com.

[Mom Mam]

so should i open a feature request?

вт, 14 июн. 2022 г. в 20:27, Mom Mam <zfrh...@gmail.com>:

[Nils Bruin]

On Tuesday, 14 June 2022 at 10:39:05 UTC-7 zfrh...@gmail.com wrote:

Check out:
import sage.numerical.gauss_legendre

looks interesting but i cant find where to put the limits

Raw implementations of numerical integration methods usually standardize to an integration interval [-1,1] or [0,1]; see documentation. In those cases, transforming the integrand to that path is left to the user. Interfacing code can provide a more pleasing interface.

 

it doesnt seems to support vectors, maybe it can be useful but i want very easy usage and generic,

i dont want separate function for each object of specific usecase.
like the integral function, it looks so great and easy, works with vectors, double integrals, you can specify the integral parameters, limits… basically perfect
the only problem with that integral is that i wait for hours and nothing happens, i think the integral is just too difficult, so i use the numerical integral instead.

Yes, there is a reason for that: to get good performance you need assumptions on the integrand. In particular the adaptive methods that numerical_integral uses to get decent heuristical results for many integrands require thought to generalize to vectors. In general it's probably better to use the adaptive techniques on each component separately. The gauss-legendre implementation above was motivated in cases where the vector of values can be evaluated at relatively little additional cost to one individual component (because of a lot of shared expressions between the components)

 Interface design for general, performant numerical methods is hard. Generally, the key is to have the fundamental functionality available (we basically have that), then see what people need and then wait until someone knowledgeable and willing to put the time in to program something of a decent level of generality that is convenient.

In the mean time, the best course of action is to program your own convenience wrappers to make the raw functionality accessible to your own purposes. Perhaps the experiences of that lead you to program something that can convince others that it's of sufficient generality and convenience to include in sagemath. The raw functionality is the most important part, though. Writing convenience wrappers is part of using expert software.