An Autonomous Linear Nonhomogeneous System

Here is an example of a linear, autonomous, nonhomogeneous 2X2 systems of ordinary differential equations.

This notebook was created by T. Shelton, Southwestern University to supplement a handout for  Differential Equations II, updated for Mathematica 11.1.

reference  Differential Equations: Theory, Technique, and Practice, Simmons & Krantz, McGraw-Hill Co, Inc. 2007

dx/dt=f[x,y]
dy/dt=g[x,y]

Nullclines, Direction Field, Equilibria

Clear[f,g]
f[x_,y_]:= 2x - 2y + 10
g[x_,y_]:=11x-8y+49

nullclines

Solve[f[x,y]==0,y]

Solve[g[x,y]==0,y]

equilibria

Solve[{f[x,y]==0,g[x,y]==0}]

graphing window containing interesting places

Clear[xmin,xmax,ymin,ymax]
xmin=-5.0;
xmax=1.0;
ymin=-2.0;
ymax=5.0;

plots of nullclines and equilibria

Clear[np1,np2,np,pts]
np1 =Plot[5+x,{x, xmin, xmax},PlotStyle->{Red,Dashed,Thick}];
np2 = Plot[(1/8) (49+11 x),{x, xmin, xmax},PlotStyle->{Blue,Dashed,Thick}];
np = Show[{np1,np2},PlotRange->{{xmin,xmax},{ymin,ymax}}];
pts = Graphics[{PointSize[Large],Point[{-3,2}]}];
Show[np,pts]

the direction field

Clear[vPlot]
vPlot = VectorPlot[{f[x,y],g[x,y]},{x,xmin,xmax},{y,ymin,ymax}];

all together

Show[np,pts,vPlot]

Solution Technique 1:  Use the Eigensystem Method and a nonhomogeneous solution technique

dx/dt=f[x,y], dy/dt=g[x,y]  where f[x_, y_] := 2 x - 2 y + 10, g[x_, y_] := 11 x - 8 y + 49.
In vector form,  (({{x'(t)}, {y'(t)}}))  =  (({{2, -2}, {11, -8}})) (({{x(t)}, {y(t)}}))  + (({{10}, {49}}))

complementary solution

Eigensystem[{{2,-2},{11,-8}}]

check for linear independence of the first solution, written as a row vector,...

Expand[ (E^((-3-(Sqrt[3]))t)){(8/11)+(1/11) (-3-(Sqrt[3])),1}]

... and the second solution, written as a row vector,...

Expand[(E^((-3+(Sqrt[3]))t)){(8/11)+(1/11) (-3+(Sqrt[3])),1}]

here is the matrix whose determinant is the Wronskian of the two solutions

MatrixForm[Transpose[ {{(5/11) (E^((-3-(Sqrt[3])) t))-(1/11) (Sqrt[3]) (E^((-3-(Sqrt[3])) t)),(E^((-3-(Sqrt[3])) t))},{(5/11) (E^((-3+(Sqrt[3])) t))+(1/11) (Sqrt[3]) (E^((-3+(Sqrt[3])) t)),(E^((-3+(Sqrt[3])) t))}}]]

Det[ Transpose[ {{(5/11) (E^((-3-(Sqrt[3])) t))-(1/11) (Sqrt[3]) (E^((-3-(Sqrt[3])) t)),(E^((-3-(Sqrt[3])) t))},{(5/11) (E^((-3+(Sqrt[3])) t))+(1/11) (Sqrt[3]) (E^((-3+(Sqrt[3])) t)),(E^((-3+(Sqrt[3])) t))}}]]

since this non-zero constant multiple of an exponential function is never zero, our two solutions are linearly independent.  Here is our complementary solution:

Expand[c1 (E^((-3-(Sqrt[3]))t)){(8/11)+(1/11) (-3-(Sqrt[3])),1} + c2 (E^((-3+(Sqrt[3]))t)){(8/11)+(1/11) (-3+(Sqrt[3])),1}]

particular solution

Undetermined Coefficients seems the simplest because our initial guess is a constant vector for our particular solution, which is linearly independent from any portion of our (exponential-based) complementary solution.

(({{xp'(t)}, {yp'(t)}}))  =  (({{2, -2}, {11, -8}})) (({{xp(t)}, {yp(t)}}))  + (({{10}, {49}}))  where (({{xp(t)}, {yp(t)}}))=(({{a1}, {a2}})), so  (({{xp'(t)}, {yp'(t)}}))  = (({{0}, {0}})).

(({{2, -2}, {11, -8}})) (({{xp(t)}, {yp(t)}}))  + (({{10}, {49}}))= (({{0}, {0}})) or (({{2, -2}, {11, -8}})) (({{xp(t)}, {yp(t)}}))  = - (({{10}, {49}}))
or   (({{xp(t)}, {yp(t)}}))  =(((({{2, -2}, {11, -8}})))^(-1))(({{-10}, {-49}})) if the matrix inverse exists

Det[{{2,-2},{11,-8}}]

