Clear Filters
Clear Filters

How to solve swallowtail integral

18 views (last 30 days)
洋斌 马
洋斌 马 on 2 May 2024
Commented: David Goodmanson on 29 May 2024 at 20:44
What do we do with this caustic integral?The simple way of writing an integral doesn't seem to compute it well.

Answers (1)

David Goodmanson
David Goodmanson on 3 May 2024
Edited: David Goodmanson on 4 May 2024
Hello YM,
You won't be able to solve this as a closed form function of X,Y,Z but you can find the integral for given values of X,Y,Z. It seems advantageous to split the integrand into cos(...) + i*sin(...) and do the integrals separately (the code below does not make the final step of multiplying by i).
The integrals involve both positive and negative s. For negative s one can substitute s --> -s so that
Int{-inf,0} f(s) ds --> Int{0,inf} f(-s) ds
and the domain of integration is always s>=0 in exchange for having two functions, f(s) and f(-s). Now divide the integrals into two regions, 0<=s<=b and b<=s<=inf. In the first region, do the integration per usual. In the second region make the substitution
s = u^(1/5), ds = (1/5) u^(-4/5)
For s<= b and b small enough the oscillations with exp(i*s^5 ...) are acceptable, and for s>b the oscillations go like exp(i*u ...) and are under control but one must stay away from the origin.
The values of b for which this method works appear to be fairly restricted. Also for upper limit for the integration, inf does not work but 1e5 does work. You will probably have to play around with the values. I am not sure how wide ranging the values of X,Y,Z can be, but smaller values seem to work all right. Sometimes you get a warning message when the tolerance is still quite acceptable. It would probably pay to increase the tolerance from the default but I did not do that.
If not for the Y term, (s.^5 + Z*s.^3 + Y*s.^2 + X*s) would be an odd function and the sine of that argument would be too. Consequently if Y=0, the final resuting sine integral = 0.
For the integration just from 0 to inf, there is a check available when X=Y=Z=0. Plugging in those values with b = 1, upper limit = 1e5 results in
Isinp =
ans =
0.283729745105399 % gamma(1/5)*sin(pi/10)/5
Icosp =
ans =
0.873230365517819 % gamma(1/5)*cos(pi/10)/5
with disagreement in the fifth decimal point for sine and the sixth decimal point for cosine.
X = 1;
Y = 2;
Z = 3;
a = 0;
b = 1;
c = 1e5;; % this limit applies to u, not s
% suffix p and m are for integrals that came from +s and -s respectively
% prefix f and g are for regions 0<=s<=b and b<=s<=inf respectively
fcosp = @(s) cos(s.^5 + Z*s.^3 + Y*s.^2 + X*s);
fcosm = @(s) cos(s.^5 + Z*s.^3 - Y*s.^2 + X*s);
gcosp = @(u) (1/5)*cos(u +Z*u.^(3/5) +Y*u.^(2/5) +X*u.^(1/5))./u.^(4/5);
gcosm = @(u) (1/5)*cos(u +Z*u.^(3/5) -Y*u.^(2/5) +X*u.^(1/5))./u.^(4/5);
Icosp = integral(fcosp,a,b) +integral(gcosp,b^5,c);
Icosm = integral(fcosm,a,b) +integral(gcosm,b^5,c);
Icos = Icosp + Icosm
fsinp = @(s) sin(s.^5 + Z*s.^3 + Y*s.^2 + X*s);
fsinm = @(s) -sin(s.^5 + Z*s.^3 - Y*s.^2 + X*s);
gsinp = @(u) (1/5)*sin(u +Z*u.^(3/5) +Y*u.^(2/5) +X*u.^(1/5))./u.^(4/5);
gsinm = @(u) -(1/5)*sin(u +Z*u.^(3/5) -Y*u.^(2/5) +X*u.^(1/5))./u.^(4/5);
Isinp = integral(fsinp,a,b) +integral(gsinp,b^5,c);
Isinm = integral(fsinm,a,b) +integral(gsinm,b^5,c);
Isin = Isinp + Isinm
Icos = 0.8397
Isin = -0.0751
% % check when X=Y=Z=0
% format long
% format compact
% Isinp
% gamma(1/5)*sin(pi/10)/5
% Icosp
% gamma(1/5)*cos(pi/10)/5
% format
洋斌 马
洋斌 马 on 29 May 2024 at 11:28
This is indeed a method. But it seems to have a large deviation from the correct result
David Goodmanson
David Goodmanson on 29 May 2024 at 20:44
could you reply with at least one example with values of X,Y,Z, and the correct result?

Sign in to comment.




Community Treasure Hunt

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

Start Hunting!