restart:with(LinearAlgebra):with(linalg):
n:=20:
A:=Matrix(1..n,1..n,0): Di:=Matrix(1..n,1..n,0): U:=Matrix(1..n,1..n,0): L:=Matrix(1..n,1..n,0): b:=Vector(1..n, 0):
Digits:=5:
b[1]:=3: b[n]:=3:for i from 1 by 1 while i <= n do
A[i, i]:=4;
Di[i, i]:=4;
if (i < n) then A[i, i+1]:=-1; U[i, i+1]:=-1; end if:
if (i < n) then A[i+1, i]:=-1; L[i+1, i]:=-1; end if:
end do:SOR:= proc( A , B , eps, ome)
local Tj,maxi,X,Xold,i,j,k,m,n;
n,m := Dimension(A) :
Tj:=Matrix(1..n,1..n,0):
for j from 1 by 1 while j <= n do
Tj[j,j]:= A[j,j] :
B[j]:=B[j]/A[j,j] :
end do :
for i from 1 by 1 while i <= n do # Matrix mit 0 in der Diagonale
for k from 1 by 1 while k <= n do
Tj[i,k]:= (Tj[i,k] - A[i,k]) / A[i,i] :
end do :
end do :
maxi:=1000: X:=Vector(1..n,0) : Xold:=Vector(1..n,0) :
for j from 1 by 1 while (j <= 25 and maxi>=eps) do
for i from 1 by 1 while i<=n do
X[i]:= B[i] :
for k from 1 by 1 while k<=n do
X[i]:= X[i] + Tj[i,k]*X[k] : # Matrix * Vektor
end do:
X[i]:=(1 - ome)*Xold[i]+ome*X[i]; # zus\344tzl. zu Seidel
end do:
print( j , evalf(Transpose(X)) ) ;
maxi:=0 :
for m from 1 by 1 while m <= n do # Vektorenelemente der beiden Vektoren werden miteinander vergleichen
maxi:= max(abs(X[m]-Xold[m]), maxi ) :
Xold[m]:= X[m] :
end do :
end do :
if (j<=25) then " SOR converged " else " SOR NOT converged ! " end if ;
end proc:SOR(A, b, 10^(-5), 1):;# Omega BerechnungTg:=Multiply(Transpose(Di-L), U):
spectral:=max(eigenvalues(Tg));
myOmega:=evalf(2/(1+sqrt(1-spectral^2))); # Optimales omega:=1;# Alle omega zwischen 0 und 2 sind ok weil (Numerik S. 276):
# - A ist invertierbar (nicht singul\344r)\243
# - A[i, i] > 0 f\374r alle i = 1..n
# - max[1 <= k, j <=n] |A[k, j] | <= max[1 <= i <= n] | A[j, j] |. z. dt: das abs(gr\366sste) Element der Matrix liegt auf der Diagonalen
# - (A[i, j]^2) < A[i, i]*A[j, j], f\374r alle i != jTTdSMApJNFJUQUJMRV9TQVZFLzEwMTQ2MzZYKiUpYW55dGhpbmdHNiI2IltnbCEkJSEhISI1IjUkIiYrXSghIiYkIiZdKD1GKSQiJnZvJSEiJyQiCiY+PCJGLiQiJihISCEiKCQiJlVLKCEiKSQiJjYkPUY2JCImd2QlISIqJCImVzkiRjskIiY1J0chIzUkIiZFOighIzYkIiYiKXkiRkMkIiYKLlolISM3JCImdzYiRkgkIiZTeiMhIzgkIiZcKXAhIzkkIiZpdSJGUCQiJmNPJSEjOiQiJjk0IkZVRidGJgo=TTdSMApJNFJUQUJMRV9TQVZFLzEwMTUzNTZYKiUpYW55dGhpbmdHNiI2IltnbCEkJSEhISI1IjUkIiYpb3ohIiYkIiYlNEBGKSQiJmtjJiEiJyQiCiZbWSJGLiQiJl8lUSEiKCQiJnIrIkYzJCImQGojISIpJCImbCdvISIqJCImIil5IkY7JCImI1xZISM1JCImcT8iRkAkIiYjSEohIzYkIiYKRDUpISM3JCImYjQjRkgkIiZMVCYhIzgkIiZxUiJGTSQiJjtnJCEjOSQiJm9GKiEjOiQiJl0oPUYpRidGJgo=TTdSMApJNFJUQUJMRV9TQVZFLzEwMTYyNzZYKiUpYW55dGhpbmdHNiI2IltnbCEkJSEhISI1IjUkIiZ0LSkhIiYkIiZnOSNGKSQiJjd0JiEiJyQiCiYqRzpGLiQiJlQyJSEiKCQiJlYzIkYzJCImRClHISIpJCImS2woISIqJCImJkg/RjskIiZjUCYhIzUkIiZAVSJGQCQiJnp2JCEjNiQiJicKPSoqISM3JCImXWgjRkgkIiZuKW8hIzgkIiY8Ij1GTSQiJjd3JSEjOSQiJnZvJUYuJCImJTRARilGJ0YmCg==TTdSMApJNVJUQUJMRV9TQVZFLzIzMTY3MDU2WColKWFueXRoaW5nRzYiNiJbZ2whJCUhISEiNSI1JCImbC4pISImJCImQzojRikkIiZLdyYhIickCiImRmEiRi4kIiZ4NyUhIigkIiZTNSJGMyQiJjgmSCEiKSQiJmQpeSEiKiQiJmU1I0Y7JCImLGkmISM1JCImISpcIkZAJCImYSpSISM2JCIKJlUxIkZFJCImRiRHISM3JCImWmAoISM4JCImRisjRk0kIiY+PCJGLiQiJmtjJkYuJCImZzkjRilGJ0YmCg==TTdSMApJNVJUQUJMRV9TQVZFLzIzMTY4NDE2WColKWFueXRoaW5nRzYiNiJbZ2whJCUhISEiNSI1JCImIlEhKSEiJiQiJk86I0YpJCImKHBkISInCiQiJmNhIkYuJCImKzklISIoJCImKTM2RjMkIiYicEghIikkIiYkXHohIiokIiZ5NyNGOyQiJlZwJiEjNSQiJk5fIkZAJCImWjIlISM2JCIKJiYqMyJGRSQiJkAiSCEjNyQiJjR5KCEjOCQiJihISEYzJCImW1kiRi4kIiY3dCZGLiQiJkM6I0YpRidGJgo=TTdSMApJNVJUQUJMRV9TQVZFLzIzMTY5OTM2WColKWFueXRoaW5nRzYiNiJbZ2whJCUhISEiNSI1JCImJVEhKSEiJiQiJlE6I0YpJCImNXgmISInCiQiJmphIkYuJCImRzklISIoJCImKjQ2RjMkIiZPKEghIikkIiZmJ3ohIiokIiZROCNGOyQiJmFyJiEjNSQiJjJgIkZAJCImIyo0JSEjNiQKIiZ3NCJGRSQiJiZRSCEjNyQiJlZLKEY4JCImXyVRRjMkIiYqRzpGLiQiJkt3JkYuJCImTzojRilGJ0YmCg==TTdSMApJNVJUQUJMRV9TQVZFLzI1ODAzOTEyWColKWFueXRoaW5nRzYiNiJbZ2whJCUhISEiNSI1JCImJlEhKSEiJiQiJlI6I0YpJCImOHgmISInCiQiJmthIkYuJCImTjklISIoJCImLTYiRjMkIiZaKEghIikkIiYsKHohIiokIiZhOCNGOyQiJjdzJiEjNSQiJkdgIkZAJCImazUlISM2JCIKJiw1IkZFJCImOCQ9RjgkIiZyKyJGMyQiJlQyJUYzJCImRmEiRi4kIiYocGRGLiQiJlE6I0YpRidGJgo=TTdSMApJNVJUQUJMRV9TQVZFLzI1ODA1NjcyWColKWFueXRoaW5nRzYiNiJbZ2whJCUhISEiNSI1JCImJlEhKSEiJiQiJlI6I0YpJCImOHgmISInCiQiJmthIkYuJCImTzklISIoJCImLjYiRjMkIiZcKEghIikkIiY3KHohIiokIiZlOCNGOyQiJkdzJiEjNSQiJkxgIkZAJCImJTNUISM2JCIKJicpZSVGOyQiJkNqI0Y4JCImVjMiRjMkIiZ4NyVGMyQiJmNhIkYuJCImNXgmRi5GKkYnRiYKTTdSMApJNVJUQUJMRV9TQVZFLzI0ODk3MjgwWColKWFueXRoaW5nRzYiNiJbZ2whJCUhISEiNSI1JCImJlEhKSEiJiQiJlI6I0YpJCImOXgmISInCiQiJmthIkYuJCImTzklISIoJCImLjYiRjMkIiZdKEghIikkIiY5KHohIiokIiZmOCNGOyQiJktzJiEjNSQiJk5gIkZAJCImYj0iRjskIiYKdihvRjskIiZHKUdGOCQiJlM1IkYzJCImKzklRjMkIiZqYSJGLiQiJjh4JkYuRipGJ0YmCg==TTdSMApJNVJUQUJMRV9TQVZFLzIxMTYyMDgwWColKWFueXRoaW5nRzYiNiJbZ2whJCUhISEiNSI1JCImJlEhKSEiJiQiJlI6I0YpJCImOXgmISInCiQiJmthIkYuJCImUDklISIoJCImLjYiRjMkIiZdKEghIikkIiY6KHohIiokIiZnOCNGOyQiJkxzJiEjNSQiJllSJUZAJCImI0g9RjskIiYKVW0oRjskIiY7JkhGOCQiJikzNkYzJCImRzklRjNGLyQiJjh4JkYuRipGJ0YmCg==