disp("Rozwiazanie kolokwium ON 7-maj-2008 gr 1015");;

clear all
function y=G(s)
#funkcja do ktorej odwrotna szukamy
# dla s\in[0,1)
 y=2*(-s-2*log(1-s));
endfunction
global YY

function y=GY(s)
#Funkcja y=exp(-G(s))-exp(-YY)
# exp(-G(s)) latwiejsza dla s bliskiego 1  
 global YY 
 y=exp(0.5*s).*(1-s)-exp(-YY/4);
end

function y=GYY(s)
#funkcja y=G(s)-YY pomocnicza
 global YY 
 y=G(s)-YY;
end

z=linspace(0,100,1e3);
N=length(z);
S=zeros(1,N);
S(1)=0;
SS=S;
for k=2:N,
 YY=z(k);
 S(k)=fsolve("GY",S(k-1));
 #SS(k)=fsolve("GYY",SS(k-1));
endfor

#sprawdzamy czy rzeczywiscie odwrotna
printf("Blad w normie inf : z-G(S): %g \n",\
norm(G(S)-z,"inf"));

plot(z,S,";Funkcja odwrotna do 2*(-s-2*log(1-s));");
#czy S=G^{-1}(z)<1?
printf("czy S<1?(0 nie 1 tak) %d\n", max(S-1)<0); 
#ewentualnie
printf("Inaczej czy S<1: %d\n",sum(S>=1)==0); 

#aby narysowa tylko wykres mozna inaczej
#zh = G^{-1}(100)
YY=100;
zh=fsolve("GYY",1-1e-9);
printf("Czy G(zh)=100? G(zh)-100=%g\n",G(zh)-100);
#kompletnie zle czyli tak sie nie da....
#ale metoda zgadniecie (albo bisekcji jak ktos wiec co to jest...
zh=1-8.17e-12;
printf("zgadujemy G(%g)=%g ? czy = 100?\n",zh,G(zh));
#blisko choc nie idealnie
#czy uzywajac GY
zh=fsolve("GY",0.99);
#sprawdzamy
printf("rozwiazujac rownanie z exp(-G(%g))=exp(-100) \n G(%g)=%g\n",zh,zh,G(zh));
s=linspace(0,zh);
hold on
plot(G(s),s,";wykres f odwrotnej do G inaczej;");
#tak czy siak trzeba uzyc exp(-G(s)) z fsolve....
#w ostatecznosci zgadnac (czyli reczna bisekcja) ewent zaprogramowac bisekcje
#ale to tylko da wykres ...



#Zadanie 2
##############################################################
disp("Zad 2");
global CC=1e-9;
function y=f(x)
global CC
y=CC*cos(2008*sqrt(x));
endfunction
N=1e7;
t=rand(N,1);
#wazne aby wektorowo wywolac f(t) 
mc=sum(f(t))/N

#jeszcze 2razy (to dodatkowo)
t=rand(N,1);mc1=sum(f(t))/N
t=rand(N,1);mc2=sum(f(t))/N
#rzad ten sam i znak ale nawet 1 cyfry dokl nie ma
#trzeba wziasc wiecej liczb

#a quad
c=quad("f",0,1)
#niby dobrze ale te 1e-9 przed funkcja niepokoi
#zmienmy tolerancje
c1=quad("f",0,1,[0 1e-8])
#czyli znow klopoty
#scalkujmy bez stalej 
CC=1;
c2=quad("f",0,1)*1e-9
#tu juz quad wychwycil ze ma klopoty

#mozna podzielic odcinek na podocinki
N=10;
#dobieramy N recznie az nie bedzie komunikatu o bledzie
t=linspace(0,1,N);
c3=0;
CC=1;
for k=2:N,
  c3+=quad("f",t(k-1),t(k));
endfor
c3=c3*1e-9
#no tu i znak i rzad jak w MCarlo 
hold off



