using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace Methods_First_Task
{
#region Matrix_Class
class Matrix
{
double[,] data;
int n, k;
public Matrix(int n, int k)
{
this.n = n;
this.k = k;
data = new double[n, k];
for (int i = 0; i < n; ++i)
for (int j = 0; j < k; ++j)
data[i, j] = 0;
}
public Matrix(Matrix c)
{
this.n = c.n;
this.k = c.k;
data = new double[n, k];
for (int i = 0; i < n; ++i)
for (int j = 0; j < k; ++j)
data[i, j] = c.data[i, j];
}
public int getColumns()
{
return k;
}
public int getRows()
{
return n;
}
public void makeIdentity()
{
if (n != k)
{
Console.WriteLine("MAKE IDENTITY ERROR: N != K");
return;
}
for (int i = 0; i < n; ++i)
for (int j = 0; j < k; ++j)
{
data[i, j] = ((i == j) ? 1 : 0);
}
}
public Matrix(double[,] d)
{
n = d.GetLength(0);
k = d.GetLength(1);
data = new double[n, k];
for (int i = 0; i < n; ++i)
for (int j = 0; j < n; ++j)
data[i, j] = d[i, j];
}
public void setValue(double v, int i, int j)
{
data[i, j] = v;
}
public double getValue(int i, int j)
{
return data[i, j];
}
static public Matrix operator *(Matrix inp1, Matrix inp2)
{
if (inp1.k != inp2.n)
{
Console.WriteLine("MATRIX MULT ERROR");
return new Matrix(1, 1);
}
Matrix rs = new Matrix(inp1.n, inp2.k);
for (int i = 0; i < inp1.n; ++i)
{
for (int j = 0; j < inp2.k; ++j)
{
double tmp = 0;
for (int c = 0; c < inp1.k; ++c)
tmp += inp1.data[i, c] * inp2.data[c, j];
rs.setValue(tmp, i, j);
}
}
return rs;
}
static public Matrix operator *(Matrix inp1, double alpha)
{
Matrix rs = new Matrix(inp1.n, inp1.k);
for (int i = 0; i < inp1.n; ++i)
{
for (int j = 0; j < inp1.k; ++j)
{
rs.setValue(inp1.data[i, j] * alpha, i, j);
}
}
return rs;
}
static public Matrix operator -(Matrix inp1, double alpha)
{
Matrix rs = new Matrix(inp1.n, inp1.k);
for (int i = 0; i < inp1.n; ++i)
{
for (int j = 0; j < inp1.k; ++j)
{
rs.setValue(inp1.data[i, j] - alpha, i, j);
}
}
return rs;
}
static public Matrix operator -(Matrix inp1, Matrix inp2)
{
if (inp1.n != inp2.n || inp1.k != inp2.k)
{
Console.WriteLine("Operator- ERROR: Matr size not equal");
return new Matrix(1, 1);
}
Matrix rs = new Matrix(inp1.n, inp1.k);
for (int i = 0; i < inp1.n; ++i)
{
for (int j = 0; j < inp1.k; ++j)
{
rs.setValue(inp1.data[i, j] - inp2.data[i, j], i, j);
}
}
return rs;
}
static public Matrix operator +(Matrix inp1, Matrix inp2)
{
if (inp1.n != inp2.n || inp1.k != inp2.k)
{
Console.WriteLine("Operator+ ERROR: Matr size not equal");
return new Matrix(1, 1);
}
Matrix rs = new Matrix(inp1.n, inp1.k);
for (int i = 0; i < inp1.n; ++i)
{
for (int j = 0; j < inp1.k; ++j)
{
rs.setValue(inp1.data[i, j] + inp2.data[i, j], i, j);
}
}
return rs;
}
/// Tensor multiplication of vectors
static public Matrix operator ^(Matrix inp1, Matrix inp2)
{
if (inp1.k != inp2.n || inp1.n != inp2.k)
{
Console.WriteLine("Operator^ ERROR: Matr size not equal");
return new Matrix(1, 1);
}
Matrix rs = new Matrix(inp1.n, inp2.k);
for (int i = 0; i < inp1.n; ++i)
{
for (int j = 0; j < inp2.k; ++j)
{
rs.setValue(inp1.data[i, 0] * inp2.data[0, j], i, j);
}
}
return rs;
}
//транспонированная матрица
public Matrix getTranspozed()
{
Matrix ret = new Matrix(k, n);
for (int i = 0; i < n; ++i)
{
for (int j = 0; j < k; ++j)
{
ret.setValue(data[i, j], j, i);
}
}
return ret;
}
public void transpoze()
{
double[,] ret = new double[k, n];
for (int i = 0; i < n; ++i)
{
for (int j = 0; j < k; ++j)
{
ret[j, i] = data[i, j];
}
}
data = ret;
}
public Matrix getSubMatrix(int si, int sj, int ei, int ej)
{
if (si >= n || sj >= k || ei >= n || ej >= k ||
si < 0 || sj < 0 || ei < 0 || ej < 0)
{
Console.WriteLine("getSubMatrix RANGE ERROR");
return new Matrix(1, 1);
}
Matrix retMtr = new Matrix(ei - si + 1, ej - sj + 1);
int fi = 0;
int fj = 0;
for (int i = si; i <= ei; ++i)
{
fj = 0;
for (int j = sj; j <= ej; ++j)
{
retMtr.setValue(data[i, j], fi, fj);
++fj;
}
++fi;
}
return retMtr;
}
public bool isZeroMatrix()
{
for (int i = 0; i < n; ++i)
for (int j = 0; j < k; ++j)
if (data[i, j] != 0)
return false;
return true;
}
public void printMatrix()
{
for (int i = 0; i < n; ++i)
{
for (int j = 0; j < k; ++j)
{
Console.Write("{0, 3:f} ", data[i, j]);
}
Console.WriteLine();
}
Console.WriteLine(); Console.WriteLine(); Console.WriteLine();
}
public double getMaximumAboveDiag(ref int ii, ref int jj)
{
double val = -1000 * 1000 * 1000;
for (int i = 0; i < n; ++i)
{
for (int j = i + 1; j < k; ++j)
{
if (val < Math.Abs(data[i, j]))
{
val = Math.Abs(data[i, j]);
ii = i;
jj = j;
}
}
}
return val;
}
}
#endregion
class Methods
{
public double[,] rightSide;
int str = 0;
int stl = 0;
public Methods(int str, int stl)
{
rightSide = Methods.rightSideMatrix(str, stl);
this.str = str;
}
// единичная матрица
public static double[,] rightSideMatrix(int str, int stl)
{
double[,] rightSide = new double[str, stl];
for (int i = 0; i < str; i++)
{
for (int j = 0; j < stl; j++)
{
if (i == j)
rightSide[i, j] = 1;
else
rightSide[i, j] = 0;
}
}
return rightSide;
}
//public double[,] equationMatrix = new double[3, 3] { { 3, 45, 1 }, { 2, 15, 3 }, { 9, 8, 12 } };
// матрица A
public double[,] equationMatrix = new double[3, 3] { { 5, 1, 2 }, { 1, 4, 1 }, { 2, 1, 3 } };
// Метод Гаусса
public void gauss_solver()
{
int n = str;//this.rightSide.Length;
for (int i = 0; i < n; ++i)
for (int j = i + 1; j < n; ++j)
{
double value = equationMatrix[i, i];
double koeff = equationMatrix[j, i] / value;
for (int c = 0; c < n; ++c)
equationMatrix[j, c] -= koeff * equationMatrix[i, c];
// для всех элементов строки
for (int k = 0; k < n; k++)
{
rightSide[j, k] -= koeff * rightSide[i, k];
}
}
for (int i = n - 1; i >= 0; --i)
for (int j = i - 1; j >= 0; --j)
{
double value = equationMatrix[i, i];
double koeff = equationMatrix[j, i] / value;
for (int c = 0; c < n; ++c)
equationMatrix[j, c] -= koeff * equationMatrix[i, c];
for (int k = 0; k < n; k++)
{
rightSide[j, k] -= koeff * rightSide[i, k];
}
}
for (int i = 0; i < n; ++i)
{
for (int k = 0; k < n; k++)
{
rightSide[i, k] = rightSide[i, k] / equationMatrix[i, i];
}
}
}
static double getNorm(Matrix a)
{
return Math.Sqrt((a.getTranspozed() * a).getValue(0, 0));
}
//public double
//Matrix[] triada = new Matrix[3];
// метод скалярных произведений
public void smAlgorithm(Matrix A, Matrix y, ref Matrix x, ref double v, bool secondIter = false)
{
int n = A.getColumns();
double EPS = 0.000001;
double s = (y.getTranspozed() * y).getValue(0, 0);
double yNorm = Math.Sqrt(s);
Matrix lx = y * (1.0f / yNorm);
double lamdda = 0;
double lastLamda = 0;
for (int iter = 0; iter < 100000; ++iter)
{
y = A * lx;
if (secondIter)
y = y - x * ((x.getTranspozed() * y).getValue(0, 0));
s = (y.getTranspozed() * y).getValue(0, 0);
double t = (lx.getTranspozed() * y).getValue(0, 0);
yNorm = Math.Sqrt(s);
lx = y * (1.0 / yNorm);
lastLamda = lamdda;
lamdda = s / t;
if (Math.Abs(lastLamda - lamdda) < EPS)
break;
}
v = lamdda;
x = lx;
}
public void findSecondBySm(Matrix A)
{
Matrix y = new Matrix(A.getColumns(), 1);
Matrix x1 = new Matrix(A.getColumns(), 1);
double vl = 0;
for (int i = 0; i < A.getColumns(); ++i) y.setValue(1, i, 0);
smAlgorithm(A, y, ref x1, ref vl);
Console.WriteLine("Максимально собственное значение: {0}: ", vl);
Console.WriteLine("Собственный вектор максимального числа:");
x1.printMatrix();
for (int i = 0; i < A.getColumns(); ++i) y.setValue(123, i, 0);
y = y - x1 * ((x1.getTranspozed() * y).getValue(0, 0));
smAlgorithm(A, y, ref x1, ref vl, true);
//Console.WriteLine("adfasdf");
Console.WriteLine("Второе по величине собственное число: {0}", vl);
//Console.WriteLine(vl);
Console.WriteLine("Его собственный вектор:");
x1.printMatrix();
}
// Метод Якоби
public static void Jacobi(int n, double[][] A, double[] F, double[] X)
{
double EPS = 0.000001;
double[] TempX = new double[n];
double norm;
for (int i = 0; i < n; i++)
{
TempX[i] = X[i];
}
do
{
for (int i = 0; i < n; i++)
{
TempX[i] = F[i];
for (int g = 0; g < n; g++)
{
if (i != g)
TempX[i] -= A[i][g] * X[g];
}
TempX[i] /= A[i][i];
}
norm = Math.Abs(X[0] - TempX[0]);
for (int h = 0; h < n; h++)
{
if (Math.Abs(X[h] - TempX[h]) > norm)
norm = Math.Abs(X[h] - TempX[h]);
X[h] = TempX[h];
}
}
while (norm > EPS);
}
// Метод Зейделя
static bool ZeidelEnd(double[] xk, double[] xp, int n, double EPS)
{
double norm = 0;
for (int i = 0; i < n; i++)
{
norm += (xk[i] - xp[i]) * (xk[i] - xp[i]);
}
if (Math.Sqrt(norm) >= EPS)
return false;
return true;
}
public static void Zeidel(int n, double[,] a, double[] x, double[] p, double[] rightSide)
{
do
{
for (int i = 0; i < n; i++)
{
p[i] = x[i];
}
for (int i = 0; i < n; i++)
{
double var = 0;
for (int j = 0; j < i; j++)
{
var += (a[i, j] * x[j]);
}
for (int j = i + 1; j < n; j++)
{
var += (a[i, j] * p[j]);
}
x[i] = (rightSide[i] - var) / a[i,i];
}
}
while(!ZeidelEnd(x, p, n, 0.0001));
}
}
class Program
{
static void Main(string[] args)
{
double[,] test = Methods.rightSideMatrix(3, 3);
Methods a = new Methods(3, 3);
double[,] test2 = a.rightSide;
a.gauss_solver();
Console.WriteLine("Gauss:");
for (int i = 0; i < 3; i++)
{
for (int j = 0; j < 3; j++)
{
Console.Write(a.rightSide[i, j].ToString("0.00"));
Console.Write(" ");
}
Console.WriteLine();
}
double[,] res = new double[3, 3] { { 5, 1, 2 }, { 1, 4, 1 }, { 2, 1, 3 } };
Matrix eqMatrix = new Matrix(res);
a.findSecondBySm(eqMatrix);
// Метод Якоби
//double[][] JA = new double[3][] { new double[3]{ 1, 1, 1 }, new double[3]{4, 2, 1 }, new double[3]{9,3,1 } };
//double[] JFree = new double[] { 0, 1, 3 };
//double[] JRes = new double[3] {1,1,1};
//Methods.Jacobi(3, JA, JFree, JRes);
//Console.WriteLine("{0} {1} {2}", JRes[0], JRes[1], JRes[2]);
//Метод Зейделя
double[,] ZA = new double[3, 3] { { 1, 1, 1 }, { 4, 2, 1 }, { 9, 3, 1 } };
double[] ZFree = new double[] { 0, 1, 3 };
double[] ZRes = new double[3] { 0, 0, 0 };
double[] Zpred = new double[3];
Methods.Zeidel(3, ZA, ZRes, Zpred, ZFree);
Console.ReadLine();
}
}
}