0

In the following code

%% 
clc
clear 
%% 
M=4.5
lambda_0=0.332;
%% 
alpha2=0.988;
beta2=-0.05;
omega2=6;
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=aisecond(eta02);
airy_INT2=integral(@(n) airy(n),eta02,eta_inf2);
eta2=@(y) ((i*alpha2*lambda_0)^(1/3))*y+eta02;

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);

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)-Gi2(eta02)/airy(eta02)*airy(eta));

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

It is clear that, when eta=eta02, aa2*(Gi2(eta02)-Gi2(eta02)/airy(eta02)*airy(eta02))=0 by definition. However, it is not and I think it is due to the precision that MATLAB uses. What would be best in order to overcome this problem?

Peter
  • 71
  • 6
  • What exactly are you trying to accomplish? [This is likely relevant.](https://stackoverflow.com/q/686439/52738) – gnovice Jul 11 '17 at 16:09
  • Apparently, the operation that I indicated isn't zero, but should be. Since I am dealing with very large numbers, I wonder if the problem has to do with the accuracy of Gi2(eta02) and Gi2(eta02)/airy(eta02)*airy(eta02), so that, when substracting, the results isn't zero. How could I solve this problem? – Peter Jul 11 '17 at 16:16

0 Answers0