Dear all,

I am trying to use Maple for Finite Element calculations. I have a 2d setup with linear basis functions and a 2d gaussian kernel that can rotate with respect to the axes. Attached please find the work sheet I am using.

Basis_function:

B := (x1,y1,x,y) -> max(0, 1-abs(x-x1))*max(0, 1-abs(y-y1))

transmissibility function:

t_hat:= (x1,y1,x,y) -> A*exp(-a*(x-x1)^2-2*b*(x-x1)*(y-y1)-c*(y-y1)^2)

where A and a,b,c are positive constants. a,b,c are calculated based on an angle phi and the two variances of the gaussian function.

I want to calculate the following function for different points (x1,y1) , (x2,y2):

trans := (x1, y1, x2, y2) -> int(int(B(x1, y1, xz, yz)*(int(int(t_hat(xz, yz, xp, yp)*(B(x2, y2, xz, yz)-B(x2, y2, xp, yp)), xp = x2-10*sigma1 .. x2+10*sigma1), yp = y2-10*sigma2 .. y2+10*sigma2)), xz = x1-hx .. x1+hx), yz = y1-hy .. y1+hy);

this integral in the form that is in the work sheet, works well for phi=0 and the results are what I want (numbers that go to zero as we move points 1 and 2 away from each other). for other values for phi it either gives an error (too many levels of recursion) or it returns expressions that seem unreasonable when I evaluate them (they don't go to zero).

for example, it doesn't work for phi = 0.5 at all. for phi = Pi/4 it will calculate some expression,

but as you move away from a point (e.g. trans(0,0,100,100)) the value does not become smaller than a certain value, but they should go to zero.

It seems that what I am trying to do is very sensitive to a,b,c, but actually it shouldn't be so different. I like to avoid exact integration, and just get a number, but I have no idea how to do this numerically. and I don't know how to write the problem in a way that would work for every angle phi.

any ideas?

thanks in advace,2d_maple_primes.mw

Download 2d_maple_primes.mw