Exercises on Systems of Differential Equations

$Post:=If[MatrixQ[#],MatrixForm[#],#]&

Exercise 1

M=Table[i+j,{i,4},{j,4}]

P=Transpose[Eigensystem[M][[2]]//Together]; P//TableForm

T=(Inverse[P]//Together).M.P//Together;
T//TableForm

z=Table[Z[i][x]/.DSolve[Z[i]'[x]==T[[i,i]]Z[i][x],Z[i][x],x] \
                                        /.{C[1]->K[i]},{i,4}]

y=P.z//Together//Simplify

y//N

Solution check:

D[y,x]-M.y//Expand

Exercise 2

M=Table[Mod[i+j,2],{i,4},{j,4}]

P=Transpose[Eigensystem[M][[2]]]

T=Inverse[P].M.P

z=Table[Z[i][x]/.DSolve[Z[i]'[x]==T[[i,i]]Z[i][x],Z[i][x],x] \
                                        /.{C[1]->K[i]},{i,4}]

y=P.z

D[y,x]-M.y//Expand

Exercise 3

M={{1,-2,-2,-2,-2},   \
   {2,6,4,4,4},       \
   {-3,-6,-6,-10,-10},\
   {1,2,3,6,1},       \
   {1,2,3,4,9}}

P=Transpose[Eigensystem[M][[2]]]

T=Inverse[P].M.P

z=Table[Z[i][x]/.DSolve[Z[i]'[x]==T[[i,i]]Z[i][x],Z[i][x],x] \
                                        /.{C[1]->K[i]},{i,5}]

y=P.z

D[y,x]-M.y//Expand

Exercise 4

M={{2,0,-1,-1,-1},   \
   {3,8,8,7,7},      \
   {-3,-6,-6,-7,-8}, \
   {1,2,3,5,4},      \
   {0,0,0,0,2}}

EV=Eigensystem[M][[2]]

Since we have only 3 linearly independent, nonzero eigenvectors,
M is not diagonalizable.

NullSpace[EV]

P=Transpose[Join[{EV[[1]]},{EV[[2]]},{EV[[4]]},%]]

T=Inverse[P].M.P

<<LinearAlgebra`Master`
B=SubMatrix[T,{4,4},{2,2}]

EVB=Eigensystem[B][[2]]

PB=Transpose[EVB]

Table[If[i<4 || j<4,IdentityMatrix[5][[i,j]],PB[[i-3,j-3]]], \
                                                {i,5},{j,5}]

P=P.%

T=Inverse[P].M.P

z=Array[0&,5];
Do[z[[6-i]]=Z[6-i][x]/.DSolve[Z[6-i]'[x]==T[[6-i,6-i]]Z[6-i][x] \
                            +Sum[T[[6-i,j]]z[[j,1]],{j,7-i,5}], \
                            Z[6-i][x],x]/.{C[1]->K[6-i]},{i,5}]
z//Expand

y=P.z//Expand

D[y,x]-M.y//Expand

Exercise 5

M={{-205,-409,-613,-817,-1021,-1225},   \
   {494,986,1479,1972,2465,2958},       \
   {-662,-1323,-1986,-2648,-3310,-3972},\
   {539,1078,1618,2156,2695,3234},      \
   {-272,-544,-816,-1087,-1360,-1632},  \
   {70,140,210,280,351,421}}

EV=Eigensystem[M][[2]]

NullSpace[EV]

P=Transpose[Join[{EV[[1]]},{EV[[3]]},{EV[[5]]},%]]

T=Inverse[P].M.P; T//TableForm

B=Array[T[[3+#1,3+#2]]&,{3,3}]

EVB=Eigensystem[B][[2]]

PB=Transpose[EVB]

Table[If[i<4 || j<4,IdentityMatrix[6][[i,j]],PB[[i-3,j-3]]], \
                                                {i,6},{j,6}];
%//TableForm

P=P.%%; P//TableForm

T=Inverse[P].M.P; T//TableForm

z=Array[0&,6];
Do[z[[7-i]]=Z[7-i][x]/.DSolve[Z[7-i]'[x]==T[[7-i,7-i]]Z[7-i][x] \
                            +Sum[T[[7-i,j]]z[[j,1]],{j,8-i,6}], \
                            Z[7-i][x],x]/.{C[1]->K[7-i]},{i,6}]
z//Expand

y=P.z//Expand

D[y,x]-M.y//Expand

Problems on Differential Equations

Problem 1

M={{0,1},  \
   {-2,-1}}

P=Transpose[Eigensystem[M][[2]]]

T=Inverse[P].M.P//Together

k={0,E^x}

IPk=Inverse[P].k//Together

z=Table[Z[i][x]/.DSolve[Z[i]'[x]==T[[i,i]]Z[i][x]+IPk[[i]], \
                        Z[i][x],x]/.{C[1]->K[i]},{i,2}]

u=P.z//Together//Simplify

yRe=Re[u[[1]]]

 Mathematica does not know if the constants are real

yRe=(E^x-E^(-x/2)Cos[Sqrt[7]x/2](K[1]-K[2])+       \
     Sqrt[7]E^(-x/2)Sin[Sqrt[7]x/2](K[1]+K[2]))/4

Solution check:

D[D[yRe,x],x]+D[yRe,x]+2yRe//Expand

Problem 2

M={{1/75,-1/150},     \
   {1/75,-1/300}}

P=Transpose[Eigensystem[M][[2]]]

T=Inverse[P].M.P//Together

Clear[K]
z=Table[Z[i][t]/.DSolve[Z[i]'[t]==T[[i,i]]Z[i][t],Z[i][t],t] \
                                        /.{C[1]->K[i]},{i,2}]

r=(P.z)[[1]]

f=(P.z)[[2]]

a)

Solve[(r/.{t->0})-1000==0,K[2]]//Expand;
Solve[(f/.{t->0}/.%)-1000==0,K[1]]//Expand;
Re[r/.{t->{10,20,100}}/.%%/.%//N]

Re[f/.{t->{10,20,100}}/.%%%/.%%//N]

b)

Solve[(r/.{t->0})-2000==0,K[2]]//Expand;
Solve[(f/.{t->0}/.%)-1000==0,K[1]]//Expand;
Re[r/.{t->{10,20,100}}/.%%/.%//N]

Re[f/.{t->{10,20,100}}/.%%%/.%%//N]

Problem 3

Since both w1/g and w2/g are equal to 1, we will just
skip writing them.

k1=12; k2=10; k3=12;
M={{0,1,0,0},      \
   {-k1-k2,0,k2,0},\
   {0,0,0,1},      \
   {k2,0,-k2-k3,0}}

P=Transpose[Eigensystem[M][[2]]]//Together

T=Inverse[P].M.P

Clear[K]
z=Table[Z[i][t]/.DSolve[Z[i]'[t]==T[[i,i]]Z[i][t],Z[i][t],t] \
                                        /.{C[1]->K[i]},{i,4}]

u=P.z

Clear[y]
Solve[(u[[1]]/.{t->0})-10==0,K[1]]//Expand;
Solve[(u[[2]]/.{t->0}/.%)==0,K[2]]//Expand;
Solve[(u[[3]]/.{t->0}/.%%/.%)-12==0,K[3]]//Expand;
Solve[(u[[4]]/.{t->0}/.%%%/.%%/.%)==0,K[4]]//Expand;
y[1]=u[[1]]/.%%%%/.%%%/.%%/.%//ComplexExpand

y[2]=u[[3]]/.%%%%%/.%%%%/.%%%/.%%//ComplexExpand

Exercises on the Matrix Exponential Function

Exercise 1

J3={{2,1,0},  \
    {0,2,1},  \
    {0,0,2}}

MatrixExp[J3 t]

Exercise 2

A=Array[1&,{4,4}]

MatrixExp[A t]

Exercise 3

A=Table[If[i>j,i,0],{i,4},{j,4}]

MatrixExp[A t]//TableForm

Exercise 4

J[n_]:=Table[Which[j>i+1,0,j==i+1,1,j==i,2,j<i,0],{i,n},{j,n}]

J[4]

MatrixExp[J[4]t]

J[5]

MatrixExp[J[5]t]

J[6]

MatrixExp[J[6]t]

Problems on the Matrix Exponential

Problem 1

A=Table[If[i<j,i,0],{i,4},{j,4}]

MatrixExp[A t]//Expand

Solution check:

D[%,t]-A.%

Actually, since every column of E^(At) is a solution of the
differential equation, we get a general solution from:

y=Transpose[{Sum[Transpose[%%][[i]]K[i],{i,4}]}]//Expand

D[y,t]-A.y//Expand

Problem 2

A={{0,1},   \
   {-13,6}}

MatrixExp[A t];
u=Transpose[{Array[Sum[%[[#,j]]K[j],{j,2}]&,2]}]//Expand

y=K[4]+Integrate[K[3]+Integrate[u[[1,1]]-4/13,t],t]//Expand

D[y,{t,4}]-6D[y,{t,3}]+13D[y,{t,2}]+4//Expand

Problem 3

M={{1/75,-1/150},     \
   {1/75,-1/300}}

A=MatrixExp[M t];
Clear[K]
y=Transpose[{Array[Sum[A[[#,j]]K[j],{j,2}]&,2]}]//Expand;
r=y[[1]]

f=y[[2]]

a)

LinearSolve[A/.{t->0},{1000,1000}];
K[1]=%[[1]]

K[2]=%%[[2]]

Re[r/.{t->{10,20,100}}]//N

Re[f/.{t->{10,20,100}}]//N

b)

LinearSolve[A/.{t->0},{2000,1000}];
K[1]=%[[1]]

K[2]=%%[[2]]

Re[r/.{t->{10,20,100}}]//N

Re[f/.{t->{10,20,100}}]//N

Problem 4

k1=12; k2=10; k3=12;
M={{0,1,0,0},      \
   {-k1-k2,0,k2,0},\
   {0,0,0,1},      \
   {k2,0,-k2-k3,0}}

Clear[K]
A=MatrixExp[M t];
u=Transpose[{Array[Sum[%[[#,j]]K[j],{j,4}]&,4]}]//Expand;
b={10,0,12,0}

LinearSolve[A/.{t->0},b];
Do[K[i]=%[[i]],{i,4}]
y={u[[1]],u[[3]]}//ComplexExpand

Problems on Linear Recurrence Systems

<<DiscreteMath`Master`

Problem 1

Unprotect[M];
M=1/600{{6,-3},   \
        {3,-2}}

P=Transpose[Eigensystem[M][[2]]]//Together

T=Inverse[P].M.P//Together

r=Table[R[i][n]/.RSolve[R[i][n+1]==(T[[i,i]]//N)R[i][n],R[i][n],n],{i,2}]

p=P.r

a)

Solve[(p[[1]]/.{n->0})-1000000.0==0,R[1][0]]//Expand;
Solve[(p[[2]]/.{n->0}/.%)-1000000.0==0,R[2][0]]//Expand;
p[[1]]/.{n->{1,2,3,4,5,6,7,8,9,10}}/.%%/.%//N//Expand

p[[2]]/.{n->{1,2,3,4,5,6,7,8,9,10}}/.%%%/.%%//N//Expand

b)

Looking at the expressions for p[[1]] and p[[2]] as functions
of n above, we see that both are exponentials of basis < 1 in n
and therefore the populations will go extinct quite rapidly,
no matter the initial numbers.

Problem 2

M={{1,1,0},   \
   {-1,2,1},  \
   {-1,-1,4}}

EV=Eigensystem[M][[2]]

NullSpace[EV]

P=Transpose[Join[{EV[[1]]},{EV[[3]]},%]]

T=Inverse[P].M.P//Together

Having now upper triangularized our system of recurrence 
equations, we "RSolve" it in a do-loop in the order 3,2,1.

r=Array[0&,3];
Do[r[[4-i]]=R[4-i][n]/.RSolve[R[4-i][n+1]==T[[4-i,4-i]]R[4-i][n] \
                             +Sum[T[[4-i,j]]r[[j,1]],{j,5-i,3}], \
                                             R[4-i][n],n],{i,3}]
r//Expand

p=P.r//Expand

Clear[R]; R[1][0]:=1; R[2][0]:=1; R[3][0]:=1;
p/.{n->0}

p/.{n->10}

Clear[R]; R[1][0]:=1; R[2][0]:=1; R[3][0]:=-1;
p/.{n->0}

p/.{n->10}

As one can guess from the above tests, the future is dependent 
on the initial populations. For instance, for p[0]={3,1,5}, the
populations would subsequenly grow indefinitely, whereis for
p[0]={1,5,3} they go extinct.

Exercises on Quadratics

Exercise 1

Clear[x]; 
X=Array[Subscripted[x[#]]&,2]

Clear[y]; 
Y=Array[Subscripted[y[#]]&,2]

a)

M={{3,-4},{-4,-3}}

X.M.X//Expand

S=Solve[%-1==0,X[[2]]]//Simplify

Clear[f]
f[i_][x1_]:=(X[[2]]/.S[[i]])/.{X[[1]]->x1}
P1=Plot[{f[1][x1],f[2][x1]},{x1,-2,-Sqrt[12.001]/10}]
P2=Plot[{f[1][x1],f[2][x1]},{x1,Sqrt[12.001]/10,2}]
Show[P1,P2]

P=Transpose[Eigensystem[M][[2]]]

T=Inverse[P].M.P

Print[Y.T.Y,"=1"]

b)

M={{0,1/2},{1/2,0}}

X.M.X

S=Solve[%-1==0,X[[2]]]

Plot[f[1][x1],{x1,-2,2}]

P=Transpose[Eigensystem[M][[2]]]

T=Inverse[P].M.P

Print[Y.T.Y,"=1"]

c)

M={{5,3},{3,3}}

X.M.X//Expand

S=Solve[%-8==0,X[[2]]]//Simplify

Plot[{f[1][x1],f[2][x1]},{x1,-2,2}]

P=Transpose[Eigensystem[M][[2]]]

T=Inverse[P].M.P//Expand

Print[Y.T.Y,"=8"]

d)

M={{-1,0},    \
   {0,1}}

k={-2,0};
X.M.X+X.k

Solve[%-1==0,X[[2]]]//Simplify

Plot[{1+x1,-1-x1},{x1,-3,1}]

In order to write this one in the standard form, we need
the translation Y=X+{1,0}:

X=Y-{1,0}

Print[(X.M.X+X.k-1)//Simplify,"=",1-1]

Exercise 2

X=Array[Subscripted[x[#]]&,3]
Y=Array[Subscripted[y[#]]&,3]

a)

M=Table[If[i==j,1,-1],{i,3},{j,3}]

X.M.X//Expand

S=Solve[%-3==0,X[[3]]]//Simplify

f[i_][x1_,x2_]:=(X[[3]]/.S[[i]])/.{X[[1]]->x1}/.{X[[2]]->x2}
P1=Plot3D[f[1][x1,x2],{x1,-1/2,1/2},{x2,-3/2,3/2}]
P2=Plot3D[f[2][x1,x2],{x1,-1/2,1/2},{x2,-3/2,3/2}]
Show[P1,P2]

P=Transpose[Eigensystem[M][[2]]]

T=Inverse[P].M.P

Print[Y.T.Y,"=3"]

b)

M={{23,0,-36},   \
   {0,50,0},     \
   {-36,0,2}}

X.M.X//Expand

S=Solve[%-50==0,X[[3]]]//Simplify

P1=Plot3D[f[1][x1,x2],{x1,-0.1,0.1},{x2,-1,1}]
P2=Plot3D[f[2][x1,x2],{x1,-0.1,0.1},{x2,-1,1}]
Show[P1,P2]

P=Transpose[Eigensystem[M][[2]]]

T=Inverse[P].M.P

Print[Y.T.Y,"=50"]

c)

M={{68,0,24},   \
   {0,50,0},    \
   {24,0,82}}

X.M.X//Expand

S=Solve[%-500==0,X[[3]]]//Simplify

P1=Plot3D[f[1][x1,x2],{x1,-2,2},{x2,-2,2}]
P2=Plot3D[f[2][x1,x2],{x1,-2,2},{x2,-2,2}]
Show[P1,P2]

P=Transpose[Eigensystem[M][[2]]]

T=Inverse[P].M.P

Print[Y.T.Y,"=500"]

d)

M={{9,-2,0},   \
   {-2,6,0},   \
   {0,0,5}}

X.M.X//Expand

S=Solve[%-500==0,X[[3]]]//Simplify

P1=Plot3D[f[1][x1,x2],{x1,-3,3},{x2,-3,3}]
P2=Plot3D[f[2][x1,x2],{x1,-3,3},{x2,-3,3}]
Show[P1,P2]

P=Transpose[Eigensystem[M][[2]]]

T=Inverse[P].M.P

Print[Y.T.Y,"=500"]