Problem of integration, the code run but the solution is not satisfying

Hi all! I coded the Eq.C.4 to find η(t) and the code works, but, as you can see from the fig. attached, my solution does not match with the one proposed by a scientific paper for the same equation and input data. Any suggestions?

Roi=916;
Row=1024; 
R=Row/Roi;

r=50; 
d=50; 

%column geometry 
H0=50; 
l=5; 
b=5;   

a=(2*l*H0/pi)^0.5; 


Area=pi*a^2; 
Vol=Area*H0; 

%Dinamic of the block
Omega=1; 

va=(3*(pi)^(3/2))/(16*sqrt(2))*(Omega*sqrt(l*H0))/(1+pi/4*R*(l/b)); 

%Surfaceelevation
E=zeros(1,301);
t=0:0.1:30;


for i=1:301
O=@(x,y)(besselj(0,x.*y.*a/d)).*((1-y.^2).^0.5).*y.*besselj(0,x*r/d).*((x.*tanh(x)).^0.5).*sin((((9.81/d*t(i).^2).*x).*tanh(x)).^0.5).*x;
E(i)=integral2(O,0,2000,0,1,'Method','iterated','AbsTol',1e-10,'RelTol',1e-10)
end


E1=-2*(a.^3)/(d^3.5)/(9.81^0.5)*va*d*E;

[Are Mjaavatten]

"Split and conquer" is always a good idea for a complex problem like this. When trying to do everything at once, the chances of getting it right are rather small. Split into smaller parts and test to make sure they make sense before moving on. Here is a possible start:

Instead of using a double integral, try for instance to first calculate the B function:

function B = attiliB(alpha,x)
B = zeros(size(x));
for i = 1:length(x)
  fun = @(y) sqrt(1-y.^2).*besselj(0,alpha*x(i).*y).*y;
  B(i) = integral(fun,0,1);
end

Test it by:

>> b = 25.2;
>> h = 50;
>> x = linspace(0,1);
>> z = attiliB(a,x);
>> plot(x,z)

Now you can use attiliB inside a single integral for the I_ksi function. Doing this instead of a double integral may perhaps be slower, but the chances of getting it right are certainly better.

It turns out that B(x) in this case is very well approximated by a second order polynomial, which you can use as a slightly less accurate, but certainly much faster, substitute for attiliB:

>> pB = polyfit(x,z,2);
>> plot(x,z,x,polyval(pB,x),'.')

Good luck with the rest!

[Tommaso Attili]

Thanks for that, I will try it!! Just a question: why x is defined between 0 and 1?

[Are Mjaavatten]

It seems that I mixed up the range of x and of y. I also forgot a y in line 4 of attiliB.m. (I have now corrected that in the previous comment.)

The good news is that there exists an analytical expression for the integral:

a = b*x/(2*h);
B = (sin(a)-a.*cos(a))./a.^3;

This will hopefully speed up the calculations significantly.

[Tommaso Attili]

Many Thanks! But the solution is still not good, any suggestions for the other integral?

[Are Mjaavatten]

I am not able to reproduce the figure from the paper either. Your curve at least resembles the original. Mine is totally different. Here is my attempt:

h = 50;
b = 25.2;
Va = 6.21;
g = 9.81;
r = 50;
a = b/2/h;
int1 = @(x) besselj(0,r/h*x).*(sin(a*x)-a*x.*cos(a*x))/a^3./x.^2;
int2 = @(x,t) sqrt(x.*tanh(x)).*sin(g*t^2/h*x.*tanh(x));
integrand = @(x,t) int1(x).*int2(x,t);
N = 200;
t = linspace(0,5,N);
I = zeros(1,N);
for i = 1:N
  fun = @(x) integrand(x,t(i));
  I(i) = integral(fun,0,200);
end
eta = -2*a^3*Va/sqrt(g*h)*I;
figure; plot(t,eta)

Increasing the upper integration limit dampens out the wiggles, but the general shape is unchanged. There seems to be two possible explanations:

1) I have made some mistake (misunderstanding or bug).

2) The figure in the original paper is based on different parameters and/or equations.

I leave it to you to try to figure this out.

[Tommaso Attili]

Thanks! Yes I plotted your solution. At this point I really don't have idea. Thanks for that.