Теперь приведу свои листинги программ на Бейсике и Паскале.
1) для системы c вырожденной матрицей:
Код: Выбрать все
'Linear Equation System for Gilbert matrix and unity roots
'Program was written by Developer for E-Science Portal, 12.06.08
cls: n=12
print "Порядок матрицы N=";n: dim a#(n,n),b#(n),c#(n)
for i=1 to n
s#=0#
for j=1 to n
a#(i,j)=1/(i+j-1)
print using"#.###";a#(i,j);
print " ";
s#=s#+a#(i,j)
next j
b#(i)=s#: print using"#.###";b#(i)
next i 'Gilbert matrix
print "Epsilon=";: input t# 'Approximation accuracy
for j=1 to n
for k=j to n
s#=0: for i=1 to n
s#=s#+a#(i,j)*a#(i,k)
next i
c#(k)=s#
next k
c#=0: for i=1 to n
c#=c#+a#(i,j)*b#(i)
next i
for i=j to n: a#(i,j)=c#(i): next i
c#(j)=c#
next j
p#=sqr(t#*n)/2
for i=1 to n: a#(i,i)=a#(i,i)+p#: next i
a#(1,1)=sqr(a#(1,1))
for j=2 to n: a#(1,j)=a#(j,1)/a#(1,1): next j
for i=2 to n
s#=0: for k=1 to i-1
s#=s#+a#(k,i)*a#(k,i)
next k
a#(i,i)=sqr(a#(i,i)-s#)
for j=i+1 to n
s#=0: for k=1 to i-1: s#=s#+a#(k,i)*a#(k,j): next k
a#(i,j)=(a#(j,i)-s#)/a#(i,i)
next j
next i
c#(1)=c#(1)/a#(1,1)
for i=2 to n
s#=0: for k=1 to i-1
s#=s#+a#(k,i)*c#(k)
next k
c#(i)=(c#(i)-s#)/a#(i,i)
next i
c#(n)=c#(n)/a#(n,n)
for i=n-1 to 1 step -1
s#=0: for k=i+1 to n: s#=s#+a#(i,k)*c#(k): next k
c#(i)=(c#(i)-s#)/a#(i,i)
next i
s#=0
for i=1 to n
print "Корень X(";i;")=";
print using "#.#####";c#(i)
s#=s#+(c#(i)-1)^2
next i
s#=sqr(s#)/n
print "CKO=";s#;
end2) для системы по методу Гаусса:
Код: Выбрать все
cls: N=13
Print "Порядок матрицы N=";n: dim a#(n,n),b#(n),s#(n),c#(n),z%(n)
For i=1 to N
For j=1 to N
Print "a(";i;",";j;")=";
a#(i,j)=1/(i+j-1)
Print a#(i,j)
Next j
Next i 'Gilbert matrix
For i=1 to n
Print "B(";i;")=";
s#=0#: for j=1 to n
s#=s#+a#(i,j)
Next j
b#(i)=s#
Print b#(i)
Next i 'Gilbert matrix
det#=1# '-Determinant variable, initial value
For i=1 To N: z%(i)=i '-The numerating of the rows of parent matrix
Next i
For i=1 To N '-Big cycle of the inversing of parent matrix
k=i: y#=a#(i,i)
if (i+1)<=n then '-Selection of main elements}
for j=i+1 to n
w#=a#(i,j)
if (abs(w)>abs(y)) then
k=j: y#=w#
End If
Next j
End If
det#=det#*y# '-Determinant variable, new value
if (abs(y)<Eps) then
'-Standart algorithm terminates calculation if matrix is bad
Print "Matrix is bad (near to singular)" 'Stop
End If
y#=1/y#
for j=1 to n
c#(j)=a#(j,k): a#(j,k)=a#(j,i): a#(j,i)=-c#(j)*y#
b#(j)=a#(i,j)*y#: a#(i,j)=b#(j)
Next j
j=z%(i): z%(i)=z%(k): z%(k)=j: a#(i,i)=y#
for k=1 to n
if (k<>i) then
for j=1 to n
if (j<>i) then a#(k,j)=a#(k,j)-b#(j)*c#(k)
Next j
End If
Next k 'for k
Next i 'Big cycle For
for i=1 to n
Do
k =z%(i)
for j=1 to n
w#=a#(i,j): a#(i,j)=a#(k,j): a#(k,j)=w#
Next j
p=z%(i): z%(i)=z%(k): z%(k)=p
Loop Until(k=i)
Next i
Print "Determinant is: ";det#
Print "Inversed matrix:"
for i=1 to n
for k=1 to n
Print a#(i,k);: Print " ";
Next k
Print
Next i
cko#=0
for i=1 to n
s#(i)=0
for j=1 to n
s#(i)=s#(i)+b#(j)*a#(i,j)
next j
Print "Корни X(";i;")=";s#(i)
cko#=cko#+(s#(i)-1)^2
next i
cko#=sqr(cko#)/n
Print "CKO=";cko#
'Print "Определитель исходной матрицы D=",d#
End
Ha Паскале для решения уравнений c вырожденной матрицей (часть кода опущена), но программу можно "гонять" Turbo-Pascal 5.5 (6.0, 7.0):
Код: Выбрать все
Program LinEqSys; {-IDE Turbo-Pascal 6.0, Developer, http://e-science.ru}
{-Finds the roots of linear equation system (LES) AX=B with singular matrix A
by method of t-approximation: LES AX=B replace LES CX=D,
where |A-C|<=t, |B-D|<=t, 0<t<<1.
[A'A+aE][X]=[A'B], where
A - parent matrix,
A' - transposed matrix,
E - unitary matrix,
B - right part of LES,
a=0.5sqrt(n)t - parameter of approximation}
Uses Crt,Dos,Printer;
Const
MaxSize=50; {-Greate Gauss didn't solve equations more...}
Eps: Double=1e-6; {-Standard epsilon to detect of the singular matrix}
Type
Massive = Array[1..MaxSize] Of Extended{Double};
Matrix = Array[1..MaxSize, 1..MaxSize] Of Extended{Double};
Var
i,j,k,n,m: Byte; {-Counters of cycles and indexies of the elements of
arrays and matrixes}
b,g,x,roots: Massive; {-Right part of LES and roots}
a,c : Matrix; {-Parent and inverse matrixes}
d,s,t : Extended{Double}; {-Sumes, temporary storages for variable's values}
Procedure In_Put(Size,gMethod:Byte);
{-Reads data from disk file or generates them as Gilbert matrix or in a random way}
Var
i,j: Byte;
Procedure ReadDisk(Example:Byte); {-Reads data from disk file}
Begin {-To be implement if need}
Case Example of
0: Begin
N:=2;
a[1,1]:=999999999; a[1,2]:=-999999998; b[1]:=1;
a[2,1]:=999999998; a[2,2]:=-999999997; b[2]:=1;
End;
1: Begin
N:=3;
a[1,1]:=1; a[1,2]:=2; a[1,3]:=3; b[1]:=6;
a[2,1]:=2; a[2,2]:=1; a[2,3]:=0; b[2]:=3;
a[3,1]:=3; a[3,2]:=1; a[3,3]:=2; b[3]:=6
End;
2: Begin
N:=3;
a[1,1]:=-2; a[1,2]:=3; a[1,3]:=4; b[1]:=21;
a[2,1]:=1; a[2,2]:=5; a[2,3]:=7; b[2]:=-3;
a[3,1]:=4; a[3,2]:=2; a[3,3]:=1; b[3]:=-5
End;
3: Begin
N:=3;
a[1,1]:=2; a[1,2]:=1; a[1,3]:=3; b[1]:=6;
a[2,1]:=3; a[2,2]:=6; a[2,3]:=1; b[2]:=10;
a[3,1]:=8; a[3,2]:=-4;a[3,3]:=2; b[3]:=6
End;
4: Begin
N:=4;
a[1,1]:=1; a[1,2]:=2; a[1,3]:=3; a[1,4]:=4; b[1]:=10;
a[2,1]:=2; a[2,2]:=-2;a[2,3]:=-1;a[2,4]:=1; b[2]:=0;
a[3,1]:=1; a[3,2]:=1; a[3,3]:=1; a[3,4]:=1; b[3]:=4;
a[4,1]:=7; a[4,2]:=-1;a[4,3]:=1; a[4,4]:=-5;b[4]:=1
End;
5: Begin
N:=5;
a[1,1]:=1; a[1,2]:=2; a[1,3]:=3; a[1,4]:=4; a[1,5]:=5; b[1]:=15;
a[2,1]:=2; a[2,2]:=-2;a[2,3]:=-1;a[2,4]:=1; a[2,5]:=7; b[2]:=7;
a[3,1]:=1; a[3,2]:=1; a[3,3]:=1; a[3,4]:=1; a[3,5]:=1; b[3]:=5;
a[4,1]:=7; a[4,2]:=-1;a[4,3]:=1; a[4,4]:=-5;a[4,5]:=-1;b[4]:=1;
a[5,1]:=1; a[5,2]:=3; a[5,3]:=5; a[5,4]:=7; a[5,5]:=11;b[5]:=27
End;
6: Begin
N:=6;
a[1,1]:=1; a[1,2]:=2; a[1,3]:=3; a[1,4]:=4; a[1,5]:=5; a[1,6]:=6; b[1]:=21;
a[2,1]:=2; a[2,2]:=-2;a[2,3]:=-1;a[2,4]:=1; a[2,5]:=7; a[2,6]:=3; b[2]:=10;
a[3,1]:=1; a[3,2]:=1; a[3,3]:=1; a[3,4]:=1; a[3,5]:=1; a[3,6]:=5; b[3]:=10;
a[4,1]:=7; a[4,2]:=-1;a[4,3]:=1; a[4,4]:=-5;a[4,5]:=-1;a[4,6]:=-1;b[4]:=0;
a[5,1]:=1; a[5,2]:=3; a[5,3]:=5; a[5,4]:=7; a[5,5]:=11;a[5,6]:=13;b[5]:=40;
a[6,1]:=-1;a[6,2]:=2; a[6,3]:=-3;a[6,4]:=4; a[6,5]:=-5;a[6,6]:=6; b[6]:=3
End;
End
End;{ReadDisk}
Procedure RightPart; {-Calculates the right part of LES}
Var
i,j: Byte;
s : Extended{Double};
Begin
For i:=1 To Size Do {-Calculates the right part of LES}
Begin
s:=0;
For j:=1 To Size Do s:=s+a[i,j];
b[i]:=s
End
End;{RightPart}
Procedure gGilbert; {-Generates data as Gilbert matrix}
Var
i,j: Byte;
s : Extended{Double};
Begin
For i:=1 To Size Do {-Generates Gilbert matrix}
For j:=1 To Size Do a[i,j]:=1/(i+j-1);
RightPart
End;{gGilbert}
Procedure gRandom;{-Generates data in a random way}
Var
i,j: Byte;
Begin
Randomize;
For i:=1 To Size Do
For j:=1 To Size Do a[i,j]:=Random;
RightPart
End;{gRandom}
Begin {Main In_Put}
Case gMethod of
1: gGilbert;
2: gRandom
End{Case}
End;{In_Put}
Procedure Out_Put;
{-Displays the results of calculation on screen, disk file or printer}
Procedure WriteDisk; {-Writes data to disk file}
Begin {-To be implement if need}
End;{WriteDisk}
Procedure DisplayScreen; {-Displays data on screen}
Var
i,j : Byte; {-Local indexies, counters}
s,valSD : Extended{Double}; {-Sumes, temporary storages for variable's values}
Begin
s:=0;
For i:=1 To N Do s:=s+Sqr(Abs(roots[i])-1);
valSD:=Sqrt(s)/N;
Writeln;
Writeln('Standard deviation is:',valSD)
End;{DisplayScreen}
Procedure DisplayPrinter; {-Displays data on printer}
Begin {-To be implement if need}
End;{DisplayPrinter}
Begin{Main Out_Put}
DisplayScreen
End;{Out_Put}
Procedure Conversion;
{-Converses the parent matrix to the transposed matrix}
Var
i,j,k : Byte; {-Local indexies, counters}
{b,}c : Massive; {-Copy of right part of LES and the roots}
w,s,p : Extended{Double}; {-Sumes, temporary storages for variable's values}
Begin
For j:=1 To N Do {-Big cycle of the inversing of parent matrix}
Begin {-Forwaed and flyback tracies}
for k:=j to n Do
Begin
s:=0;
for i:=1 to N Do s:=s+a[i,j]*a[i,k];
c[k]:=s
End;
w:=0; for i:=1 To N Do w:=w+a[i,j]*b[i];
For i:=j To N Do a[i,j]:=c[i];
c[j]:=w
End;{j}
p:=sqrt(t*N)/2;
for i:=1 to N Do a[i,i]:=a[i,i]+p;
a[1,1]:=sqrt(a[1,1]);
For j:=2 To N Do a[1,j]:=a[j,1]/a[1,1];
For i:=2 To N Do
Begin
s:=0;
For k:=1 To i-1 Do s:=s+a[k,i]*a[k,i];
a[i,i]:=sqrt(a[i,i]-s);
For j:=i+1 To N Do
Begin
s:=0; For k:=1 To i-1 Do s:=s+a[k,i]*a[k,j];
a[i,j]:=(a[j,i]-s)/a[i,i]
End;{j}
End;{i}
c[1]:=c[1]/a[1,1];
for i:=2 to N Do
Begin
s:=0; For k:=1 To i-1 Do s:=s+a[k,i]*c[k];
c[i]:=(c[i]-s)/a[i,i]
End;
c[N]:=c[N]/a[N,N];
for i:=N-1 DownTo 1 Do
Begin
s:=0; For k:=i+1 To N Do s:=s+a[i,k]*c[k];
c[i]:=(c[i]-s)/a[i,i] {[C] - roots of LES}
End; roots:=c
End;{Conversion}
Procedure Decision;
{-Calculates the roots of LES}
Var
i,j : Byte;
g,s : Massive; {-Right part of LES and roots}
a_ : Matrix; {-Parent and inverse matrixes}
Begin
g:=b; a_:=a; {-Creates copies parent matrix and right part of LES}
Writeln('System roots are:');
for i:=1 to n Do
Begin
s[i]:=0; for j:=1 to n Do s[i]:=s[i]+g[j]*a_[i,j];
Write('X(',i:1,')=',s[i]:2:17,' ');
End; roots:=s
End;{Decision}
Begin{Main}
ClrScr;
N:=13{50};
In_Put(N,1);
Write('t=');
Readln(t);
{ N:=2;
a[1,1]:=999999999; a[1,2]:=-999999998; b[1]:=1;
a[2,1]:=999999998; a[2,2]:=-999999997; b[2]:=1;}
{ Writeln('Parent matrix:',' ':20,'Right part of:');
For i:=1 To N Do
Begin
For j:=1 To N Do
Begin
Write(a[i,j]:2:4); Write(' ')
End;
Writeln(' ',b[i]:2:4)
End;}
Conversion;
for i:=1 to N do writeln('X(',i,')=',roots[i]:2:5);
{Decision;}
Out_Put
End.{LinEqSys}
Всякие отладочные и промежуточные "штучки" я из листинга не исключил, - они взяты в фигурные скобки комментариев...
