Problem nr. 5

There are some nice solutions like ... Nobody seems to have an elegant solution for the direction field. After writing the text for L3 I found out the "correct solution" that is elegant, symmetric and efficient. (!! -:) )

File dir2.m

% Direction field for autonomous 2 by 2 systems m-file (script) 
% 3.2.1998 by HA (modified from Bäckström's 1 eq. case
% Requires the system m-file to be written to work for an initial value
% matrix (with 2 rows) instead of a column vector. This form is perfectly
% good for the ode-solvers
%
clear variables; hold off;
xmax=10; ymax=3; delx=xmax/20; dely=ymax/20;   % *** problem dependent
x=-4:delx:xmax; y=-ymax:dely:ymax;
[X,Y]=meshgrid(x,y);
Z=sys5(0,[X(:),Y(:)]');
DX=Z(1,:);DY=Z(2,:);
DX=reshape(DX,size(X));DY=reshape(DY,size(Y));
  L=sqrt(DX.^2+DY.^2); % Initial lengths
  L=L+1e-10;
  DX1=DX./L; DY1=DY./L; % Unit length for all arrows
figure(1); quiver(X,Y,DX1,DY1,'g'); grid on;title('Direction field');
% axis([x1 x2 y1 y2]);

File sys5.m

The trick is to write the file defining the system in a slightly more general form than is needed for the ode45 etc. functions. Instead of requiring the initial values as a one column vector of length 2 we will allow it to be given as several column vectors, i.e. (2 by "several") matrix. Of course it has to be written in an array smart fashion. (This doesn's show in the code below as there are no "mixed" terms.)

After having played with meshgrid, reshape and Mat(:) you should now be able to read the above dir2-code as an open book.

function u=sys5(t,y)
%
%global k c;
k=1;c=1/2; % Make these local while testing 
u(1,:)=y(2,:);
u(2,:)=-c*y(2,:)-k*sin(y(1,:));

Combining direction field and interactive soltution curves

Having done the dir2-commands in Matlab you can do many nice things like:

hold on
init=ginput(1);[T,Y]=ode45('sys5',[0,10],init);plot(Y(:,1),Y(:,2),'g')
axis([5 7 -1 1])
init=ginput(1);[T,Y]=ode45('sys5',[0,20],init);plot(Y(:,1),Y(:,2),'r')

for i=1:5;
  init=ginput(1);plot(init(1),init(2),'o');[T,Y]=ode45('sys5',[0,20],init);
  plot(Y(:,1),Y(:,2),'r')
end;

Write dir2 as a function

Now it is an easy matter to do this. I let you do it though. Then one can easily change the grid density, lengths of arrows etc.