How can I define an operator (eg integrator) in Python operating on a multi-variable function? My problem is that when I define an integrator function numint to numerically integrate a multivariable function along one of its variables that variable should be introduced at the first place, while I need it to be introduced by the user when the operator is called as it can change from one line of the code to another. One way is to not use an operator and implant the integration formula whenever it is needed with the variable over which integration should be computed, but that way the code would get quite cumbersome, so better to avoid it as much as is possible. The code that works wrongly for me is the following. Any idea to make it work properly?
var('x,y,z,t,x1,y1,z1,t1')
U=function('U0',x,y,z,t,x1,y1,z1,t1)
U=sin(x*y*z*t - x1*y1*z1*t1)^5
num=2
def numint(f, u, h):
integ = 0.5*h*(f(u=0) + f(u=1))
for i in range(1,num):
integ = integ + h * f(u=i*h)
return integ
print numint(U,t1,1/num)+numint(U,t,1/num)
Now the variable over which integration is taken is 'u', not found in thefunction 'U' at all, so that the result would be:
2*sin(t*x*y*z - t1*x1*y1*z1)^5
while I was expecting it to once integrate U WRT t and once WRT t1 then adding them together.
As much as I read the problem is that of the local and global variables. Maybe using closure and nested function is a suitable way to define an operator as is done in the example provided here, but the example is useful for single variable functions and I couldn't use it for multivariable functions.
Update. The following code does what I want (see the last print command, the previous prints are for experiment) but I should use a free argument of the function (like y) and tries with a dummy variable (like u in the example code above) fails:
var('x,y,z,t,x1,y1,z1,t1')
R=[x,y,z,t]
R1=[x1,y1,z1,t1]
U=function('U0',*(R+R1))
U0(x,y,z,t,x1,y1,z1,t1)=sin(x-x1)*sin(y-y1)*sin(z-z1)*sin(t-t1)
# if write U0=... instead of U0(...)=... the order of arguments of U0 is not
# specified and e.g. sin(x-x1) might become sin(x-y) from the system's viewpoint
num=2
def numint(func, h):
#integ = 0.5*h*(func(x=x,y=0,z=z,t=t,x1=x1,y1=y1,z1=z1,t1=t1) + func(x=x,y=1,z=z,t=t,x1=x1,y1=y1,z1=z1,t1=t1))
#for i in range(1,num):
# integ = integ + h * func(x=x,y=i*h,z=z,t=t,x1=x1,y1=y1,z1=z1,t1=t1)
#integ = 0.5*h*(func(x,0,z,t,x1,y1,z1,t1) + func(x,1,z,t,x1,y1,z1,t1))
#for i in range(1,num):
# integ = integ + h * func(x,i*h,z,t,x1,y1,z1,t1)
integ = 0.5*h*(func(y=0) + func(y=1))
for i in range(1,num):
integ = integ + h * func(y=i*h)
return integ
print numint(U,h=1/num),'\n'
print numint(U0,h=1/num),'\n\n'
print U0(y=z,z=y),'\n'
print numint(U(y=z,z=y),h=1/num),'\n'
print numint(U0(y=z,z=y),h=1/num),'\n\n'
print numint(U0(y=t,t=y),h=1/num).substitute(t=y)+numint(U0(y=t1,t1=y),h=1/num).substitute(t1=y)
the result is:
0.250000000000000*U0(x, 0, z, t, x1, y1, z1, t1) + 1/2*U0(x, 1/2, z, t,
x1, y1, z1, t1) + 0.250000000000000*U0(x, 1, z, t, x1, y1, z1, t1)
1/2*sin(-y1 + 1/2)*sin(z - z1)*sin(x - x1)*sin(t - t1) +
0.250000000000000*sin(-y1 + 1)*sin(z - z1)*sin(x - x1)*sin(t - t1) +
0.250000000000000*sin(z - z1)*sin(x - x1)*sin(t - t1)*sin(-y1)
sin(-y1 + z)*sin(y - z1)*sin(x - x1)*sin(t - t1)
0.250000000000000*U0(x, z, 0, t, x1, y1, z1, t1) + 1/2*U0(x, z, 1/2, t,
x1, y1, z1, t1) + 0.250000000000000*U0(x, z, 1, t, x1, y1, z1, t1)
1/2*sin(-z1 + 1/2)*sin(-y1 + z)*sin(x - x1)*sin(t - t1) +
0.250000000000000*sin(-z1 + 1)*sin(-y1 + z)*sin(x - x1)*sin(t - t1) +
0.250000000000000*sin(-y1 + z)*sin(x - x1)*sin(t - t1)*sin(-z1)
1/2*sin(-t1 + 1/2)*sin(z - z1)*sin(y - y1)*sin(x - x1) +
0.250000000000000*sin(-t1 + 1)*sin(z - z1)*sin(y - y1)*sin(x - x1) +
0.250000000000000*sin(t - 1)*sin(z - z1)*sin(y - y1)*sin(x - x1) +
1/2*sin(t - 1/2)*sin(z - z1)*sin(y - y1)*sin(x - x1) +
0.250000000000000*sin(z - z1)*sin(y - y1)*sin(x - x1)*sin(t) +
0.250000000000000*sin(z - z1)*sin(y - y1)*sin(x - x1)*sin(-t1)
Note how I was forced to use .substitute(). Here integration was uni-variable but when the dimension of integration increases this way of coding can become confusing. Any idea to do this cleaner and more directly?
var('x,y,z,t,x1,y1,z1,t1,var1')
U=function('U0',x,y,z,t,x1,y1,z1,t1)
U0(x,y,z,t,x1,y1,z1,t1)=sin(x-x1)*sin(y-y1)*sin(z-z1)*sin(t-t1)
# if write U0=... instead of U0(...)=... the order of arguments of U0 is not
# specified and e.g. sin(x-x1) might become sin(x-y) from the system's viewpoint
num=2.0
def numint(func, h):
#global var1
integ = 0.5*h*(func(var1=0) + func(var1=1))
for i in range(1,num):
integ = integ + h * func(var1=i*h)
return integ
print numint(U0(x=var1),1/num)+numint(U0(y1=var1),1/num)