How to solve swallowtail integral
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
What do we do with this caustic integral?The simple way of writing an integral doesn't seem to compute it well.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1684646/image.png)
0 Kommentare
Antworten (1)
David Goodmanson
am 3 Mai 2024
Bearbeitet: David Goodmanson
am 4 Mai 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 =
0.283749732315802
ans =
0.283729745105399 % gamma(1/5)*sin(pi/10)/5
Icosp =
0.873231080653766
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
2 Kommentare
David Goodmanson
am 29 Mai 2024
could you reply with at least one example with values of X,Y,Z, and the correct result?
Siehe auch
Kategorien
Mehr zu Numerical Integration and Differentiation finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!