Why the integration can not be calculated out correctly?

As can be seen in the following codes, jr is the integration which need to be calculated. TR is the integrand, and x and y are the variables. Here, we want to calculate a curve of jr versus v, as can be seen in the main program.

The first problem is that, when v is from 0.46 to 1, we found that the curve is not smooth, and the curve of jr and its differential curve are shown in the figures below.

Since the integrand TR is continuously derivative, the curve of jr versus v should be smooth.

The second problem is that, when v is from 0.01 to 0.45, the integration jr can not be worked out or the value is obviously wrong. The preliminary reason may be that the results of Airy function in TR are calculated as 0 or infinity, but in fact they are not. How to solve the above two problems? Many thanks!

main program:

clc
clear
%fid=fopen('D:\2v=0.5-1.txt','wt');
iv = 0;
for v = 0.5: 0.01: 1
    iv = iv + 1;
    jr(iv) = integral(@(y) integral(@(x) TR(x,y,v), 9.578251022212591 - v, 9.658875554160838), 0, 9.658875554160838,'ArrayValued',1);
    %fprintf(fid,'%g %g\n',v,jr);
end
%disp(jr);
plot(0.5:0.01:1,jr)

subprogram:

function T=TR(x,y,v)
k1 = (51.2317 * (x + 9.57825102221259));
k2 = (51.2317 * (x + v - 9.57825102221259));
rho1 = 23.5875 * (v - 0.4120716294548327)^(1/3);
lambda01 = 4.7175048762838045 * (11.630653268960174 - x) / (v - 0.4120716294548327)^(2/3);
lambdad1 = 4.7175048762838045 * (12.042724898415006 - v - x) / (v - 0.4120716294548327)^(2/3);
Ai_lambda01 = airy(0, lambda01);
AiP_lambda01 = airy(1, lambda01);
Bi_lambda01 = airy(2, lambda01);
BiP_lambda01 = airy(3, lambda01);
Ai_lambdad1 = airy(0, lambdad1);
AiP_lambdad1 = airy(1, lambdad1);
Bi_lambdad1 = airy(2, lambdad1);
BiP_lambdad1 = airy(3, lambdad1);
term1 = rho1 .* AiP_lambda01 .* BiP_lambdad1 ./ k1 + k2 .* Ai_lambda01 .* Bi_lambdad1 ./ rho1 - rho1 .* AiP_lambdad1 .* BiP_lambda01 ./ k1 - k2 .* Ai_lambdad1 .* Bi_lambda01 ./ rho1;
term2 = -Ai_lambda01 .* BiP_lambdad1 + k2 .* AiP_lambda01 .* Bi_lambdad1 ./ k1 - k2 .* Ai_lambdad1 .* BiP_lambda01 ./ k1 + AiP_lambdad1 .* Bi_lambda01;
c1 = 4 .* pi^(-2) ./ (term1.^2 + term2.^2);
traverse = @(x, y) (x + y > 9.658875554160838 - v) & (x + y < 9.658875554160838);
T=double((k2.*abs(c1).^2./k1).*traverse(x,y));
end

[Torsten]

Is your 2d region of integration a triangle ? Try "integral2".

[Henan Fang]

@Torsten Thanks for your comment. However, we have tried "integral2", the two problems still exist, and are even more serious.

[Torsten]

lambda01 = 4.7175048762838045 * (11.630653268960174 - x) / (v - 0.4120716294548327)^(2/3);
lambdad1 = 4.7175048762838045 * (12.042724898415006 - v - x) / (v - 0.4120716294548327)^(2/3);

You divide by (v - 0.4120716294548327)^(2/3). This explains at least that you get problems when v is smaller or equal 0.4120716294548327 because results will be complex-valued or you get a division by 0.

And since you restrict integration to (x + y > 9.658875554160838 - v) & (x + y < 9.658875554160838), you define a discontinuous integrand for which you cannot expect high-precision results.

Did you try "vpaintegral" for your task ?

Answer 1

[Torsten]

That's the best you can get, I guess:

V = 0.46:0.01:1;
for i = 1:numel(V)
   v = V(i);
   lower_x = 9.578251022212591 - v;
   upper_x = 9.658875554160838 - v;
   lower_y = @(x)-x + 9.658875554160838 - v;
   upper_y = @(x)-x + 9.658875554160838;
   jr1 = integral2(@(x,y)arrayfun(@(x,y)TR(x,y,v),x,y),lower_x,upper_x,lower_y,upper_y,'AbsTol',1e-25,'RelTol',1e-25);
   lower_x = upper_x;
   upper_x = 9.658875554160838;
   lower_y = 0;
   upper_y = upper_y;
   jr2 = integral2(@(x,y)arrayfun(@(x,y)TR(x,y,v),x,y),lower_x,upper_x,lower_y,upper_y,'AbsTol',1e-25,'RelTol',1e-25);
   jr(i) = jr1 + jr2;
end
plot(V,jr)

function T = TR(x,y,v)
  k1 = (51.2317 * (x + 9.57825102221259));
  k2 = (51.2317 * (x + v - 9.57825102221259));
  rho1 = 23.5875 * (v - 0.4120716294548327)^(1/3);
  lambda01 = 4.7175048762838045 * (11.630653268960174 - x) / (v - 0.4120716294548327)^(2/3);
  lambdad1 = 4.7175048762838045 * (12.042724898415006 - v - x) / (v - 0.4120716294548327)^(2/3);
  Ai_lambda01 = airy(0, lambda01);
  AiP_lambda01 = airy(1, lambda01);
  Bi_lambda01 = airy(2, lambda01);
  BiP_lambda01 = airy(3, lambda01);
  Ai_lambdad1 = airy(0, lambdad1);
  AiP_lambdad1 = airy(1, lambdad1);
  Bi_lambdad1 = airy(2, lambdad1);
  BiP_lambdad1 = airy(3, lambdad1);
  term1 = rho1 .* AiP_lambda01 .* BiP_lambdad1 ./ k1 + k2 .* Ai_lambda01 .* Bi_lambdad1 ./ rho1 - rho1 .* AiP_lambdad1 .* BiP_lambda01 ./ k1 - k2 .* Ai_lambdad1 .* Bi_lambda01 ./ rho1;
  term2 = -Ai_lambda01 .* BiP_lambdad1 + k2 .* AiP_lambda01 .* Bi_lambdad1 ./ k1 - k2 .* Ai_lambdad1 .* BiP_lambda01 ./ k1 + AiP_lambdad1 .* Bi_lambda01;
  c1 = 4 .* pi^(-2) ./ (term1.^2 + term2.^2);
  T = double(k2.*abs(c1).^2./k1);
end

[Henan Fang]

@Torsten Thanks very much for your nice answer. It well solves the first problem.