Phase diagrams for x' = Ax where a is 2x22008-04-01restart;with(LinearAlgebra):with(DEtools):with(plots):
<Text-field style="Heading 1" layout="Heading 1">Linear algebra with Maple</Text-field>Row vector:<1 | 2 | 3>;Column vector:<1, 2, 3>;A matrix is a "column" of row vectors:< <1 | 2 | 3>, <4 | 5 | 6> >;The same matrix may be viewed as a "row" of column vectors:< <1, 4> | <2, 5> | <3, 6> >;Maple performs multiplication on vectors and matrices by the "." (dot) operator which is not the same as the "*" (asterisk) operator.The "*" operator is reserved for commutative products, such as a*b = b*a where a and b are numbers, polynomials, etc.The "." operator is reserved for non-commutative product, such as A.B where A and B are matrices or A . x where A is a matrix and x is a vector.For exampleA := < <1, 2> | <3, 4> >;B := < <5, 6> | <7, 8> >;A . B;B . A;x := <9, 10>;A . x;The determinant:Determinant(A);The inverse:A^(-1);The transpose:A^%T;OK, we are done with this intro. Here I am going to "undefine" the variable A, B, x so that they won't interfere with the rest of the worksheet:A := 'A'; B := 'B'; x := 'x';
<Text-field style="Heading 1" layout="Heading 1">Saddle</Text-field>Here we view the coefficient matrix A a given by its diagonal factorization, that is A = P L P^-1. We specify L and P more or less arbitrarily and then compute A.This is sort of backward. Normally things are done the other way around, that is, A is given, and one computes L and P, but then the eigenvalues and eigenvectors generally turn out to be "ugly" numbers. The purpose of this backward approach is to assure that A has "nice" eigenvalues and eigenvectors.L := DiagonalMatrix([-2,3]);P := <<1,1>|<-1,2>>;A := P . L . P^(-1);Eigenvectors(A);tmp := A . < x(t), y(t) >;de1 := diff(x(t),t) = tmp[1];de2 := diff(y(t),t) = tmp[2];ic := [x(0)=0, y(0)=1], [x(0)=0, y(0)=2];DEplot([de1, de2], [x(t), y(t)], t=0..1, [ic], x=-5..5, y=-5..5, stepsize=0.1, linecolor=blue, thickness=1);ic := seq([x(0)=0, y(0)=h], h=-5..5, 0.5);DEplot([de1, de2], [x(t), y(t)], t=-1..1, [ic], x=-5..5, y=-5..5, stepsize=0.1, linecolor=blue, thickness=1);ic := seq([x(0)=0, y(0)=h], h=-5..5, 0.5), seq([x(0)=h, y(0)=0], h=-5..5, 0.5);DEplot([de1, de2], [x(t), y(t)], t=-1..1, [ic], x=-5..5, y=-5..5, stepsize=0.1, linecolor=blue, thickness=1, arrows=none); p1 := %:Here we plot two lines along the eigenvector directions. We express the lines in parametric form. Maple interprets [a(t), b(t), t=t1..t2] as the parametric specification of a curve.plot([[t,t,t=-5..5], [-t,2*t, t=-5..5]], scaling=constrained, thickness=2); p2 := %:display([p1,p2]);
<Text-field style="Heading 1" layout="Heading 1">Stable node, example 1</Text-field>Here we produce a matrix A = P L P^-1 by specifying L and P. See the remarks under the Saddle heading for the rationale behind this.The eigenvalues are both negative therefore we will have a stable node. The eigenvectors, that is the column vectors of P, are chosen here to be almost (but not quite) 90 degrees apart, so that the resulting phase portrait is not too distorted from what one would see in the (y1,y2) plane.Further down in example 2 we will see a case where the angle between the eigenvectors is smallish, therefore the phase portrait is severely distorted relative to the one in the (y1, y2) plane.L := DiagonalMatrix([-2,-3]);P := <<1,1>|<-1,2>>;A := P . L . P^(-1);Eigenvectors(A);tmp := A . < x(t), y(t) >;de1 := diff(x(t),t) = tmp[1];de2 := diff(y(t),t) = tmp[2];ic := [x(0)=1, y(0)=5], [x(0)=-1, y(0)=5];DEplot([de1, de2], [x(t), y(t)], t=0..2, [ic], x=-5..5, y=-5..5, stepsize=0.1, linecolor=blue, thickness=1);ic := seq([x(0)=h, y(0)=5], h=-5..5, 1.0), seq([x(0)=h, y(0)=-5], h=-5..5, 1.0),
seq([x(0)=-5, y(0)=h], h=-5..5, 1.0), seq([x(0)=5, y(0)=h], h=-5..5, 1.0);DEplot([de1, de2], [x(t), y(t)], t=0..2, [ic], x=-5..5, y=-5..5, stepsize=0.1, linecolor=blue, thickness=1); p3 := %:plot([[t,t,t=-5..5], [-t,2*t, t=-5..5]], scaling=constrained, thickness=2, color=red); p4 := %:display([p3,p4]);
<Text-field style="Heading 1" layout="Heading 1">Stable node, example 2</Text-field>See the remark under example 1.L := DiagonalMatrix([-1,-3]);P := <<2,1>|<1,2>>;A := P . L . P^(-1);Eigenvectors(A);tmp := subs(x=x(t), y=y(t), A . <x,y>);de1 := diff(x(t),t) = tmp[1];de2 := diff(y(t),t) = tmp[2];ic := [x(0)=1, y(0)=5], [x(0)=-1, y(0)=5];DEplot([de1, de2], [x(t), y(t)], t=0..3, [ic], x=-5..5, y=-5..5, stepsize=0.1, linecolor=blue, thickness=1);ic := seq([x(0)=h, y(0)=5], h=-5..5, 1.0), seq([x(0)=h, y(0)=-5], h=-5..5, 1.0),
seq([x(0)=-5, y(0)=h], h=-5..5, 1.0), seq([x(0)=5, y(0)=h], h=-5..5, 1.0);DEplot([de1, de2], [x(t), y(t)], t=0..2, [ic], x=-5..5, y=-5..5, stepsize=0.1, linecolor=blue, thickness=1, arrows=none); p5 := %:plot([[2*t,t,t=-5..5], [t,2*t, t=-5..5]], scaling=constrained, thickness=2, view=[-5..5, -5..5], color=red); p6 := %:display([p5,p6]);
<Text-field style="Heading 1" layout="Heading 1">Spiral</Text-field>L := < <-2, 3> | <-3, -2>>;P := <<1,1>|<-1,2>>;A := P . L . P^(-1);(e,v) := Eigenvectors(A);map(Re, Column(v,1));map(Im, Column(v,1));Q := < map(Im, Column(v,1)) | map(Re, Column(v,1)) >;Q^(-1) . A . Q;tmp := A . < x(t), y(t) >;de1 := diff(x(t),t) = tmp[1];de2 := diff(y(t),t) = tmp[2];ic := [x(0)=0, y(0)=1], [x(0)=0, y(0)=2];DEplot([de1, de2], [x(t), y(t)], t=-1..1, [ic], x=-5..5, y=-5..5, stepsize=0.1, linecolor=blue, thickness=1);ic := seq([x(0)=h, y(0)=0], h=-5..5, 0.5);DEplot([de1, de2], [x(t), y(t)], t=0..1, [ic], x=-5..5, y=-5..5, stepsize=0.1, linecolor=blue, thickness=1);DEplot([de1, de2], [x(t), y(t)], t=-2..2, [ic], x=-5..5, y=-5..5, stepsize=0.05, linecolor=blue, thickness=1, arrows=none); p7 := %:plot([[3*t,0*t,t=-1..1], [-1*t,5*t, t=-1..1]], scaling=constrained, thickness=2, color=red); p8 := %:display([p7,p8]);
<Text-field style="Heading 1" layout="Heading 1">Center</Text-field>We choose the eigenvalues to be pure imaginary, therefore we get a center.L := < <0, 1> | <-1, 0>>;P := <<1,1>|<-1,2>>;A := P . L . P^(-1);(e,v) := Eigenvectors(A);map(Re, Column(v,1));map(Im, Column(v,1));Q := < map(Im, Column(v,1)) | map(Re, Column(v,1)) >;Q^(-1) . A . Q;tmp := A . < x(t), y(t) >;de1 := diff(x(t),t) = tmp[1];de2 := diff(y(t),t) = tmp[2];ic := [x(0)=0, y(0)=1], [x(0)=0, y(0)=2];DEplot([de1, de2], [x(t), y(t)], t=-2..2, [ic], x=-5..5, y=-5..5, stepsize=0.1, linecolor=blue, thickness=1);ic := seq([x(0)=h, y(0)=0], h=-5..5, 0.5);DEplot([de1, de2], [x(t), y(t)], t=-Pi/2..Pi/2, [ic], x=-5..5, y=-5..5, stepsize=0.1, linecolor=blue, thickness=1, arrows=none); p7 := %:plot([[3*t,0*t,t=-1..1], [-1*t,5*t, t=-1..1]], thickness=4, color=red); p8 := %:display([p7,p8]);