% inte08.m
% fd, 13 octobre 2008, 04 mars 2025  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
printf("calcul d une quatrieme integrale singuliere\n");
function  A=ss(x);   A = 1/sqrt(x) ;   endfunction 
function  A=ff(x);   A = 1/sin(sqrt(x+x**4+x**6)) ;   endfunction 
function  A=gg(x);   A = 1/sqrt(x) +  sqrt(x) / 6  ... 
+ 7*x*sqrt(x) / 360 - 269*x*x*sqrt(x)/540   ;   endfunction 
function  A=g1(x);   A = ff(x) - gg(x) ;   endfunction 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear  ; format long

printf("tentative de calcul plus precis de la comtribution analytique\n");
format long

n = 100 ;
dx = 1 / n ; 

xmin = 10**(-2) ; 
for j=1:n  ; 
yy(j) = j*xmin ; uu(j) = ff(yy(j)) ; 
vv(j) = -129600 * (ff(yy(j)) - gg(yy(j))) / ( 10819  *sqrt(yy(j)) * yy(j)**3)  ; 
ww(j) = g1 (yy(j))  / ( sqrt(yy(j)) * yy(j)**3)  ;  end 

% plot ( yy , uu  , '*b;;' ,  yy , vv , '@r;;' ) ;  pause
plot ( yy , vv , 'or;vv;' , yy , ww , '*k;xx;' ) ; 
title ('fonctions d erreur') ; pause (.7) 
print('08-graphe-fonction-ecart', '-djpg','-S891,630'); 

jjmax = 4
for jj = 1:jjmax 						 

if jj==1 ; n = 10 , end 
if jj==2 ; n = 100 , end 
if jj==3 ; n = 1000 , end 
if jj==4 ; n = 10000 , end 

%%% n = 2000  ,  dx = 1 / n ; 
dx = 1 / n ; 

xx = zeros(n+1,1) ;
for j=1:n+1 ;  xx(j) = (j-1) * dx  ; end 

%                                                    trapezes  
zi4 = 37357 / 18900  ; 
zi4 = zi4 + dx * ( g1(xx(2)) +  g1(xx(n+1))  )  / 2    ; 
for j=3:n ;  
zi4 = zi4 + dx * g1(xx(j)) ; % zi4  
end  
%                                                     Simpson 
zi5 =  37357 / 18900  ;  
zi5 = zi5 +  dx * ( g1(xx(2)) + g1(xx(n+1)) ) / 6    ; 
for j=3:n ;   zi5 = zi5 + dx * g1(xx(j)) / 3 ; end  
for j=2:n ;   
xxx = (xx(j)+xx(j+1))/2 ; 
zi5 = zi5 + 2 * dx * g1(xxx) / 3 ; end  

printf("trapezes "); zi4   
printf("Simpson  "); zi5

end 
