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/CAHBzMo4GLc5iWTZO4muEucpOGtLaEFX3%2BtVjHtU3WZSL3KN_1A%40mail.gmail.com.
