

#solvery nieliniowe octave'a - kilka przykladow 

g=@cos; #uchwyt (wsk) do funkcji cos

x=fzero(g,[1,2])
y=fsolve(@cos,1.5)
x-pi/2
y-pi/2
x-y

#f anonimowa

f=@(x) [x'*x-1;[1,1]*x-0.5]

x= fsolve(f,-[1;1.1])
f(x)


#odwracanie funkcji f(x) na al be
#f=cos na [eps,pi*(1-eps)]
function y=f(x)
 y=cos(x);
end

eps=0.01;
al=eps;
be=pi*(1-eps);
a=f(be);
b=f(al);
global yy;

function z=fmy(x)
global yy
z=f(x)-yy;
endfunction

y=linspace(a,b,N=100);
x=zeros(size(y));
x(1)=be;
x(N)=al;
for k=2:N-1,
 #rozwiazujemy rownanie f(x)-y(k)=0 startujac z x(k-1)
 #ktore rozwiazuje f(x)-y(k-1)=0
 yy=y(k);
 x(k)=fsolve(@fmy,x(k-1)); 
endfor
norm(acos(y)-x,'inf')

plot(y,acos(y),";acos;",y,x,"r;nasz acos;");
pause(1);
plot(y,x-acos(y),";blad;")
pause(1);
#dodatkowo  przez odwrocenie osi
s=linspace(al,be);
plot(cos(s),s,";odw osie;",y,x,";nasz acos;",y,acos(y),";acos;")
pause(1);
#jeszcze raz przez odwracanie osi blisko lewego konca odcinka ale na [0,pi]
fr=0.5;
K=round(fr*N);
plot(cos(s(1:K)),s(1:K),"+",y(N-K+1:N),acos(y(N-K+1:N)),"*;acos;")
#widac ze odwracanie osi daje troche inne punkty.... 


#mozemy teraz narysowac wykres f(x)=1./(sin(1)+eps-sin(x)*x) na [f(eps),f(1)] 
global eps
eps=1e-2;
function y=f(x)
 global eps
 y=1./(sin(1)+eps-sin(x).*x);
end
al=eps;
be=1;
a=f(al);
b=f(be);
global yy;

function z=fmy(x)
global yy
z=f(x)-yy;
endfunction

N=500;
x=zeros(1,N);
if(b<a)
 y=linspace(b,a,N);
 x(1)=be;x(N)=al;
else 
 y=linspace(a,b,N);
 x(N)=be;x(1)=al;
endif

for k=2:N-1,
 #rozwiazujemy rownanie f(x)-y(k)=0 startujac z x(k-1)
 #ktore rozwiazuje f(x)-y(k-1)=0
 yy=y(k);
 x(k)=fsolve(@fmy,x(k-1)); 
endfor
norm(f(x)-y,'inf')
s=linspace(al,be,N);
fs=f(s);
plot(y,x,"r;f^{-1};",fs,s);
pause(1);

#blisko lewego konca 
K1=round(1);K2=round(0.03*N);
a1=x(K1);a2=x(K2);
plot(fs(K1:K2),s(K1:K2),"+;f(s) vs s;");
pause(1);
plot(y(K1:K2),x(K1:K2),"*;nasz f^{-1};");
pause(1);
#policzmy przyblizenie pochodnej 
M=length(s);h=s(2)-s(1);
df=(fs(2:M)-fs(1:M-1))/h;[max(abs(df)),min(abs(df))]
h=y(2)-y(1);dfi=(x(K1+1:K2)-x(K1:K2-1))/h;[max(abs(dfi)),min(abs(dfi))]


#blisko prawego konca 
K1=round(N*0.95);K2=round(N);
a1=x(K1);a2=x(K2);
plot(fs(K1:K2),s(K1:K2),"+;f(s) vs s;");
pause(1);
plot(y(K1:K2),x(K1:K2),"*;nasz f^{-1};");
pause(1);
#policzmy przyblizenie pochodnej 
M=length(s);h=s(2)-s(1);
df=(fs(2:M)-fs(1:M-1))/h;[max(abs(df)),min(abs(df))]
h=y(2)-y(1);dfi=(x(K1+1:K2)-x(K1:K2-1))/h;[max(abs(dfi)),min(abs(dfi))]

