so should i open a feature request? вт, 14 июн. 2022 г. в 20:27, Mom Mam <zfrhv2...@gmail.com>:
> 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 > <https://doc.sagemath.org/html/en/reference/calculus/sage/symbolic/integration/integral.html> > 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>: > >> 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) >> >> -- >> 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+unsubscr...@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 >> <https://groups.google.com/d/msgid/sage-devel/2fc49d91-6f52-4676-bb0f-4177a6badac4n%40googlegroups.com?utm_medium=email&utm_source=footer> >> . >> > -- You received this message because you are subscribed to the Google Groups "sage-devel" group. To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/sage-devel/CAHBzMo7wmM%3D40KsbyF_H7sGidFqhQFrbX1tZ6eD8398xRTjuOg%40mail.gmail.com.
