Why the integration can not be calculated out correctly?

2 次查看（过去 30 天）
Henan Fang
Henan Fang 2025-5-7
评论： Henan Fang 2025-5-14
Ran in:
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
  3 个评论
显示 1更早的评论
Henan Fang
Henan Fang 2025-5-9
@Torsten Thanks for your comment. However, we have tried "integral2", the two problems still exist, and are even more serious.
Torsten
Torsten 2025-5-9
编辑：Torsten 2025-5-9
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 ?

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2025-5-9
移动：Torsten 2025-5-9
Ran in:
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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by