Iterative solutions

1153 days ago by ncr006

Initializing $A$ and $\mathbf{b}$

A = matrix(RR, 3,3, [2,0,-3, -1,-1,-1,1,2,-4 ])*matrix(RR,3,3,[1/2,0,0, 0,1/3,0,0,0,-1/7 ])*matrix(RR, 3,3, [2,0,-3, -1,-1,-1,1,2,-4 ]).inverse() b = matrix(RR,3,1,[1,0,-1]) show(A) show(b) 
       

                                
                            

                                
A.eigenvalues() 
       
__main__:1: UserWarning: Using generic algorithm for an inexact ring,
which will probably give incorrect results due to numerical precision
issues.
[0.500000000000000, 0.333333333333333, -0.142857142857143]
__main__:1: UserWarning: Using generic algorithm for an inexact ring, which will probably give incorrect results due to numerical precision issues.
[0.500000000000000, 0.333333333333333, -0.142857142857143]
#Actual solution show(A.inverse()*b) 
       

                                
                            

                                

Richardson Iteration

This is done with the matrix $Q$ being the identity.

Q = identity_matrix(RR,3) xx0 = matrix(RR,3,1,[1,1,1]) # an initial estimate of the solution show(Q) show(xx0) 
       

                                
                            

                                
#Let's try with omega = 1 omega=1 show(Q) show(omega*A) show(Q-omega*A) 
       

                                
                            

                                
omega = 1 xx1 = (Q-omega*A)*xx0+omega*b show(xx1) 
       

                                
                            

                                
xx2 = (Q-omega*A)*xx1+omega*b show(xx2) 
       

                                
                            

                                
xx3 = (Q-omega*A)*xx2+omega*b show(xx3) 
       

                                
                            

                                
#Let's try with omega = 1 # what happens as I change the range? omega = 1 xx = matrix(RR,3,1,[1,1,1]) for i in range(200): xx = (Q-omega*A)*xx + omega*b show(xx) 
       

                                
                            

                                
#Let's try with omega = 2 omega = 2 xx = matrix(RR,3,1,[1,1,1]) for i in range(200): xx = (Q-omega*A)*xx + omega*b show(xx) 
       

                                
                            

                                
#Let's try with omega = 1/6 omega = 1/6 xx = matrix(RR,3,1,[1,1,1]) for i in range(1000): xx = (Q-omega*A)*xx + omega*b show(xx) 
       

                                
                            

                                

Jacobi iteration

This is where you choose $Q$ to be the diagonal of $A$.

# initialize new input B = matrix(RR,3,3,[1,2,3,-1,0,1,1,1,2]) A = B*diagonal_matrix(RR,[1/4,1/4,1/5])*B.inverse() show(A) b = matrix(RR,3,1,[1,0,-1]) show(b) show(A.inverse()*b) Q = matrix(AA,3,3,[A[0,0],0,0,0,A[1,1],0,0,0,A[2,2]]) show(Q) 
       
Traceback (click to the left of this block for traceback)
...
TypeError: Illegal initializer for algebraic number
Traceback (most recent call last):    show(A.inverse()*b)
  File "", line 1, in <module>
    
  File "/tmp/tmpWurWLo/___code___.py", line 9, in <module>
    Q = matrix(AA,_sage_const_3 ,_sage_const_3 ,[A[_sage_const_0 ,_sage_const_0 ],_sage_const_0 ,_sage_const_0 ,_sage_const_0 ,A[_sage_const_1 ,_sage_const_1 ],_sage_const_0 ,_sage_const_0 ,_sage_const_0 ,A[_sage_const_2 ,_sage_const_2 ]])
  File "/nfs/software/sage-6.3/local/lib/python2.7/site-packages/sage/matrix/constructor.py", line 729, in _matrix_constructor
    return matrix_space.MatrixSpace(ring, nrows, ncols, sparse=sparse)(entries)
  File "/nfs/software/sage-6.3/local/lib/python2.7/site-packages/sage/matrix/matrix_space.py", line 466, in __call__
    return self.matrix(entries, coerce, copy)
  File "/nfs/software/sage-6.3/local/lib/python2.7/site-packages/sage/matrix/matrix_space.py", line 1364, in matrix
    return MC(self, x, copy=copy, coerce=coerce)
  File "matrix_generic_dense.pyx", line 109, in sage.matrix.matrix_generic_dense.Matrix_generic_dense.__init__ (build/cythonized/sage/matrix/matrix_generic_dense.c:2856)
  File "parent.pyx", line 1089, in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:8902)
  File "coerce_maps.pyx", line 95, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4203)
  File "coerce_maps.pyx", line 90, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4110)
  File "/nfs/software/sage-6.3/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 693, in _element_constructor_
    return AlgebraicReal(x)
  File "/nfs/software/sage-6.3/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 4448, in __init__
    AlgebraicNumber_base.__init__(self, AA, x)
  File "/nfs/software/sage-6.3/local/lib/python2.7/site-packages/sage/rings/qqbar.py", line 3121, in __init__
    raise TypeError("Illegal initializer for algebraic number")
TypeError: Illegal initializer for algebraic number
#Actual solution show(A.inverse()*b) 
       

                                
                            

                                

Recall that the iterative step is of the form $Q\mathbf{x}^{(k+1)} = (Q-\omega A)\mathbf{x}^{(k)} + \omega\mathbf{b}.$  This means that each time we iterate I have to solve for $x^{(k+1)}$.  In the code below I'm cheating a little by just multiplying both sides by $Q^{-1}$.

# start to see convergence at 1000 iterations xx = matrix(RR,3,1,[1,1,1]) omega = 1/100 Qi = Q.inverse() for i in range(100): xx = Qi*((Q-omega*A)*xx + omega*b) show(xx) 
       

                                
                            

                                

Gauss-Seidel

For this algorithm we choose $Q$ to be the lower triangular part of $A$.

Q = matrix(RR,3,3,[A[0][0],0,0,A[1][0],A[1][1],0,A[2,0],A[2,1],A[2,2]]) show(A) show(Q) 
       

                                
                            

                                
# start to see convergence at 1000 iterations xx = matrix(RR,3,1,[1,1,1]) omega = 1/100 Qi = Q.inverse() for i in range(100): xx = Qi*((Q-omega*A)*xx + omega*b) show(xx)