Numerical integrals and MATLAB precision

The following code

    clc
    clear
    M=4.5
    lambda_0=0.332;
    betats=0.05;
    betan=betats/(lambda_0^(5/4));
    f=@(x,y) y^(1/3)*(y^2+x^2)-1.001*sqrt((1-M^2)*y^2 + x^2);
    alphan=real(fsolve(@ (y) y^(1/3)*(y^2+betan^2)-1.001*sqrt((1-M^2)*y^2 + betan^2),3));
    alphats=alphan*(lambda_0^(5/4));
    omegan=2.299*(alphan^(2/3));
    omegats=(lambda_0^(3/2))*omegan;
    gammats=sqrt((M^2-1)*alphats^2-betats^2);
    gammats=-complex(0,1)*gammats;
    beta1=0;
    alpha1=1;
    omega1=6;
    gamma1=sqrt((M^2-1)*alpha1^2-beta1^2);
    alpha2=alpha1-alphats;
    beta2=beta1-betats;
    omega2=omega1-omegats;
    gamma2=sqrt((M^2-1)*alpha2^2-beta2^2);
    N=25;
    YMAX=150;
    eta02=-i*omega2/((i*alpha2*lambda_0)^(2/3));
    eta_inf2=((i*alpha2*lambda_0)^(1/3))*YMAX+eta02;
    syms lu
    aiprime(lu) = airy(1,lu);
    aisecond(lu)= diff(airy(1,lu));

    airy_D2=airy(1,eta02);
    airy_DD2=vpa(aisecond(eta02));
    airy_INT2=integral(@(n) airy(n),eta02,eta_inf2);
    aa2=2*pi*complex(0,1)*gamma2*beta2*lambda_0*airy_D2/(alpha2*(alpha2^2+beta2^2)*airy_INT2-gamma2*lambda_0*airy_D2*(i*alpha2*lambda_0)^(2/3));
    bb2=-aa2*airy(2,eta02)*airy_INT2/airy(eta02);

    ai2=@(x) x*airy(x);
    bi2=@(x) x*airy(2,x);

    eta2=@(y) ((i*alpha2*lambda_0)^(1/3))*y+eta02;
    Gi2=@(x) -(airy(2,x)*integral(@(n) airy(n),eta_inf2,x)-airy(x)*integral(@(n) airy(2,n),eta02,x));
    W2=@(eta) aa2*Gi2(eta)+bb2*airy(eta);

    W2VEC=arrayfun(W2,arrayfun(eta2,0:0.01:N));
    plot(abs(W2VEC),0:0.01:N)

produces a noisy `W2VEC`. I think there is something very wrong either in the numerical integration or in the Airy functions.

The constants `aa2` and `bb2` have been chosen so that the two terms cancel each other at `0`, so a good result will verify `W2VEC(1)=0`.

What is going on here? Is it the accuracy of the integrals?

[carlos g]

Any idea? It seems to me that it has to do with the precision for the integral calculation being lower than the magnitude of the numbers involved (10^15). Is this right? If so, how could I overcome such problem?

[Karan Gill]

The code is hard to read. If there was some explanation, that would help.

• Try checking the value of intermediate expressions
• High-precision integration is available using the vpasolve function in Symbolic Math Toolbox.