% Matlab-skripti: ..v/3/01/matlab/suuntak1.m
% Kts. myös exasuuntak.m
% HA 4.9.2001 
%
% Aihe: Diffyhtälön  y'=f(x,y)  suuntakenttä ja isokliinit.
%
% Havainnollinen versio, jossa jananpätkät piirretään for-
% silmukassa. Oikeaoppisempi tehoversio: suuntak.m
%
clf
f=inline('y-x','x','y')     % f(x,y)=y-x
%f=inline('x.*y','x','y')   % Toinen esimerkki: f(x,y)=xy

h=0.5;a=-1;b=1;c=0;d=2;   % Muuta tarpeen mukaan
x=a:h:b;y=c:h:d;                          

[X,Y]=meshgrid(x,y);
Z=f(X,Y);

 flipud(Z)   % havainnollistukseen, tosi käytössä kommentti.

 plot(X(:),Y(:),'x')   % Jos haluat merkitä ristillä pisteet, joissa
                        % derivaatta-arvot on laskettu, poista kommentti
hold on
xx=X(:); yy=Y(:); n=length(xx);
for i=1:n
   v=[1,f(xx(i),yy(i))];v=0.5*h*v/norm(v);
   alkup=[xx(i),yy(i)]-v;loppup=[xx(i),yy(i)]+v;
   plot([alkup(1),loppup(1)],[alkup(2),loppup(2)])
end;

% Jos haluat vain suuntakentän, lopeta tähän.
% 
% Piirretään edellä olevia korkeusarvoja vastaavat isokliinit.

x=linspace(a,b,30);y=linspace(c,d,30);
[X,Y]=meshgrid(x,y);
ZZ=f(X,Y);
 contour(x,y,ZZ,Z(:))  % Korkeusvektoriksi valitaan Z(:), tällöin           
axis equal             %  saadaan kaikki suuntakenttäpiirroksen arvot,
axis([a,b,c,d])        %    ovelaa!
grid
title('y''=f(x,y)')     % Jos haluat oikean otsikon, editoi tähän
shg