Решаване на системи от линейни алгебрични уравнения (СЛАУ) - 
A.x = b, x = ?

Общ преглед

Пример

разгънат вид:
123 (Subscript[x, 1]) + 0.2 (Subscript[x, 2])   = 34
-(Subscript[x, 1]) + 34 (Subscript[x, 2]) + (Subscript[x, 3]) = -21
	0.9 (Subscript[x, 2]) -211(Subscript[x, 3]) = -0.7
(Subscript[x, 1]) =?,	(Subscript[x, 2]) =?,(Subscript[x, 3]) =?
матричен вид:
A.x = b, x = ?
където
А = (({{123, 0.2, 0}, {-1, 34, 1}, {0, 0.9, -211}})), x  = (({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}})), b = (({{34}, {-21}, {-0.7}}))

Обща постановка на задачата

разгънат вид:
(Subscript[a, 11])(Subscript[x, 1])+ (Subscript[a, 12])(Subscript[x, 2])+ ... (Subscript[a, 1n])(Subscript[x, n]) = (Subscript[b, 1])
(Subscript[a, 21])(Subscript[x, 1])+ (Subscript[a, 22])(Subscript[x, 2])+ ... (Subscript[a, 2n])(Subscript[x, n]) = (Subscript[b, 2])
.............................................
(Subscript[a, n1])(Subscript[x, 1])+ (Subscript[a, n2])(Subscript[x, 2])+ ... (Subscript[a, nn])(Subscript[x, n]) = (Subscript[b, n])
(Subscript[x, 1]) =?,	(Subscript[x, 2]) =?,...,(Subscript[x, n]) =?
матричен вид:
A.x = b, x = ?
където
A = (({{(Subscript[a, 11]), (Subscript[a, 12]), ..., (Subscript[a, 1n])}, {(Subscript[a, 21]), (Subscript[a, 22]), ..., (Subscript[a, 2n])}, {.., .., .., ..}, {(Subscript[a, n1]), (Subscript[a, n2]), ..., (Subscript[a, nn])}})), x = (({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {.}, {(Subscript[x, n])}})), b = (({{(Subscript[b, 1])}, {(Subscript[b, 2])}, {.}, {(Subscript[b, n])}}))

Възможности на Wolfram Mathematica

A = (({{123, 0.2, 0}, {-1, 34, 1}, {0, 0.9, -211}}))
(*LinearSolve*)

A

A//MatrixForm

A[[1,2]]

A[[1,2]]

A[[1]]

2*A[[1]]

решаване на СЛАУ с вграден метод във Wolfram:

A = (({{123, 0.2, 0}, {-1, 34, 1}, {0, 0.9, -211}})); b = {34,-21,-0.7};
LinearSolve[A,b]

Общa идея на итерационните методи за решаване на СЛАУ

A.x = b 
<=> конструиране на метода
x = C.x + d 
=> преминаване към итерационен процес
 (x^((k+1))) = C.(x^((k))) + d, k=0,1,2,....

Метод на Якоби (проста итерация)  за решаване на СЛАУ

Конструиране на метода

получаване на матрицата C и вектора d

за матрицата C:
(Subscript[c, ij]) = ((-(Subscript[a, ij]))/(Subscript[a, ii])), i ≠j, i,j = 1,2,...,n
(Subscript[c, ii]) = 0, i = 1,2,...,n
за вектора d:
(Subscript[d, i]) = ((Subscript[b, i])/(Subscript[a, ii])), i = 1,2,...,n

конструиране на метода за конкретния пример

разгънат вид:
123 (Subscript[x, 1]) + 0.2 (Subscript[x, 2])   = 34
-(Subscript[x, 1]) + 34 (Subscript[x, 2]) + (Subscript[x, 3]) = -21
	0.9 (Subscript[x, 2]) -211(Subscript[x, 3]) = -0.7

 за i-тия ред прехвърлям всичко, което е без (Subscript[x, i]) от дясната с обратен знак:
 123 (Subscript[x, 1])  =- 0.2 (Subscript[x, 2])  + 34
+ 34 (Subscript[x, 2])   =(Subscript[x, 1])  - (Subscript[x, 3])-21
-211(Subscript[x, 3]) =-0.9 (Subscript[x, 2]) -0.7

за i-тия ред делим на коефициента пред  	(Subscript[x, i])
(Subscript[x, 1])  =((- 0.2)/123) (Subscript[x, 2])  + (34/123)
(Subscript[x, 2])   =(1/34)(Subscript[x, 1])  -(1/34) (Subscript[x, 3])-(21/34)
(Subscript[x, 3]) =((-0.9)/(-211)) (Subscript[x, 2]) -(0.7/(-211))

в матричен вид:
x = C.x + d, където
C = (({{0, ((- 0.2)/123), 0}, {(1/34), 0, -(1/34)}, {0, ((-0.9)/(-211)), 0}})), d = (({{(34/123)}, {-(21/34)}, {-(0.7/(-211))}}))

програмна реализация на конструирането на метода за конкретния пример

