[Octave-bug-tracker] [bug #34266] Four issues with "residue" and residuez"

[Bernard Grung]

Follow-up Comment #13, bug #34266 (project octave):

Jordi:

Sorry for the long delay -- I have been very busy.

I have just uploaded a file named "reply02.pdf" which includes some answers to
your questions. Also, it includes a discussion of the following two new Octave
test functions for investigating error tolerances using "bg_residue" and
"residue”.

Bernie 


%%% Error calculations using "bg_residue" 
z1=+7.0372976777e+06;
p1=-3.1415926536e+09;
p2=-4.9964813512e+08;
r1=-(1+z1/p1)/(1-p1/p2)/p2/p1  ;
r2=-(1+z1/p2)/(1-p2/p1)/p2/p1  ;
r3=(1+(p2+p1)/p2/p1*z1)/p2/p1  ;
r4=z1/p2/p1                    ; 
r=[r1;r2;r3;r4];
p=[p1;p2;0;0];
k=[];
e=[1;1;1;2];
b=[1,z1];
a=[1,-(p1+p2),p1*p2,0,0];
[br, ar] = bg_residue (r, p, k, e);
%%-------------------------------------------------------
fprintf('Start of error calculationsn');
fprintf('         eps = %+22.16en',eps);
assert (abs(ar-a)         < eps); %absolute error
assert (abs(1-a(1)/ar(1)) < eps); %relative error
%%-------------------------------------------------------
assert (abs(br(1)-b(1))   < eps*10);  %absolute error
assert (abs(1-b(1)/br(1)) < eps*10);  %relative error
fprintf('        b(1)=%+22.16en',b(1));
fprintf('       br(1)=%+22.16en',br(1));
fprintf('  b(1)-br(1)=%+22.16e  ',b(1)-br(1));
fprintf(       '<  eps*10=  %+22.16e for absolute errorn',eps*10);
fprintf('1-b(1)/br(1)=%+22.16e  ',1-b(1)/br(1));
fprintf(       '<  eps*10=  %+22.16e for relative errorn',eps*10);
%%-------------------------------------------------------
assert (abs(br(2)-b(2))   < eps*10^7);  %absolute error
assert (abs(1-b(1)/br(1)) < eps*10  );  %relative error
fprintf('        b(2)=%+22.16en',b(2));
fprintf('       br(2)=%+22.16en',br(2));
fprintf('  b(2)-br(2)=%+22.16e  ',b(2)-br(2));
fprintf(       '<  eps*10^7=%+22.16e for absolute errorn',eps*10^7);
fprintf('1-b(2)/br(2)=%+22.16e  ',1-b(2)/br(2));
fprintf(       '<  eps*10  =%+22.16e for relative errorn',eps*10  );
fprintf('End of error calculationsn');

%%% Error calculations using "residue" 
z1=+7.0372976777e+06;
p1=-3.1415926536e+09;
p2=-4.9964813512e+08;
r1=-(1+z1/p1)/(1-p1/p2)/p2/p1  ;
r2=-(1+z1/p2)/(1-p2/p1)/p2/p1  ;
r3=(1+(p2+p1)/p2/p1*z1)/p2/p1  ;
r4=z1/p2/p1                    ; 
r=[r1;r2;r3;r4];
p=[p1;p2;0;0];
k=[];
e=[1;1;1;2];
b=[1,z1];
a=[1,-(p1+p2),p1*p2,0,0];
[br, ar] = residue (r, p, k, e);
%%-------------------------------------------------------
fprintf('Start of error calculationsn');
fprintf('          eps=%+22.16en',eps);
ax=[a(2),a(3),a(4),a(5)]; % trim first coefficient
assert (abs(ar-ax)         < eps); %absolute error
assert (abs(1-ax(1)/ar(1)) < eps); %relative error
%%-------------------------------------------------------
bx=[z1]; % trim first coefficient
assert (abs(br(1)-bx(1))   < eps*10^9);  %absolute error
assert (abs(1-bx(1)/br(1)) < eps*100);  %relative error
fprintf('        bx(1)=%+22.16en',bx(1));
fprintf('        br(1)=%+22.16en',br(1));
fprintf('  bx(1)-br(1)=%+22.16e ',bx(1)-br(1));
fprintf(         '< eps*10^9=%+22.16e for absolute errorn',eps*10^9);
fprintf('1-bx(1)/br(1)=%+22.16e ',1-bx(1)/br(1));
fprintf(         '< eps*100 =%+22.16e for relative errorn',eps*100);
fprintf('End of error calculationsn');




    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?34266>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/