########################################################
#
#
# Rozwiazanie (przykladowe) zadan z kolokwium z Obliczen
#  nakowych 26/04/2006 830-10 - czas 60min
#  polecenia rysujace sa wykomentowane
#
#
########################################################

disp("rozwiazanie kolokwium obl naukowe 05-06");


#Definiujemy pole wektorowe
#dx/dt=(2- x-2y)*x
#dy/dt=(2-2x-y)*y

global A;
A=[1 2;2 1];

function y=g(x,t)
#pole wektorowe
global A;
y=([2;2]- A*x).*x;
end

#Zad 1 
#Rozwiazac zag poczatkowe dla t w [0,3] z h=0.05 
#Narysowac wykresy (t,x(t)) i (t,y(t) na 1 rysunku
#oraz wykres trajektorii (x,y) na kolejnym rysunku
#Znalezc max x i y

#definujemy dyskretny wektor czasu na [0.3]
#z h=0.05
t=0:0.05:3;
#Rozwiazyjemy z war poczatkowym x(0)=2*y(0)=0.2
X=lsode("g",[0.2;0.1],t);

#Obliczamy max obliczonego x i y na [0.3]
max(X)

#wykresy wzgl czasu na 1 rysunku
#plot(t,X,[";x(t);";";y(t);"])
#pause(5);

#a teraz trajektoria - krzywa fazowa tzn wykres (x,y)
#plot(X(:,1),X(:,2),";x vers y;")

#Zad 2:
# kilka trajektorii na 1 rysunku dla danych war poczatkowych
t=0:0.05:10;X1=lsode("g",[0.05;0.2],t);
X2=lsode("g",[0.3;0.05],t);
X3=lsode("g",[0.2;0.3],t);
X4=lsode("g",[0.05;0.04],t);
#pause(5);
#plot(X1(:,1),X1(:,2),";war pocz x=0.05 y=0.2;",X2(:,1),X2(:,2),X3(:,1),X3(:,2),X4(:,1),X4(:,2));

#Zad 3
#Szukamy pkt stacjonarnych (poza osiami)
#i potem obliczamy Jakobian w tym punkcie - ze wzru albo 
#przyblizony jakobian przez roznice dzielone z krokiem delta=1e-4

#Mozna rozwiazujac uklad rownan liniowych (bo F=([2;2]-Ax).*x
x0=A\[2;2]
#albo wprost korzystajac z funkcji fsolve:
#potrzebujemy funkcji y=F(x) takiej ze F(x)=f(x,t)
function y=F(x)
y=g(x,0);
endfunction

#rozwiazujemy rownanie F(x)=0 startujac z [0.5;1]
y0=fsolve("F",[0.5;1]);

#i porownujemy
printf("Blad fsolve kontra f-Ax=0 %g\n",norm(x0-y0,2));

#Jakbina z wzoru (na kartce liczymy wzory)
DF=[-x0(1) -2*x0(2); -2*x0(1) -x0(2)]

#a przyblizony przez roznice dzielone
delta=1e-4;
k1=(F(x0+delta*[1;0])-F(x0))/delta;
k2=(F(x0+delta*[0;1])-F(x0))/delta;
dF=[k1 k2];
printf("Blad DF wzor - dF- roznice dzielone dla delta=%g  %g\n",delta,norm(DF-dF,"inf"));


#Zad 4
#obliczmy wartosci wlasne i wektory wlasne tego jakobianu
[V,D]=eig(DF)
[v,d]=eig(dF)
#sa male roznice ale pomijalne

#Zad 5
#to jest siodlo! (pkt niestabilny)
