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...