# A matrix with complex eigenvalues
# Using Octave
A = [1 -1;1 0]
eig(A)
# Create an initial starting point and iterate
x1(1) = 6
x2(1) = 4
for k=1:10
x = [x1(k);x2(k)];
xnew = A*x;
x1(k+1) = xnew(1);
x2(k+1) = xnew(2);
x = xnew;
end
plot(x1,x2,'+')
A =
1 -1
1 0
ans =
(0.5,0.866025)
(0.5,-0.866025)
x1 = 6
x2 = 4
6
|----------------------------------------------------------------------|
| + + + + +
|
|
|
4 |-+
+-|
|
|
|
|
|
|
2 |-+ A
+-|
|
|
|
|
0 |-+
+-|
|
|
|
|
-2 |-+ A
+-|
|
|
|
|
|
|
-4 |-+
+-|
|
|
| + + + + +
|
-6
|----------------------------------------------------------------------|
-6 -4 -2 0 2 4
6
A =
1 -1
1 0
ans =
(0.5,0.866025)
(0.5,-0.866025)
x1 = 6
x2 = 4
6 |----------------------------------------------------------------------|
| + + + + + |
| |
4 |-+ +-|
| |
| |
| |
2 |-+ A +-|
| |
| |
0 |-+ +-|
| |
| |
-2 |-+ A +-|
| |
| |
| |
-4 |-+ +-|
| |
| + + + + + |
-6 |----------------------------------------------------------------------|
-6 -4 -2 0 2 4 6
|
A = [1 -0.95;1 0]
eig(A)
x1(1) = 6
x2(1) = 4
for k=1:2000
x = [x1(k);x2(k)];
xnew = A*x;
x1(k+1) = xnew(1);
x2(k+1) = xnew(2);
x = xnew;
end
figure(1)
plot(x1,'b')
plot(x2,'r')
plot(x1,x2,'x')
A =
1 -0.95
1 0
ans =
(0.5,0.83666)
(0.5,-0.83666)
x1 =
6 2 -4 -6 -2 4 6 2 -4 -6 -2
x2 =
4 6 2 -4 -6 -2 4 6 2 -4 -6
6
|----------------------------------------------------------------------|
| + + + + B +
|
|
|
4 |-+ B
+-|
| B
|
| B B
|
| B B B B B BB B
|
2 |-+ B BBB B B
+-|
| BB BBB B
|
| B BBBB BB BBB B B
|
0 |-+ B B BB BBBBBBBB B BB B
+-|
| B B BBB BB BBBBB
|
| B BB BB
|
-2 |-+ B BB BB B B
+-|
| B B B B B B
|
| B
|
| B B
|
-4 |-+ B
+-|
| B
|
| + B + + + +
|
-6
|----------------------------------------------------------------------|
-6 -4 -2 0 2 4
6
A =
1 -0.95
1 0
ans =
(0.5,0.83666)
(0.5,-0.83666)
x1 =
6 2 -4 -6 -2 4 6 2 -4 -6 -2
x2 =
4 6 2 -4 -6 -2 4 6 2 -4 -6
6 |----------------------------------------------------------------------|
| + + + + B + |
| |
4 |-+ B +-|
| B |
| B B |
| B B B B B BB B |
2 |-+ B BBB B B +-|
| BB BBB B |
| B BBBB BB BBB B B |
0 |-+ B B BB BBBBBBBB B BB B +-|
| B B BBB BB BBBBB |
| B BB BB |
-2 |-+ B BB BB B B +-|
| B B B B B B |
| B |
| B B |
-4 |-+ B +-|
| B |
| + B + + + + |
-6 |----------------------------------------------------------------------|
-6 -4 -2 0 2 4 6
|
# let's try it in Sage to get a better picture
A = matrix(QQ,[[1, -0.95],[1, 0]])
pts = []
npts = 20
x = matrix([[6],[4]])
for k in range(npts):
xnew = A*x
pts.append((xnew[0][0],xnew[1][0]))
x = xnew
end
for k in range(npts):
print pts[k]
points(pts)
(11/5, 6) (-7/2, 11/5) (-559/100, -7/2) (-453/200, -559/100) (6091/2000, -453/200) (20789/4000, 6091/2000) (92161/40000, 20789/4000) (-210669/80000, 92161/40000) (-3857749/800000, -210669/80000) (-3712787/1600000, -3857749/800000) (36169361/16000000, -3712787/1600000) (5715267/1280000, 36169361/16000000) (741598891/320000000, 5715267/1280000) (-1231554043/640000000, 741598891/320000000) (-26405919359/6400000000, -1231554043/640000000) (-29412311901/12800000000, -26405919359/6400000000) (207589348811/128000000000, -29412311901/12800000000) (974012623741/256000000000, 207589348811/128000000000) (5795928610001/2560000000000, 974012623741/256000000000) (-6914382631077/5120000000000, 5795928610001/2560000000000) (11/5, 6) (-7/2, 11/5) (-559/100, -7/2) (-453/200, -559/100) (6091/2000, -453/200) (20789/4000, 6091/2000) (92161/40000, 20789/4000) (-210669/80000, 92161/40000) (-3857749/800000, -210669/80000) (-3712787/1600000, -3857749/800000) (36169361/16000000, -3712787/1600000) (5715267/1280000, 36169361/16000000) (741598891/320000000, 5715267/1280000) (-1231554043/640000000, 741598891/320000000) (-26405919359/6400000000, -1231554043/640000000) (-29412311901/12800000000, -26405919359/6400000000) (207589348811/128000000000, -29412311901/12800000000) (974012623741/256000000000, 207589348811/128000000000) (5795928610001/2560000000000, 974012623741/256000000000) (-6914382631077/5120000000000, 5795928610001/2560000000000) 
|
You can also convert scalar nth order difference equations into a nxn matrix first order difference form.
# Solve the difference equation
# y0 = 1
# y1 = 1/2
# ynp1 = 5/6 * yn - 1/6 * ynm1
format short
inits = [1;1/2]
# or use [1;1/3]
A = [ 0 1;-1/6 5/6]
w = inits
for k=1:10
w(:,k+1) = A*w(:,k)
end
plot(w(2,:))
WARNING: Output truncated!full_output.txt
inits =
1.00000
0.50000
A =
0.00000 1.00000
-0.16667 0.83333
w =
1.00000
0.50000
w =
1.00000 0.50000
0.50000 0.25000
w =
1.00000 0.50000 0.25000
0.50000 0.25000 0.12500
w =
1.000000 0.500000 0.250000 0.125000
0.500000 0.250000 0.125000 0.062500
w =
1.000000 0.500000 0.250000 0.125000 0.062500
0.500000 0.250000 0.125000 0.062500 0.031250
w =
1.000000 0.500000 0.250000 0.125000 0.062500 0.031250
0.500000 0.250000 0.125000 0.062500 0.031250 0.015625
w =
Columns 1 through 6:
1.0000000 0.5000000 0.2500000 0.1250000 0.0625000 0.0312500
0.5000000 0.2500000 0.1250000 0.0625000 0.0312500 0.0156250
Column 7:
0.0156250
0.0078125
w =
Columns 1 through 6:
1.0000000 0.5000000 0.2500000 0.1250000 0.0625000 0.0312500
0.5000000 0.2500000 0.1250000 0.0625000 0.0312500 0.0156250
...
w =
Columns 1 through 6:
1.0000000 0.5000000 0.2500000 0.1250000 0.0625000 0.0312500
0.5000000 0.2500000 0.1250000 0.0625000 0.0312500 0.0156250
Columns 7 through 9:
0.0156250 0.0078125 0.0039062
0.0078125 0.0039062 0.0019531
w =
Columns 1 through 6:
1.00000000 0.50000000 0.25000000 0.12500000 0.06250000
0.03125000
0.50000000 0.25000000 0.12500000 0.06250000 0.03125000
0.01562500
Columns 7 through 10:
0.01562500 0.00781250 0.00390625 0.00195312
0.00781250 0.00390625 0.00195312 0.00097656
w =
Columns 1 through 6:
1.00000000 0.50000000 0.25000000 0.12500000 0.06250000
0.03125000
0.50000000 0.25000000 0.12500000 0.06250000 0.03125000
0.01562500
Columns 7 through 11:
0.01562500 0.00781250 0.00390625 0.00195312 0.00097656
0.00781250 0.00390625 0.00195312 0.00097656 0.00048828
0.5
|---------------------------------------------------------------------|
| * + + + + +
|
| *
|
| *
|
0.4 |-+ *
+-|
| *
|
| *
|
| *
|
0.3 |-+ *
+-|
| *
|
| **
|
| *
|
0.2 |-+ **
+-|
| *
|
|
|
| ****
|
0.1 |-+ *
+-|
| *****
|
| *
|
| + + ************ + +
|
0
|---------------------------------------------------------------------|
0 2 4 6 8 10
12
WARNING: Output truncated!full_output.txt
inits =
1.00000
0.50000
A =
0.00000 1.00000
-0.16667 0.83333
w =
1.00000
0.50000
w =
1.00000 0.50000
0.50000 0.25000
w =
1.00000 0.50000 0.25000
0.50000 0.25000 0.12500
w =
1.000000 0.500000 0.250000 0.125000
0.500000 0.250000 0.125000 0.062500
w =
1.000000 0.500000 0.250000 0.125000 0.062500
0.500000 0.250000 0.125000 0.062500 0.031250
w =
1.000000 0.500000 0.250000 0.125000 0.062500 0.031250
0.500000 0.250000 0.125000 0.062500 0.031250 0.015625
w =
Columns 1 through 6:
1.0000000 0.5000000 0.2500000 0.1250000 0.0625000 0.0312500
0.5000000 0.2500000 0.1250000 0.0625000 0.0312500 0.0156250
Column 7:
0.0156250
0.0078125
w =
Columns 1 through 6:
1.0000000 0.5000000 0.2500000 0.1250000 0.0625000 0.0312500
0.5000000 0.2500000 0.1250000 0.0625000 0.0312500 0.0156250
...
w =
Columns 1 through 6:
1.0000000 0.5000000 0.2500000 0.1250000 0.0625000 0.0312500
0.5000000 0.2500000 0.1250000 0.0625000 0.0312500 0.0156250
Columns 7 through 9:
0.0156250 0.0078125 0.0039062
0.0078125 0.0039062 0.0019531
w =
Columns 1 through 6:
1.00000000 0.50000000 0.25000000 0.12500000 0.06250000 0.03125000
0.50000000 0.25000000 0.12500000 0.06250000 0.03125000 0.01562500
Columns 7 through 10:
0.01562500 0.00781250 0.00390625 0.00195312
0.00781250 0.00390625 0.00195312 0.00097656
w =
Columns 1 through 6:
1.00000000 0.50000000 0.25000000 0.12500000 0.06250000 0.03125000
0.50000000 0.25000000 0.12500000 0.06250000 0.03125000 0.01562500
Columns 7 through 11:
0.01562500 0.00781250 0.00390625 0.00195312 0.00097656
0.00781250 0.00390625 0.00195312 0.00097656 0.00048828
0.5 |---------------------------------------------------------------------|
| * + + + + + |
| * |
| * |
0.4 |-+ * +-|
| * |
| * |
| * |
0.3 |-+ * +-|
| * |
| ** |
| * |
0.2 |-+ ** +-|
| * |
| |
| **** |
0.1 |-+ * +-|
| ***** |
| * |
| + + ************ + + |
0 |---------------------------------------------------------------------|
0 2 4 6 8 10 12
|
x0 = [1;-5;25]
A = [ 0 1 0;0 0 1;-40 22 1]
x1 = A*x0
x2 = A*x1
x3 = A*x2
x4 = A*x3
x5 = A*x4
x0 =
1
-5
25
A =
0 1 0
0 0 1
-40 22 1
x1 =
-5
25
-125
x2 =
25
-125
625
x3 =
-125
625
-3125
x4 =
625
-3125
15625
x5 =
-3125
15625
-78125
x0 =
1
-5
25
A =
0 1 0
0 0 1
-40 22 1
x1 =
-5
25
-125
x2 =
25
-125
625
x3 =
-125
625
-3125
x4 =
625
-3125
15625
x5 =
-3125
15625
-78125
|
|
|