since the determinant is nonzero, its inverse exists.

Inverse[{{2,-2},{11,-8}}]

{{-(4/3),(1/3)},{-(11/6),(1/3)}}.{-10,-49}

Thus, our particular solution is {-3, 2}.

general solution

Expand[c1 (E^((-3-(Sqrt[3]))t)){(8/11)+(1/11) (-3-(Sqrt[3])),1} + c2 (E^((-3+(Sqrt[3]))t)){(8/11)+(1/11) (-3+(Sqrt[3])),1} + {-3,2}]

Solution Technique 2:  Use a change of variables

We decide on a change of variables, X(t) = x(t) - x0, Y(t) = y(t) - y0,  to transform our system into a homogenous one.

dx/dt=f[x,y], dy/dt=g[x,y]  where f[x, y]= 2 x - 2 y + 10, g[x, y]= 11 x - 8 y + 49.

Clear[x0,y0]

f[x-x0,y-y0]

Expand[%]

g[x-x0,y-y0]

Expand[%]

(2 x - 2 y) + ( 10 - 2 x0 + 2 y0) == (2 x - 2 y) iff  10 - 2 x0 + 2 y0 == 0
(11 x - 8 y) + (49 - 11 x0 + 8 y0) == 11 x - 8 y iff 49 - 11 x0 + 8 y0 == 0

Solve[{10-2 x0+2 y0==0,49-11 x0+8 y0==0}]

Our change of variables is x(t) = X(t) - 3, y(t) = Y(t) + 2.  Note that dX/dt = dx/dt and dY/dt = dy/dt. 
In our new coordinate system, we have dX/dt= 2 X - 2 Y,  dY/dt= 11X - 8 Y. 
In vector form,  (({{X'(t)}, {Y'(t)}}))  =  (({{2, -2}, {11, -8}})) (({{X(t)}, {Y(t)}}))

Eigensystem[{{2,-2},{11,-8}}]

Expand[c1 (E^((-3-(Sqrt[3]))t)){(8/11)+(1/11) (-3-(Sqrt[3])),1} + c2 (E^((-3+(Sqrt[3]))t)){(8/11)+(1/11) (-3+(Sqrt[3])),1}]

The solution above is X(t), Y(t).  However, x(t) = X(t) - 3, y(t) = Y(t) + 2.

Expand[c1 (E^((-3-(Sqrt[3]))t)){(8/11)+(1/11) (-3-(Sqrt[3])),1} + c2 (E^((-3+(Sqrt[3]))t)){(8/11)+(1/11) (-3+(Sqrt[3])),1} + {-3,2}]

Trajectories from the General Solution

Expand[c1 (E^((-3-(Sqrt[3]))t)){(8/11)+(1/11) (-3-(Sqrt[3])),1} + c2 (E^((-3+(Sqrt[3]))t)){(8/11)+(1/11) (-3+(Sqrt[3])),1} + {-3,2}]

The initial value of the solution at t=0

%/.{t->0}

tableOfPts = {};
Do[Do[AppendTo[tableOfPts,{x0,y0}],{x0,xmin+.5,xmax-.5,2.5}],{y0,ymin+.5,ymax-.5,2.5}]
tableOfPts

Clear[initialPts];
initialPts =  Graphics[{PointSize[Large],Point[tableOfPts]}];
Show[initialPts,Axes->True, PlotRange->{{xmin,xmax},{ymin,ymax}}]

For each initial point chosen, solve for the constants to obtain a unique solution, and add it to a list of solution curves.

Clear[numPts,tableOfFcns,rule,i];
numPts = Length[tableOfPts];
tableOfFcns={};
For[i=1,i≤numPts,i++,
rule=Solve[{-3+((5 c1)/11)-(((Sqrt[3]) c1)/11)+((5 c2)/11)+(((Sqrt[3]) c2)/11) == tableOfPts[[i,1]],2+c1+c2== tableOfPts[[i,2]]}];
AppendTo[tableOfFcns,
{-3+(5/11) c1 (E^((-3-(Sqrt[3])) t))-(1/11) (Sqrt[3]) c1 (E^((-3-(Sqrt[3])) t))+(5/11) c2 (E^((-3+(Sqrt[3])) t))+(1/11) (Sqrt[3]) c2 (E^((-3+(Sqrt[3])) t)),2+c1 (E^((-3-(Sqrt[3])) t))+c2 (E^((-3+(Sqrt[3])) t))}/.rule]]
tableOfFcns

Here are the trajectories, clearly beginning at the chosen initial points.

Show[  {
ParametricPlot[tableOfFcns,{t,0,2},PlotStyle->Thick],
initialPts}]

Trajectories, Nullclines, Equilibria, and the Direction Field

Show[  {np,pts,VectorPlot[{f[x,y],g[x,y]},{x,-5.25,2},{y,-2,8}],
ParametricPlot[tableOfFcns,{t,0,2},PlotStyle->Thick],
initialPts}, PlotRange->{{-5.25,2},{-2,8}}]