2013-Example.sagews

a=matrix(SR,[[1,0,2],[0,1,-1],[-1,1,1]])

a

[ 1 0 2] [ 0 1 -1] [-1 1 1]

ev=a.eigenvalues()

[abs(i).n() for i in ev]

[2.00000000000000, 2.00000000000000, 1.00000000000000]

The spectral radius of A is 2.

A=matrix(QQ,[[1/2,0],[1/4,1/2]])

A^2

[1/4 0] [1/4 1/4]

A^3

[ 1/8 0] [3/16 1/8]

[A^k for k in range(10)]

[ [1 0] [1/2 0] [1/4 0] [ 1/8 0] [1/16 0] [1/32 0] [0 1], [1/4 1/2], [1/4 1/4], [3/16 1/8], [ 1/8 1/16], [5/64 1/32], [1/64 0] [1/128 0] [1/256 0] [ 1/512 0] [3/64 1/64], [7/256 1/128], [ 1/64 1/256], [9/1024 1/512] ]

(A^1000).n()

[9.33263618503219e-302 0.000000000000000] [4.66631809251609e-299 9.33263618503219e-302]

[abs(e) for e in A.eigenvalues()]

[1/2, 1/2]

We say that $A$ is convergent.

Example

D=diagonal_matrix(RDF,[3,4,5])
L=-matrix(RDF,[[0,0,0],[-2,0,0],[-1,2,0]])
U=-matrix(RDF,[[0,1,-1],[0,0,1],[0,0,0]])
Tj=D^(-1)*(L+U)
Tg=(D-L)^(-1)*U
omega=1.25
Tomega=(D-omega*L)^(-1)*((1-omega)*D+omega*U)
print "T (Jacobi):", [abs(e) for e in Tj.eigenvalues()]
print "T (Gauss-Seidel):", [abs(e) for e in Tg.eigenvalues()]
print "T (SOR (omega=1.25):", [abs(e) for e in Tomega.eigenvalues()]

T (Jacobi): [0.368403149864, 0.368403149864, 0.368403149864] T (Gauss-Seidel): [0.0, 0.166666666667, 0.1] T (SOR (omega=1.25): [0.68800410828, 0.150700431987, 0.150700431987]

# Jacobi iteration
x=vector(QQ,[0,0,0])
iterations=10
jacobi=[list(x)]
for i in range(iterations):
    xold=copy(x)
    x[0]=(4-xold[1]+xold[2])/3
    x[1]=(1-2*xold[0]-xold[2])/4
    x[2]=(1+xold[0]-2*xold[1])/5
    jacobi+=[list(x)]
    print x.n(), (x-xold).norm(Infinity).n()

(1.33333333333333, 0.250000000000000, 0.200000000000000) 1.33333333333333 (1.31666666666667, -0.466666666666667, 0.366666666666667) 0.716666666666667 (1.61111111111111, -0.500000000000000, 0.650000000000000) 0.294444444444444 (1.71666666666667, -0.718055555555556, 0.722222222222222) 0.218055555555556 (1.81342592592593, -0.788888888888889, 0.830555555555556) 0.108333333333333 (1.87314814814815, -0.864351851851852, 0.878240740740741) 0.0754629629629630 (1.91419753086420, -0.906134259259259, 0.920370370370370) 0.0421296296296296 (1.94216820987654, -0.937191358024691, 0.945293209876543) 0.0310570987654321 (1.96082818930041, -0.957407407407407, 0.963310185185185) 0.0202160493827160 (1.97357253086420, -0.971241640946502, 0.975128600823045) 0.0138342335390947

# Gauss-Seidel iteration
x=vector(QQ,[0,0,0])
iterations=10
gs=[list(x)]
for i in range(iterations):
    xold=copy(x)
    x[0]=(4-xold[1]+xold[2])/3
    x[1]=(1-2*x[0]-xold[2])/4
    x[2]=(1+x[0]-2*x[1])/5
    gs+=[list(x)]
    print x.n(), (x-xold).norm(Infinity).n()

(1.33333333333333, -0.416666666666667, 0.633333333333333) 1.33333333333333 (1.68333333333333, -0.750000000000000, 0.836666666666667) 0.350000000000000 (1.86222222222222, -0.890277777777778, 0.928555555555556) 0.178888888888889 (1.93961111111111, -0.951944444444444, 0.968700000000000) 0.0773888888888889 (1.97354814814815, -0.978949074074074, 0.986289259259259) 0.0339370370370370 (1.98841277777778, -0.990778703703704, 0.993994037037037) 0.0148646296296296 (1.99492424691358, -0.995960632716049, 0.997369102469136) 0.00651146913580247 (1.99777657839506, -0.998230564814815, 0.998847541604938) 0.00285233148148148 (1.99902603547325, -0.999224903137860, 0.999495168349794) 0.00124945707818930 (1.99957335716255, -0.999660470668724, 0.999778859700000) 0.000547321689300411

# SOR iteration
x=vector(QQ,[0,0,0])
omega=1.25
iterations=10
sor=[list(x)]
for i in range(iterations):
    xold=copy(x)
    x[0]=(1-omega)*xold[0]+omega*(4-xold[1]+xold[2])/3
    x[1]=(1-omega)*xold[1]+omega*(1-2*x[0]-xold[2])/4
    x[2]=(1-omega)*xold[2]+omega*(1+x[0]-2*x[1])/5
    sor+=[list(x)]
    print x.n(), (x-xold).norm(Infinity).n()

(1.66666666666667, -0.729166666666667, 1.03125000000000) 1.66666666666667 (1.98350694444444, -1.06716579861111, 1.02164713541667) 0.337999131944444 (2.04112865306713, -1.01567868833189, 1.01270972357856) 0.0576217086226852 (2.00154634169590, -1.00101858009527, 0.997718444576970) 0.0395823113712268 (1.99908717485612, -0.998461853191563, 0.999573109165570) 0.00255672690370432 (1.99940944060144, -0.999882033692251, 0.999900099705093) 0.00142018050068742 (2.00005686209853, -1.00003381154636, 1.00005609637154) 0.000647421497091728 (2.00002324610783, -1.00002360604691, 1.00000359045753) 0.0000525059140156560 (2.00000552034989, -0.999998670724931, 0.999999817835557) 0.0000249353219759930 (1.99999799014606, -0.999999019233669, 0.999999052694461) 7.53020382627932e-6

line(jacobi, corner_cutoff=0.99)+points(jacobi,)+line(gs,color='green', corner_cutoff=0.99)+points(gs,color='green')+line(sor,color='purple', corner_cutoff=0.99)+points(sor,color='purple')