%% Harj. 3    R/harj3ratk.m
% 28.3.2012
%% Teht. 1 Heiluri
% Muunnetaan systeemiksi
%{
  Merk: x=Theta, y=Theta'=x' =>
  x' = y
  y' = -(g/L)*sin(x)

Kirjoitetaan tiedostoon heiluri.m
-------------------------------------
function [yp] = heiluri(t,Y)
% Heiluri-systeemin m"a"arittely
%   
x=Y(1);y=Y(2);
g=9.81;L=g;
yp=[y;
    -(g/L)*sin(x)];
end
-------------------------------------
Testataan:
heiluri(0,[0;1])
%}
%% Lasketaan, ode45 
close all
[T,Y]=ode45(@heiluri,[0 40],[0 2]);
subplot(2,1,1)
plot(T,Y)
legend('Heilahduskulma','Kulmanopeus')
title('Aikataso')
shg
subplot(2,1,2)
plot(Y(:,1),Y(:,2))
axis equal;
title('Faasitaso')
%%
%figure
close all
Tvali=[0 10]
[T1,Y1]=ode45(@heiluri,Tvali,[1; 1]);
[T2,Y2]=ode45(@heiluri,Tvali,[-5; 2]);
[T3,Y3]=ode45(@heiluri,Tvali,[5; -2]);

plot(Y1(:,1),Y1(:,2),Y2(:,1),Y2(:,2),Y3(:,1),Y3(:,2))
axis equal; axis([-5 5 -3 3])
%%
hold on
Tvali=[0 40]
[T,Y]=ode45(@heiluri,Tvali,[0; 2]);
plot(Y(:,1),Y(:,2))
%%
[T,Y]=ode45(@heiluri,Tvali,[0; 2-.1]);
plot(Y(:,1),Y(:,2),'r')
[T,Y]=ode45(@heiluri,Tvali,[0; 2+.1]);
plot(Y(:,1),Y(:,2),'g')

%% Suuntakentt"a
x=-5:.5:5; y=-3:.5:3;
[X,Y]=meshgrid(x,y);
dx=Y; dy=-sin(X);   % Derivaatat heiluriyhtalosta
quiver(X,Y,dx,dy)

%%
%pplane8


%% Teht. 3 RA-teht. "shooting method"

% y''=y^2 -1
%{
 y1=y, y2=y'=y1'

  y1'=y2
  y2'=y1^2 -1

Funktio shoot.m
-----------------
function [yp] = shoot(t,y)
% Diffyhtaloryhma, kuten odexx sen haluaa.
%   
  %y1'=y2
  %y2'=y1^2 -1
yp=[y(2);y(1).^2-1];  
end
------------------
%}

%% Testataan:
[t y1]=ode45(@shoot,[0 0.5 1],[0;1.4])  
% Ollaan kiinnostuneita aikapisteesta 1.
%  Kokeilun jalkeen huomataan, etta dy(0)=1.4 ampuu lahelle
%  y:n arvoa y(1)=1
%% 
% reuna_arvo-funktio
% ----------------------------------
%{
function [y1] = reuna_arvo(dy0)
%  syote: derivaatta-arvo
%  tulos: y(1) 
[t y]=ode45(@shoot,[0 0.5 1],[0;dy0]);
y1=y(end,1);
end
%} 
%
%% Milla derivaatan alkuarvolla osutaan reuna-arvoon 1.
%  Pikku briljeeraus: Maaritellaan komentorivilla nimeton
%  funktio, jossa reuna-arvon arvosta vahennetaan 1.
%  (Oltais voitu vahentaa jo reuna_arvofunktiossa.)
%%

dy0=fzero(@(x)reuna_arvo(x)-1,1)

%% Ratkaistaan alkup. tehtava, missa derivaatan alkuarvona ylla ratkaist dy0.

[T,Y]=ode45(@shoot,[0 1],[0;dy0]);
plot(T,Y(:,1));shg
