echo on
########################################################################
#
#	Przykladowy skrypt w octavie
#
#	Calkowanie w octavie (przy okazji wprowadzimy zmienne globalne)
#		funkcja quad()
#	
#	Rozwiazywanie rownan nieliniowych f: fsolve()
#
#	Proste przyklady
#
########################################################################

#
#  na poczatek wielomianu aby funkcja miala parametr uzyjemy zmiennej 
#  globalnej
#########################################################################
global p
p=4
function y=jedn(x)
global p
y=x^p;
endfunction
#lepiej wyswielac dokladniejsze wyniki
format long e

#Liczymy calke wartozwrocic uwage na sposb podania funkcji jako 
# napis z nazwa funkcji ktora musi byc zdefiniowana jako funkcja 1 argumentu
#zwracajaca jeden argument (skalarny)
#np tak jak my powyzej
###########################################################################
quad("jedn",0,1)

#Mozemy okreslic dokladnosc (to tylko teoretyczna estymacja)
quad("jedn",0,1,1e-1)

#Mozemy dostac wiecej danych - IER - kod bledu - 0 jak OK, 1,2,etc jak zle
#NFUN -ilosc wywolan funkcji
#ERR - estymacja bledu
p=0.5;[V, IER, NFUN, ERR] = quad('jedn',0,1,1)
#Sprawdzamy
V- 1/(p+1)
#B dobrze mimo ze funkcja nawet nie Lipschitzowska
#Sprobujmy f z mocna osobliwosci±
function y=os(x)
if(x!=0) 
  y=1./x;
else 
  y=10^7;
endif  
endfunction

#Przy okazji warto nauczyc operacji logicznych
# rownosc ==, roznosc !=, ! zaprzeczenie, <,>,<=,>= - to latwo zgadnac, && -and , || - or
#np
a=10;
if((a<20)&&(a>-20)) disp("a w (-20,20)"); 
else disp("a poza (-20,20)");
endif
a=30;
if((a<20)&&(a>-20)) disp("a w (-20,20)"); 
else disp("a poza (-20,20)");
endif
#przy czym prawda jest liczba rozna od zera, falsz 0 np
1 && 0
1 || 0
!(34)
!0
#Mozna tez uzywac jako negacji ~ np
~4
2~=4
#lub w wersji skroconej & i |
2&3
2&0
#t orowniez dziala dla wektorow 

#A czy operacje logiczne dzialaja dla macierzy?
format short; 
x=[ -2 -1 3 0 4]
#sprawdzmy czy elementy dodatnie
ox=x>0
#teraz co bedzie z negacja
~ox
#ox i x
ox && x # wyszlo dziwnie
#okazuje sie ze &&, i II tylko dla liczb dla wektorow
ox & x 
ox | x

#Zadanie napisac funkcje ktora WEKTOROWO (tzn bez petli!)
#oblicza funkcje 'daszel' tzn f parzysta; f(x)=0 dla x poza [-1,1]
#i f(x)=1-x na [0,1] 


#Wracamy do calek
[y,IER]=quad("os",-1,1)
#mamy blad! Powinno wyjsc zero
# [1e-3,1e-4] - to wektor tolerancji bledu bezwzglednej i wzglednej
#Podajemy potencjalne osobliwosci
[y,IER]=quad("os",-1,1,[1e-3,1e-4],0)
#I jest OK!

#Warto pamietac o funkcji eval - wykonuje polecenia ze stringu bedacego arg np
eval("a=2;...
      a*a");


#A teraz niech Panstwo zaimplementuja w octavie 
#funkcje obliczajaca norme L2 - calkowa! (trzeba zdefiniowac funkcje w funkcji
#przy pomocy eval())
#A nastepnie zmodyfikuja ja tak
#aby liczyla norme Lp tzn nazwa funkcji i p maja byc jedynymi parametrami
#jesli tego nie bylo na wykladzie to to jest dosc trudne (sztuczkowe) ale do zrobienia!

#Zadanie: napisac funkcje liczaca dla danej funckji postaci y=f(x)
#calke \int_0^1 f(x)x^p dx dla roznych p

#Zadanie: korzystajac z poprzedniego zadania i macierzy Hilberta
#znalezc wiel. najlepszej aproksymacji stopnia n w normie L^2(0,1) 
#dla fukcji charaterystycznej
#odcinka [0,1/2] - policzyc macierz Gramma dla (1,x,x...,x^k) w L^2(0,1) 
#porownac z macierza Hilberta (generowana funkcja hilb(k+1)!


#Teraz rozwiazmy proste  rownanie nieliniowe
# x^2 - a =0
global a 
a=2;
function y=fkw(x)
global a
 y= x.^2 -a;
endfunction
#Wystartujemy z a
x0=a;
fsolve('fkw',x0)

#to ponizej powinno dzialac ale nie dziala w octave 3.0.0!
#dlatego wykomentowano
#Mozna podac pochodna  - pyt czy fsolve uzywa met Newtona? Sa tez inne metody
#wykorzystujace pochodna....
#function y=dfkw(x)
#global a
# y= 2*x - a;
#endfunction
#fsolve(["fkw";"dfkw"],a)

#No i zadanie dla Panstwa - sprawdzic jak dziala fsolve dla funkcji niewypuklej
#(ani nie wkleslej) np x.*(1.2+sin(x)) (jak mozna to definiujmy funkcje
#wektorowo)czy zbiega dla dowolnego x0 np
#-12,-4,0.5,3 itp i to w wersji z pochodna i bez. Pytanie dlaczego nie ma
#zbiezosci? (Wsk prosze narysowac wykres na [-10,10])


