echo on

######################################
#                                    #
#   WPROWADZENIE DO OCTAVE/MATLABa   #
#                                    #
# (na podstawie MATLAB User's Guide) #
#                                    #
#        Piotr Krzyzanowski          # 
#            luty 2000               #
#                                    #
######################################

# macierz z klawiatury definiujemy wierszami
A = [1 2 3; 4 5 6 ; 7 8 9]

#lub wrecz w formacie macierzowym
A = [ 1 2 3
4 5 6
7 8 9]

#mozna tez wczytac macierz z pliku
load a.dat

#albo wygenerowac za pomoca skryptu (tutaj: gena.m)
gena

#aby sprawdzic, jakie mamy obiekty w pamieci
who


#aby zobaczyc, jak ktorys z nich wyglada (np. a)
a

#mozemy usunac z pamieci to, co niepotrzebne
clear a

#wspolrzedne wektorow moga byc wyliczone z wyrazen
x = [-1.3, sqrt(3), (1+2+3)*4/5]

#mozemy takze dodawac nowe wspolrzedne (zwrocmy uwage na dodany
#automatycznie zerowy x(4)
x(5) = abs(x(1))

#wektor mozemy dolaczyc do macierzy
r = [10 11 12];
A = [A; r]

#a z macierzy mozemy wyciac podmacierz
A_small = A(1:3,:)

#uzywamy go jako kalkulatora
1900/2000
sin(pi)
e^{i*pi}

#przenosimy dlugie wyrazenia
s = 1 + 2 + 3 ...
- 4 - 5 - 6

#i sprawdzamy, co mamy w pamieci
who
whos

#zapisujemy wszystkie obiekty do pliku
save
save -ascii zmienne.dat

#a nastepnie wszystko czyscimy w pamieci
clear

#aby zaraz odtworzyc z dysku
load zmienne.dat
who
A

#zmieniamy format wydruku
format long
A
format long e
A
format short
A

#jaki jest epsilon maszynowy i czy maszyna realizuje standard arytmetyki IEEE
eps
isieee
s = 1/0

#operacje macierzowe
A = A_small
B = A'
x = [1 2 3]'
C = A + B
y = x -1
x'*y
b = A*x
size(A)
size(x)
z = A\b
z = inv(A)*b
expm(A)
det(A)
trace(A)
poly(A)

#pare slynnych macierzy
H = hilb(6)
cond(H)
V = vander(z)
cond(V)

#dalej nt notacji dwukropkowej
z = 6:-1:1
x = (0:0.2:3)';
y = exp(-x).*sin(x);
[ x y]
k = linspace(-pi,pi,4)
A(:)

#algebra liniowa cd
p = [r 13]
P = compan(p)
roots(p)
eig(P)

A = rand(4,3)
tril(A)
triu(A)
diag(A)
size(A)

A = A(1:3,1:3);
[ L U ] = lu(A)
[Q R ] = qr(A_small)
poly(A_small)
roots(ans)
svd(A_small)
eig(A_small)

#definiujemy funkcje z linii komend
function y = humps(x)
y = 1./((x-.3).^2+.01) + 1./((x-.9).^2+.04) - 6;
endfunction

#i rysujemy uzywajac gnuplota
x = -1:.01:2;
plot(x,humps(x))

#calkujemy
q = quad('humps',0,1)

#i znajdujemy mce zerowe z przybl poczatk x_0 = 1.0
fsolve('humps',1)

#zapisujemy wszystkie nasze komendy do pliku historia_ddmmyy
history -w historia_ddmmyy

# aby zakonczyc sesje piszemy: quit

#rysujemy funkcje dwoch zmiennych
#      -2x^2 + 3xy   + 4 sin(y) - 6
#       3x^2 - 2xy^2 + 3 cos(x) + 4

function y = f (x)
  y(1) = -2*x(1)^2 + 3*x(1)*x(2)   + 4*sin(x(2)) - 6;
  y(2) =  3*x(1)^2 - 2*x(1)*x(2)^2 + 3*cos(x(1)) + 4;
endfunction

#liczba wezlow siatki
N = 20;
#rownomiernie rozlozone wezly w prostokacie [0,1]x[0,5]
x = linspace(0,1,N);
y = linspace(0,5,N);
#utworzenie siatki
[xx,yy] = meshgrid(x,y);
#wartosci w wezlach
for i = 1:N
  for j = 1:N
    z= f([xx(i,j);yy(i,j)]);
    zz1(i,j) = z(1); # pierwsza wspolrzedna 
    zz2(i,j) = z(2); # druga wspolrzedna
    zero(i,j) = 0.0; # plaszczyzna zerowa
  endfor
endfor

#rysujemy jedna po drugiej
mesh(xx,yy,zz1);
mesh(xx,yy,zero);
mesh(xx,yy,zz2);

#a teraz wszystko na jednym obrazku
hold off
hold on
mesh(xx,yy,zz1);
mesh(xx,yy,zero);
mesh(xx,yy,zz2);
hold off


#rozwiazujemy uklad rownan nieliniowych
#      -2x^2 + 3xy   + 4 sin(y) = 6
#       3x^2 - 2xy^2 + 3 cos(x) = -4


[x, info] = fsolve ("f", [1; 2])
perror ("fsolve", info)

#rozwiazujemy uklad rownan zwyczajnych
##
## The 'oregonator'. Przyklad dostarczany razem z Octave.
##
## Reference:
##
##   Oscillations in chemical systems.  IV.  Limit cycle behavior in a
##   model of a real chemical reaction. Richard J. Field and Richard
##   M. Noyes, The Journal of Chemical Physics, Volume 60 Number 5,
##   March 1974.

function dx = oregonator (x, t)

  dx = zeros (3, 1);

  dx(1) = 77.27*(x(2) - x(1)*x(2) + x(1) - 8.375e-06*x(1)^2);
  dx(2) = (x(3) - x(1)*x(2) - x(2)) / 77.27;
  dx(3) = 0.161*(x(1) - x(3));

end

t = linspace (0, 1500, 1500);
x0 = [1,4,1.01];
y = lsode ("oregonator", x0, t);

#tak nie nalezy ogladac rozwiazania
y
#lepiej tak!
plot(t,y);
#ewentualnie, przeskalowawszy, normujac wszystko do jednosci
plot(t,y*diag(1./max(abs(y))));

#likwidujemy oznaczenia "line1,..." (aby pojawily sie w zamian inne,
#trzeba wpisac stosowny tekst miedzy srednikami)
plot(t,y*diag(1./max(abs(y))),';;');

#dalej, drukujemy ten obrazek w Postscripcie
gset term postscript 
gset output "wykres.ps"
replot
#rownie dobrze moze byc 
#  gset term gif/latex; aby zobaczyc liste opcji,
#nalezy napisac gset term

