# 
# Nonlinear (almost linear) systems studied by looking at the linear
# system "closest to" the nonlinear one.
# 
# Art Belmonte
# Mon, 27/May/96
# Math 308-509 [Maple V Release 4]
# Section 9.3: Almost Linear Systems
# 
# T-456/1a: Critical points of the nonlinear system.
> with(DEtools): with(linalg): with(student): readlib(mtaylor):
Warning, new definition for D
> F:=(x, y)->x - y + x*y; G:=(x, y)->3*x - 2*y - x*y;

                      F := (x, y) -> x - y + x y


                    G := (x, y) -> 3 x - 2 y - x y

> dx_dt:=F(x, y); dy_dt:=G(x,y);

                         dx_dt := x - y + x y


                       dy_dt := 3 x - 2 y - x y

> derivs:=[dx_dt, dy_dt];

               derivs := [x - y + x y, 3 x - 2 y - x y]

> eqs:=equate(derivs, 0); solve(eqs, {x, y});

            eqs := {x - y + x y = 0, 3 x - 2 y - x y = 0}


                  {y = 0, x = 0}, {y = 1/3, x = 1/4}

> deq:=equate([diff(x(t), t), diff(y(t), t)],
>      derivs(t));

        d
deq := {-- y(t) = 3 x(t) - 2 y(t) - x(t) y(t),
        dt

    d
    -- x(t) = x(t) - y(t) + x(t) y(t)}
    dt

# T-456/1b: Here's Maple's phase portrait. (Also use MATLAB for further
# exploration.)
> inits:={[0, 0.1, 0.1], [0, 0.5, 0.5]};

                 inits := {[0, .5, .5], [0, .1, .1]}

> DEplot(deq, [x(t), y(t)], t=-5..5, x=-1..1, y=-1..1, inits,
> scaling=constrained, stepsize=0.05);

> 
# T-456/1c:
# 1. Linear parts of F and G (as expressions).
# 2. Formulation of coefficient matrix A of the corresponding linear
# system.
# 3. Investigation of the nature of the critical point (0, 0) of the
# said linear system: Since the eigenvalues of A are complex conjugates
# with negative real parts, we have from Table 9.3.1 (T-451) that (0, 0)
# is an asymptotically stable spiral point (as we surmised from the
# phase portrait above). 
> LF:=mtaylor(F(x,y), [x, y], 2);

                             LF := x - y

> LG:=mtaylor(G(x,y), [x, y], 2);

                           LG := 3 x - 2 y

> A:=genmatrix([LF, LG], [x, y]);

                                 [1    -1]
                            A := [       ]
                                 [3    -2]

> lambda:=eigenvals(A); evalf(");

                                    1/2                 1/2
           lambda := - 1/2 + 1/2 I 3   , - 1/2 - 1/2 I 3


      -.5000000000 + .8660254040 I, -.5000000000 - .8660254040 I

# T-456/1d: Verification that original system is almost linear.
> g:=evalm(derivs - [LF, LG]);

                           g := [x y, -x y]

> g_polar:=subs(x=r*cos(theta), y=r*sin(theta), evalm(g));

              [ 2                          2                      ]
   g_polar := [r  cos(theta) sin(theta), -r  cos(theta) sin(theta)]

> map(Limit, evalm(g_polar / r), r=0, right);

[  lim    r cos(theta) sin(theta),   lim    -r cos(theta) sin(theta)]
[r -> 0+                           r -> 0+                          ]

> `|g(x,y)| / |(x,y)| -> 0`:=value(");

                  |g(x,y)| / |(x,y)| -> 0 := [0, 0]

# 
>