Question: Numerical Integration, too many levels of recursion

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

with(plots); with(LinearAlgebra); with(ArrayTools)

``

Transmissibility function specifications

alpha := 1;

1

 

4

 

1

 

(1/4)*Pi

(1)

 

a := (1/2)*(cos(phi)/sigma1)^2+(1/2)*(sin(phi)/sigma2)^2;

17/64

 

15/64

 

17/64

 

A = (1/8)/Pi

(2)

 

 

Transmissibility*Kernel

t_hat := proc (x1, y1, x, y) options operator, arrow; A*exp(-a*(x-x1)^2-2*b*(x-x1)*(y-y1)-c*(y-y1)^2) end proc

proc (x1, y1, x, y) options operator, arrow; A*exp(-a*(x-x1)^2-2*b*(x-x1)*(y-y1)-c*(y-y1)^2) end proc

(3)

B := proc (a, b, x, y) options operator, arrow; max(0, 1-abs(x-a))*max(0, 1-abs(y-b)) end proc

proc (a, b, x, y) options operator, arrow; max(0, 1-abs(x-a))*max(0, 1-abs(y-b)) end proc

(4)

trans := proc (x1, y1, x2, y2) options operator, arrow; 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) end proc

proc (x1, y1, x2, y2) options operator, arrow; 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) end proc

(5)

 

#####testing here######

#for phi == 0 the results are what i want, numbers that go to zero as the points go far from each other. for phi != 0 trans returns an expression and the evaluation of that expression doesn't go to zero as we move the points far apart.

NULL

trans(0, 0, 0, 1)

trans(0, 0, 1, 0)

trans(0, 0, 5, 5)

``

This should be zero for any angle

trans(0, 0, 50, 50)

 

 

 

``


Download 2d_maple_primes.mw

Please Wait...