A = (({{123, 0.2, 0}, {-1, 34, 1}, {0, 0.9, -211}})); b = {34,-21,-0.7};
n = Length[A];
(* инициализация матрицата C и вектора d  *)
c = (({{\[Placeholder], \[Placeholder], \[Placeholder]}, {\[Placeholder], \[Placeholder], \[Placeholder]}, {\[Placeholder], \[Placeholder], \[Placeholder]}}));d = {0,0,0};
For[i=1,i≤n,i++,
c[[i]]=- ((A[[i]])/(A[[i,i]]));
c[[i,i]]=0;
d[[i]] = ((b[[i]])/(A[[i,i]]))
]
c//MatrixForm
d//MatrixForm

за сравнение:
в матричен вид:
x = C.x + d, където
C = (({{0, ((- 0.2)/123), 0}, {(1/34), 0, -(1/34)}, {0, ((-0.9)/(-211)), 0}})), d = (({{(34/123)}, {-(21/34)}, {-(0.7/(-211))}}))

запис на итерационния процес за конкретния пример

((Subscript[x, 1])^((к+1)))  =((- 0.2)/123) ((Subscript[x, 2])^((к)))  + (34/123)
((Subscript[x, 2])^((к+1)))   =(1/34)((Subscript[x, 1])^((к)))  -(1/34) ((Subscript[x, 3])^((к)))-(21/34)
((Subscript[x, 3])^((к+1))) =((-0.9)/(-211)) ((Subscript[x, 2])^((к))) -(0.7/(-211))
к = 0,1, 2....

Проверка на сходимост

Ако поне една норма на матрицата C ||C||<1 то итерационния процес (x^((к+1))) = C.(x^((к))) + d, к = 0,1,2,...  ще е сходящ при всеки избор на начално приближение.

Norm[c]

Избор на начално приближение

ако е налице условието за сходимост, избираме произволно начално приближение

Уточняване на корена (извършване на итерациите)

x0 = {3,6,78};
Print["k = ", 0, " (x^((k))) = ", x0]
For[k = 1, k≤5, k++,
x1 = c.x0+d;
x0=x1;
Print["k = ", k, " (x^((k))) = ", x0]
]

пускаме итерациите при неподходящо начално приближение:

x0 = {3*(10^6),6000,-78*(10^10)};
Print["k = ", 0, " (x^((k))) = ", x0]
For[k = 1, k≤10, k++,
x1 = c.x0+d;
x0=x1;
Print["k = ", k, " (x^((k))) = ", x0]
]

за сравнение:
{0.277414, -0.609509, 0.000717734}

Оценка на грешката

Първи начин - чрез явна формула

(e^((k)))=

x0 = {3,6,78};nx0 = Norm[x0];
eps = ((Norm[c])^0)(nx0+((Norm[d])/(1-Norm[c])));
Print["k = ", 0, " (x^((k))) = ", x0, " (e^((k))) = ", eps]
For[k = 1, k≤5, k++,
x1 = c.x0+d;
x0=x1;
eps = ((Norm[c])^k)(nx0+((Norm[d])/(1-Norm[c])));
Print["k = ", k, " (x^((k))) = ", x0, " (e^((k))) = ", eps]
]

Втори начин - чрез разстоянието между две последователни приближения

x0 = {3,6,78};
Print["k = ", 0, " (x^((k))) = ", x0]
For[k = 1, k≤5, k++,
x1 = c.x0+d;
eps = Norm[x1-x0];
x0=x1;
Print["k = ", k, " (x^((k))) = ", x0, " (e^((k))) = ", eps]
]

Извършване на итерациите до достигане на предварително зададена точност

Първи начин - чрез явна формула

x0 = {3,6,78};nx0 = Norm[x0];
prec = (10^(-6));
eps = ((Norm[c])^0)(nx0+((Norm[d])/(1-Norm[c])));
Print["k = ", 0, " (x^((k))) = ", x0, " (e^((k))) = ", eps]
For[k = 1, eps>prec, k++,
x1 = c.x0+d;
x0=x1;
eps = ((Norm[c])^k)(nx0+((Norm[d])/(1-Norm[c])));
Print["k = ", k, " (x^((k))) = ", x0, " (e^((k))) = ", eps]
]

Втори начин - чрез разстоянието между две последователни приближения:

x0 = {3,6,78};
prec = (10^(-6));eps=1;
Print["k = ", 0, " (x^((k))) = ", x0]
For[k = 1, eps>prec, k++,
x1 = c.x0+d;
eps = Norm[x1-x0];
x0=x1;
Print["k = ", k, " (x^((k))) = ", x0, " (e^((k))) = ", eps]
]

Намиране на минимален брой итерации за достигане на предварително зададена точност

(e^((k)))<e => k> ((lg((e/(||(x^((0)))||+((||d||)/(1-||c||))))))/(lg(||c||)))

x0 = {3,6,78};prec = (10^(-6));
((Log10[(prec/(Norm[x0]+((Norm[d])/(1-Norm[c]))))])/(Log10[Norm[c]]))

Извод: Т.к. k е брой итерации, то взимаме най-малкото цяло число, което надвишава получения резултат => k=6

x0 = {3*(10^6),6000,-78*(10^10)};prec = (10^(-6));
((Log10[(prec/(Norm[x0]+((Norm[d])/(1-Norm[c]))))])/(Log10[Norm[c]]))