Linear Differential Equations (continued):

General Strategy from last time:

Example: Solve the system  ({{{, ({{x'=x+y, (1)}, {y'=4x+y, (2)}})}}) using DSolve[..].

By the analysis done in class last time, we know that (0,0) is a saddle. The eigenvalues of the coefficient matrix A=(({{1, 1}, {1, 4}}))  are 3 and -1 (real and with opposite signs), with corresponding eigenvectors v=(({{1}, {2}}))  and w=(({{-1}, {2}})), so the general solution of the system of differential equations is
(({{x(t)}, {y(t)}})) = (Subscript[c, 1]) (e^(3t))(({{1}, {2}})) + (Subscript[c, 2]) (e^(-t))(({{-1}, {2}})), where (Subscript[c, 1]) and (Subscript[c, 2]) are constant parameters. 
There are two special solutions when (Subscript[c, 1])=0, and respectively when (Subscript[c, 2])=0. The two solutions are (({{x(t)}, {y(t)}})) = (Subscript[c, 2]) (e^(-t))(({{-1}, {2}})) and (({{x(t)}, {y(t)}})) = (Subscript[c, 1]) (e^(3t))(({{1}, {2}})). For any value of t, the first solution remains of the line spanned by the eigenvector w (shown in red), while the second solution remains on the line spanned by the eigenvector v (shown in orange). These is why they look like straight lines when we plot the Phase Portrait of the differential equation using the command StreamPlot[..].

See the lecture notes for the code used to produce the Phase Portrait shown below:

DSolve

Since a linear system with constant coefficients can be solved explicitly, we can use the command DSolve[..]

DSolve[eqn,y,x]
solves a differential equation for the function y, with independent variable x.
DSolve[eqn,y,{x,xmin,xmax}]
solves a differential equation for the function y, with the independent variable x between xmin and xmax.
DSolve[{eqn1,eqn2,…},{y1,y2,…},…]
solves a list of differential equations for y1, y2, ... .

DSolve returns results as lists of rules. This makes it possible to return multiple solutions to an equation. For a system of equations, possibly multiple solution sets are grouped together inside curly brackets { }. You can use the rules to substitute the solutions into other calculations.

How to use the output of DSolve[...]

Clear[x,y,t]
sol1=DSolve[{x'[t]==x[t]+y[t], y'[t]==4x[t]+y[t]},{x,y},t]

When the second argument of DSolve is specified as y, the solution is returned as a transformation rule y->Function[...], where the right hand side of the rule is a pure function (a function with no name). A general solution contains arbitrary parameters labeled C[i] that can be varied to produce particular solutions for the equation.

f=Function[{t},(1/2) (E^(-t)) (1+(E^(4 t))) C[1]+(1/4) (E^(-t)) (-1+(E^(4 t))) C[2]]   (*copy paste, not advisable!*)

How to use the solution returned by DSolve[...] elegantly

g= x /. Flatten[sol1]
h= y /. Flatten[sol1]

The command Flatten[...] will get rid of the additional curly brackets. Instead of using two command lines, we can use only one:

{g,h} = {x,y} /. Flatten[sol1]

To obtain a particular solution, we can give values to the constants C[1] and C[2].

?C

{g0,h0} = {x,y} /. Flatten[sol1] /. {C[1]-> 0.1, C[2]->0.1}

We can plot the solution that we have obtained using ParametricPlot[...].

ParametricPlot[{g0[t],h0[t]}, {t,-1,0.1}, AxesOrigin->{0,0}]

Several other ways to call DSolve[...]

sol2=DSolve[{x'[t]==x[t]+y[t], y'[t]==4x[t]+y[t]},{x[t],y[t]},t]

We can also relabel the parameters generated by DSolve[...], see the example below and the Mathematica documentation on Generated Parameters.

sol3=DSolve[{x'[t]==x[t]+y[t], y'[t]==4x[t]+y[t]},{x[t],y[t]},t,GeneratedParameters->(Subscript[d,#]&)]

We can no longer use the same command as before to extract the solutions, {x,y} /. Flatten[sol2], because the replacement rule is now of the form x[t]→Rule, and not of the form x→Rule. Notice also that we defined the functions g2 and h2 using “=” and not “:=”.

Clear[t];
g2[t_]= x[t] /. Flatten[sol2]
h2[t_]= y[t] /. Flatten[sol2]

When an adequate number of initial conditions are specified, DSolve[...] returns particular solutions to the given equations.

sol4=DSolve[{x'[t]==x[t]+y[t], y'[t]==4x[t]+y[t],x[0]==2,y[0]==0},{x[t],y[t]},t]

{x1[t_],y1[t_]}= {x[t] ,y[t]}/. Flatten[sol4]

ParametricPlot[{x1[t],y1[t]},{t,-2.5,1.5}]

Nonlinear Differential Equations:

Consider the system  ({{{, ({{x'= f(x,y), (1)}, {y'= g(x,y), (2)}})}})  where f and g are functions of two variables x and y. Assume that f and g are not both linear.

The functions x and y are functions of time t and the derivatives are also taken with respect to the variable t.  ({{{, ({{x'(t)=f(x(t),y(t)), (1)}, {y'(t)=g(x(t),y(t)), (2)}})}}) .

A linear system  ({{{, ({{x'= ax+by, (1)}, {y'= cx+dy, (2)}})}})  whose coefficient matrix (({{a, b}, {c, d}})) is invertible has only one equilibrium point (critical point, steady state solution) at (x,y)=(0,0). Analyzing the stability of this equilibrium solutions gives complete information about the Phase Portrait of the ODE (the trajectories).

Unlike a linear system, a nonlinear system could have none, one, two, three, or any number of critical points. Like a linear system, however, the critical points are found by setting x’ = y’ = 0, and solve the resulting system f(x,y)=0 and g(x,y)=0. 
Since there might be multiple critical points present on the phase portrait, each trajectory could be influenced by more than one critical point. This results in a much more chaotic appearance of the phase portrait. Consequently, the type and stability of each critical point need to be determined locally (in a small neighborhood on the phase plane around the critical point in question) on a case-by-case basis. The trajectories will appear very differently than those of the linear systems, except very near the critical points.

General Strategy :

We approximate the nonlinear system by a linear system near the equilibrium points. 
Suppose that (a, b) is an equilibrium point, that is f(a, b) = g(a, b) = 0.

({{{, ({{x'= f(x,y) ≃ f(a.b) + (∂_x)f(a,b)(x-a) + (∂_y)f(a,b)(y-b)  = (∂_x)f(a,b)(x-a) + (∂_y)f(a,b)(y-b) , (1)}, {y'= g(x,y) ≃ g(a,b) + (∂_x)g(a,b)(x-a) + (∂_y)g(a,b)(y-b)  = (∂_x)g(a,b)(x-a) + (∂_y)g(a,b)(y-b) , (2)}})}})

The critical point can be translated to 0 via the coordinate change x := x - a and y := y - b and still retains all of its properties. The approximated system becomes a linear system:

({{{, ({{x'= (∂_x)f(a,b)x + (∂_y)f(a,b)y, (1)}, {y'=(∂_x)g(a,b)x + (∂_y)g(a,b)y , (2)}})}})   with coefficient matrix  (Subscript[J, (a,b)])=(({{(∂_x)f(a,b),  (∂_y)f(a,b)}, {(∂_x)g(a,b), (∂_y)g(a,b)}})) .

The matrix (Subscript[J, (x,y)])=(({{(∂_x)f(x,y),  (∂_y)f(x,y)}, {(∂_x)g(x,y), (∂_y)g(x,y)}}))is called the Jacobian matrix of the system. It contains the first order partial derivatives of f and g evaluated at the point (x,y).

Stability of equilibrium points:

We can use the eigenvalues of the matrix Jacobian matrix to decide the type and stability of the critical point (a,b). Let λ and μ be the two eigenvalues. Assume that λ and μ are different from 0.

({{Eigenvalues, Linear System, Nonlinear System}, {λ ,μ   Real, □, □}, {λ >μ>0, Nodal Source (Unstable ), Nodal Source (Unstable)}, {λ <μ<0, Nodal Sink(Stable ), Nodal Sink( Stable)}, {λ >0>μ, Saddle point (Unstable), Saddle point (Unstable)}, {λ =μ>0, Degenerate Source or Nodal Source(Unstable) 
depending of the geometric multiplicity of λ ,   Source (Degenerate, Nodal, Spiral Source)(Unstable)
depending of the geometric multiplicity of λ and the nonlinear terms}, {λ =μ<0, Degenerate Sink or Nodal Sink(Stable)
depending of the geometric multiplicity of λ, Sink (Degenerate, Nodal, Spiral Sink)(Stable)
depending of the geometric multiplicity of λ and the nonlinear terms}, {λ ,μ   Complex Conjugate, □, □}, {Re(λ)>0, Spiral Source (Unstable), Spiral Source (Unstable)}, {Re(λ)<0, Spiral Sink (Stable), Spiral Sink (Stable)}, {Re(λ)=0, Center (Stable), Center, Spiral Sink ,Spiral Source
(Stability cannot be determined based on the linear part alone)}})

Hyperbolic Equilibrium Points:

The equilibrium is said to be hyperbolic if all eigenvalues of the Jacobian matrix have non - zero real parts.

Hyperbolic equilibria are robust: Small perturbations do not change qualitatively the phase portrait near the equilibria. Moreover, local phase portrait of a hyperbolic equilibrium of a nonlinear system is equivalent to that of its linearization. This statement has a mathematically precise form known as the Hartman - Grobman Theorem.

Examples of nonlinear differential equations:

Example. Solve the system  ({{{, ({{x'= x(1-y), (1)}, {y'= y(x-1), (2)}})}})

Usually an explicit solution cannot be found, so DSolve will not work, but one can ask Mathematica to find a numeric solution using the NDSolve command.

DSolve[{x'[t]==x[t](1-y[t]), y'[t]==y[t](x[t]-1),x[0]==1,y[0]==1},{x[t],y[t]},t]

NDSolve[{x'[t]==x[t](1-y[t]), y'[t]==y[t](x[t]-1),x[0]==2,y[0]==1},{x[t],y[t]},{t,0,10}]

The most reliable information here comes from the Vector Field Plot.

StreamPlot[{x(1-y),y(x-1)},{x,-5,5},{y,-5,5},StreamScale->Automatic,StreamPoints->{{{{4,1},Red},Automatic}}]

What is the nature of the equilibrium solutions?

1. Set x' = 0 and y' = 0. In this case, it is easy to solve the equation by hand:  x(1-y)=0 and y(x-1)=0  implies (x,y)=(0,0) or (x,y)=(1,1). For more complicated cases, we will use Solve[..] or NSolve[..]

EqPoints=Solve[{x(1-y)==0, y(x-1)==0}, {x,y}]

2. Next we find the Jacobian matrix of the system. In this case, it is easy to compute it by hand: 
   (Subscript[J, (x,y)])=(({{(∂_x) (x(1-y)), (∂_y) (x(1-y))}, {(∂_x) (y(x-1)), (∂_y) (y(x-1))}}))=(({{1-y, -x}, {y, x-1}}))

In more complicated cases, we will ask Mathematica to compute the Jacobian matrix for us, as follows :

The commands D[f[x], {x}] or f'[x]  give the derivative of f with respect to x. If f is a function of several variables, D[f[x, y], {x}] or (∂_y)(f[x,y]) gives the partial derivative of f with respect to x. Try them out!

(∂_y)(x y+ y^2)

D[y^2, {y}]

The following code computes the Jacobian matrix:

Clear[f,J,x,y]
f[x_,y_]:= {x(1-y), y(x-1)}
J[f_]:=
Module[{A=Table[0,2,2]},
A[[1,1]]=D[f[x,y][[1]],{x}];
A[[1,2]]=D[f[x,y][[1]],{y}];
A[[2,1]]=D[f[x,y][[2]],{x}];
A[[2,2]]=D[f[x,y][[2]],{y}];
Return[A]]
J[f]

3. We evaluate the Jacobian matrix at the two equilibrium points and find its eigenvalues:

B=J[f] /. Flatten[EqPoints[[1]]];
Print["The Jacobian Matrix at the equilibrium point (0,0) is ", MatrixForm[B]]
Print["The eigenvalues of B are ", Eigenvalues[B]]

T=J[f] /. Flatten[EqPoints[[2]]];
Print["The Jacobian Matrix at the equilibrium point (1,1) is ",MatrixForm[ T]]
Print["The eigenvalues are ", Eigenvalues[T]]

4. Conclusions:

By analyzing the Real and Imaginary Part of the Eigenvalues, we conclude that the point (0, 0) is a saddle (the Jacobian matrix has two real eigenvalues, one smaller than 0, one greater than 0) whereas the point (1, 1) is a center (the Jacobian matrix has two complex conjugate eigenvalues with Real part equal to 0, and the Phase Portrait of the nonlinear equation depicts closed circular trajectories). This gives us information about the shape of trajectories in a vicinity of the equilibrium solutions.

Exercise :

Consider the system  ({{{, ({{x'= y, (1)}, {y'=-y-2sin(x), (2)}})}}).  Study the nature of the critical points (x, y) of the system for which -2π < x,y < 2π. Use only Mathematica commands, even when computing the critical points!