echo on
clear all;#removing all variables
#some examples of simpl calls t ononlinear solvers
#using fzero

x=fzero(@sin,[3 3.2])

[x, FVAL, INFO, OUTPUT] = fzero (@sin,[3 3.2])

f=@(x) x*x-2;
x = fzero (f,[1 2])
[x, FVAL, INFO, OUTPUT] = fzero (f,[1 2])


#using fsolve
#univariate function

x=fsolve(f,2)

#multivariable
f2d=@(x) [x'*x-1;x'*[1;-1]];
x=fsolve(f2d,[1;1])

function [y,dy]=g(x)
 y=x*x-2;
 if(nargout>1)
  dy=2*x;
  disp("computes derivative");
 endif
endfunction


fsolve(@g,2)
#computing analytic Jacobian
fsolve(@g,2,optimset("Jacobian","on"))

#same for 2d


function [y,dy]=g2d(x)
 y=[x'*x-1;x'*[1;-1]];
 if(nargout>1)
  dy=[2*x';[1 -1]];
  disp("computes 2d derivative");
 endif
endfunction

fsolve(@g2d,[2;1])
#computing analytic Jacobian
fsolve(@g2d,[2;1],optimset("Jacobian","on"))
#FVEC=F(X);INFOerror code (1 if OK) 
#OUTPUT -structure with some info FJAC= approximation of jacobian at X 
 [X, FVEC, INFO, OUTPUT, FJAC] = fsolve(@g2d,[2;1]) 



#implicit function - getting explicit ....
#we want to compute values of y(x) given by F(x,y)=0
#at discrete mesh e.g .in order to plot the graph...

function y=explicit(F,x,y0)
#getting approx of function y(x) given in an explicit way: F(x,y)=0
#at the mesh points in x
#INPUT 
#F - function handle to a two varable scalar function F(x,y) \in R
#x- mesh ordered in increasing way i.e. x(k)>x(k-1)
#y0 - value of y at x(1) 
#OUTPUT
#y - vector with approx. values of y(x)
#Example
#F=@(x,y) x*x+y*y-1
#y=explicit(F,linspace(0,1-eps/2),1);
y(1)=y0;
for k=2:length(x),
 [y(k),Fxy,info]=fsolve(@(y) F(x(k),y),y(k-1));
 if(info!=1) 
   disp("Warning: fsolve failed at x=");disp(x(k)); return; 
 endif
endfor
endfunction

#tests (trivial)

F=@(x,y) x*x+y*y-1
y=explicit(F,x=linspace(0,0.99),1);
plot(x,y,";implicit f. y(x): x^2+y^2-1=0;")
pause(1);


#inverse function  ....
#we want to compute values of x(y)  given by x->y(x)
#at discrete mesh e.g .in order to plot the graph...
#naturally x(y) is the inverse of y(x)

function x=invfunc(f,y,x0)
#getting approx of the inverse function x(y) for a known y(x)=f(x)
#at the mesh points in y
#INPUT 
#f - function handle to a univariate scalar function f(x) \in R
#y- mesh ordered in increasing way i.e. y(k)>y(k-1)
#x0 - value of x at y(1) 
#OUTPUT
#x - vector with approx. values of x(y)
#Example
#x0=0.8*pi; y0=cos(x0);x=invfunc(f=@cos,y=linspace(y0,0.9),x0);plot(y,x)
x(1)=x0;
for k=2:length(y),
 [x(k),ymfx,info]=fsolve(@(x) y(k)-f(x),x(k-1));
 if(info!=1) 
   disp("Warning: fsolve failed at y=");disp(y(k)); return; 
 endif
endfor
endfunction

#testujemy
#acos:
x0=0.8*pi; y0=cos(x0);x=invfunc(f=@cos,y=linspace(y0,0.9),x0);plot(y,x,";our acos;")
pause(1);
plot(y,acos(y)-x,";error acos- our acos;");
pause(1);
#onother function
f=@(x) 1.1*x+sin(x)
x=invfunc(f,y=linspace(0,10),0);
x1=linspace(min(x),max(x));y1=f(x1);
plot(y,x,";computed;",y1,x1,";inverse by reflection;")
pause(1)
plot(x1,y1,";f;");figure(2);plot(y1,x1,";inverse by reflection;")


echo off
