Harj. 7 teht. 3
================

"Rabbits and foxes"

r'=  2r - a r f
f'= -f +  a r f

Ensin systeemin määräävä .m-tiedosto:

rf.m
====

function yp=rf(t,y)
r=y(1);f=y(2);a=0.01;
yp=[2*r - a*r.*f;
    -f+a*r.*f];
---------------------------------

>> y0=[300;150];[T,Y]=ode23('rf',0,10,y0);plot(T,Y);figure;plot(Y(:,1),Y(:,2)
>> grid;hold on;
>> y0=ginput;[T,Y]=ode23('rf',0,10,y0);plot(Y(:,1),Y(:,2),'r')
(RTN, kun hiirikursori on kuvaikkunassa vie osoitetun (ja klikatun) paikan 
koordinaatit muuttujaan y0.)
Yllä oleva on erittäin mukava ja vuorovaikutteinan tapa piirtää hiirellä
osoitettuihin kohtiin faasitasokäyriä, erityisen mukava X:ssä.
(Windows:ssa joutuu aina kaivamaan kuvan esiin piilosta.)

Maple:ssa taspainopisteet ja linearisointi:

> f1:=2*r - a* r* f;

                                f1 := 2 r - a r f
--------------------------------------------------------------------------------
> f2:=-f+a*r*f;

                                f2 := - f + a r f
--------------------------------------------------------------------------------
> with(linalg):
--------------------------------------------------------------------------------
> J:=jacobian([f1,f2],[r,f]);

                                [ 2 - a f    - a r   ]
                           J := [                    ]
                                [   a f    - 1 + a r ]
--------------------------------------------------------------------------------
> krp:=solve({f1=0,f2=0},{r,f});

                    krp := {f = 0, r = 0}, {r = 1/a, f = 2/a}
Jos a=0.01, niin (100,200), kuten kuvastakin näkyy,

--------------------------------------------------------------------------------
> J1:=subs(krp[2],op(J));'

                                       [ 0  -1 ]
                                 J1 := [       ]
                                       [ 2   0 ]
--------------------------------------------------------------------------------
> eigenvects(J1);

                          2                           2
                [RootOf(_Z  + 2), 1, {[ 1, - RootOf(_Z  + 2) ]}]
--------------------------------------------------------------------------------
> eigenvals(J1);

                                   1/2       1/2
                                I 2   , - I 2
--------------------------------------------------------------------------------
> 
Puhtaasti imag, siis keskus
===========================