Методы решения систем линейных алгебраических уравнений

Название	Методы решения систем линейных алгебраических уравнений
Дата	22.11.2018
Размер	151.49 Kb.
Формат файла
Имя файла	Chislennye_metody_reshenie_SLAU_Bagramov_Robert_09-405_Avtosokhr.docx
Тип	Решение #57324
страница	10 из 10

1 2 3 4 5 6 7 8 9 10

Сравнение методов Якоби, Зейделя, верхней релаксации.

Эксперимент 1.

Сравним данные методы по количеству итераций.

n	Количество итераций
n	Метод Якоби	Метод Зейделя	Метод релаксации
5	12	11	1
10	73	53	1
15	153	136	1
20	345	262	2
25	484	435	3
30	845	656	5

Вывод: Из графиков видно, что метод Якоби наиболее затратнее по количеству итераций, чем метод Зейделя. Также видно как метод релаксации превосходит предыдущие два метода по скорости достижения заданной точности, при хорошем выборе веса.

Листинг.

import org.apache.poi.hssf.usermodel.HSSFWorkbook;
import org.apache.poi.ss.usermodel.Cell;
import org.apache.poi.ss.usermodel.Row;
import org.apache.poi.ss.usermodel.Sheet;
import org.apache.poi.ss.usermodel.Workbook;

import java.io.BufferedOutputStream;
import java.io.File;
import java.io.FileOutputStream;
import java.io.IOException;
import java.util.Scanner;

/**
* 11.10.2017
*
* @author Robert Bagramov.
*/
public class NumericalMethods {

public static double[] gradientDescent(double matA[][], double matB[], double y0[]) {
int n = matB.length;
double y1[] = new double[n];
double r;
double t = 0;
double den = 0;
double eps = 1.0 / (n * n * n);
int counter = 0;

do {
double[] rk = new double[n];
double[] rA = new double[n];

for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
rk[i] += matA[i][j] * y0[j];
}
rk[i] -= matB[i];
t += rk[i] * rk[i];
}

for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
rA[i] += matA[i][j] * rk[j];
}
}

for (int i = 0; i < n; i++) {
den += rA[i] * rk[i];
}

t = t / den;

for (int i = 0; i < n; i++) {
y1[i] = y0[i] - t * rk[i];
}

r = Math.abs(y1[0] - y0[0]);
for (int z = 1; z < n - 1; z++) {
if (r < Math.abs(y1[z] - y0[z])) {
r = Math.abs(y1[z] - y0[z]);
}
}

for (int i = 0; i < n; i++) {
y0[i] = y1[i];
}

t = 0;
den = 0;

counter++;

} while (r >= eps);
System.out.println(n + " " + counter + " " + eps + " " + r + " " + eps * n);

System.out.println(counter);

return y1;
}

public static double[] relaxationAlgorithm(double matA[][], double matB[], double y0[], double w, double h) {
double a = 0.0, r;
int n = matB.length;
double eps = h;
double y1[] = new double[n];

int counter = 0;

do {
for (int i = 0; i < n; i++) {

for (int j = 0; j < i; j++) {
a += matA[i][j] / matA[i][i] * y1[j];
}

for (int j = i + 1; j < n; j++) {
a += matA[i][j] / matA[i][i] * y0[j];
}

y1[i] = (1 - w) * y0[i] + w * (matB[i] / matA[i][i] - a);
a = 0.0;

}

r = Math.abs(y1[0] - y0[0]);
for (int z = 1; z < n; z++) {
if (r < Math.abs(y1[z] - y0[z])) {
r = Math.abs(y1[z] - y0[z]);
}
}

for (int i = 0; i < n; i++) {
y0[i] = y1[i];
}

counter++;
} while (r >= eps);
r = Math.abs(y0[0] - u(h));

for (int z = 1; z < n; z++) {
if (r < Math.abs(y0[z] - u((z + 1) * h))) {
r = Math.abs(y0[z] - u((z + 1) * h));
}
}
excelWork(n, counter, eps, r, h * h);
System.out.println(n+1);
System.out.println(w);
System.out.println(eps);
System.out.println(r);
System.out.println(h * h);
System.out.println(counter);
System.out.println("------------------");
return y0;
}

public static double[] seidelAlgorithm(double matA[][], double matB[], double y0[], double h) {
double a = 0, r;
int n = matB.length;
double eps = h * h * h * h * h;
double y1[] = new double[n];

int counter = 0;

do {
for (int i = 0; i < n; i++) {

for (int j = 0; j < i; j++) {

a += matA[i][j] / matA[i][i] * y1[j];
}
for (int j = i + 1; j < n; j++) {

a += matA[i][j] / matA[i][i] * y0[j];
}
y1[i] = matB[i] / matA[i][i] - a;

a = 0;
}

r = Math.abs(y1[0] - y0[0]);

for (int z = 1; z < n; z++) {
if (r < Math.abs(y1[z] - y0[z])) {
r = Math.abs(y1[z] - y0[z]);
}
}

for (int i = 0; i < n; i++) {
y0[i] = y1[i];
}

counter++;
} while (r >= eps);

r = Math.abs(y0[0] - u(h));

for (int z = 1; z < n; z++) {
if (r < Math.abs(y0[z] - u((z + 1) * h))) {
r = Math.abs(y0[z] - u((z + 1) * h));
}
}
excelWork(n, counter, eps, r, h * h);
System.out.println(n + " " + counter + " " + eps + " " + r + " " + h * h);
return y0;
}

public static double[] jacobiAlgorithm(double matA[][], double matB[], double y0[], double h) {
System.out.println("Jacobi Algorithm");
double a = 0.0, r;
int n = matB.length;
double eps = h * h * h;
double y1[] = new double[n];

int counter = 0;

do {
for (int i = 0; i < n; i++) {

for (int j = 0; j < n; j++) {
if (i == j) {
continue;
}
a += matA[i][j] / matA[i][i] * y0[j];
}

y1[i] = matB[i] / matA[i][i] - a;

a = 0.0;
}

r = Math.abs(y1[0] - y0[0]);

for (int z = 1; z < n; z++) {
if (r < Math.abs(y1[z] - y0[z])) {
r = Math.abs(y1[z] - y0[z]);
}
}

for (int i = 0; i < n; i++) {
y0[i] = y1[i];
}

counter++;
} while (r >= eps);

r = Math.abs(y0[0] - u(h));

for (int z = 1; z < n; z++) {
if (r < Math.abs(y0[z] - u((z + 1) * h))) {
r = Math.abs(y0[z] - u((z + 1) * h));
}
}
excelWork(n, counter, eps, r, h * h);
System.out.println(n + " " + counter + " " + eps + " " + r + " " + h * h);
return y0;
}

public static double[] tridiagonalMatrixAlgorithm(double matA[][], double matB[]) {
int n = matA.length;
int n1 = n - 1;

double alpha[] = new double[n];
double betta[] = new double[n];
double resY[] = new double[n];
double d = matA[0][0];

alpha[0] = -matA[0][1] / d;
betta[0] = matB[0] / d;

for (int i = 1; i < n1; i++) {
d = matA[i][i - 1] * alpha[i - 1] + matA[i][i];
alpha[i] = -matA[i][i + 1] / d;
betta[i] = (matB[i] - matA[i][i - 1] * betta[i - 1]) / d;
}

resY[n1] = (matB[n1] - matA[n1][n1 - 1] * betta[n1 - 1]) / (matA[n1][n1] + matA[n1][n1 - 1] * alpha[n1 - 1]);
for (int i = n1 - 1; i >= 0; i--) {
resY[i] = alpha[i] * resY[i + 1] + betta[i];
}

return resY;
}

public static double p(double x) {
return (1.0 + x);
}

public static double u(double x) {
return x * (1.0 - x);
}

public static double q(double x) {
return x + 1.0;
}

public static double f(double x) {
return (1.0 + 5.0 * x - x * x * x);
}

public static void excelWork(int n, int counter, double eps, double r, double hh) {
Workbook book = new HSSFWorkbook();

Sheet sheet = book.createSheet("Numerical methods");

Row row = sheet.createRow(0);

Cell nulCell = row.createCell(0);
nulCell.setCellValue(n + 1);
Cell firstCell = row.createCell(1);
firstCell.setCellValue(counter);
Cell secondCell = row.createCell(2);
secondCell.setCellValue(eps);
Cell thirdCell = row.createCell(3);
thirdCell.setCellValue(r);
Cell fifthCell = row.createCell(4);
fifthCell.setCellValue(hh);

try (BufferedOutputStream fos = new BufferedOutputStream(new FileOutputStream(new File("Table0.xls")))) {
book.write(fos);
} catch (IOException e) {
e.printStackTrace();
}

}

public static void showMatrix(double matrix[][]) {
for (int i = 0; i < matrix.length; i++) {
for (int j = 0; j < matrix.length; j++) {
System.out.print(matrix[i][j] + " ");
}
System.out.println();
}
}

public static double[] createFreeMembersVector(int n, double h) {
double matB[] = new double[n - 1];

for (int i = 0; i < n - 1 ; i++) {
matB[i] = f((i + 1) * h) * h * h;
}

return matB;
}

public static double[][] createDiagonalMatrix(int n, double h) {

double matA[][] = new double[n - 1][n - 1];

matA[0][0] = p(h) + p(2 * h) + h * h * q(h);
matA[0][1] = -p(2 * h);

for (int i = 1; i < n - 2; i++) {
matA[i][i - 1] = -p((i + 1) * h);

matA[i][i] = p((i + 1) * h) + p((i + 2) * h) + h * h * q((i + 1) * h);

matA[i][i + 1] = -p((i + 2) * h);
}

matA[n - 2][n - 2] = p((n - 1) * h) + p(n * h) + h * h * q((n - 1) * h);
matA[n - 2][n - 3] = -p((n - 1) * h);

return matA;
}

public static void main(String[] args) {
Scanner in = new Scanner(System.in);
int n = in.nextInt();

double h = 1.0 / n;

double matA[][];
double matB[];
double resY[] = new double[n - 1];

double matDiagTest[] = new double[n - 1];

matA = createDiagonalMatrix(n, h);
matB = createFreeMembersVector(n, h);

double r = 0;

for (int i = 0; i < n - 1; i++) {
for (int j = 0; j < n - 1; j++) {

if (i == j) {
continue;
}
r += matA[i][j];
}

matDiagTest[i] = matA[i][i] + r;
r = 0;
}

// System.out.println();
// System.out.println("mat diag test");
// for (int i = 0; i < n - 1; i++) {
// System.out.println(matDiagTest[i]);
// }
// System.out.println();

////---------------------------------------------------------------------------//
// resY = tridiagonalMatrixAlgorithm(matA, matB);
//
// for (int i = 0; i < n-1; i++) {
// System.out.println(resY[i]);//- u(i * h)
// }
//---------------------------------------------------------------------------//
// System.out.println();
//
double[] y0 = new double[n - 1];
for (int i = 0; i < n - 1; i++) {
y0[i] = 0.5;//matB[i]/matA[i][i];
}

//---------------------------------------------------------------------------//
//
resY = jacobiAlgorithm(matA, matB, y0, h);

// for (int i = 0; i < n - 1; i++) {
//
// System.out.println(Math.abs(resY[i] - u(i * h)));//
// }
////---------------------------------------------------------------------------//
// System.out.println();
////---------------------------------------------------------------------------//
// resY = seidelAlgorithm(matA, matB, y0, h );

// for (int i = 0; i < n; i++) {
// System.out.println(Math.abs(resY[i]));//- u(i * h)
// }
////---------------------------------------------------------------------------//
// System.out.println();
////---------------------------------------------------------------------------//

// for (double w = 1.1; w < 2; w = w + 0.1) {
// resY = relaxationAlgorithm(matA, matB, y0, w, h);
//// System.out.println(w);
// }
//
// for (int i = 0; i < n; i++) {
// System.out.println(Math.abs(resY[i]));//- u(i * h)
// }
//////---------------------------------------------------------------------------//
// System.out.println();
// resY = gradientDescent(matA, matB, y0);

// for (int i = 0; i < n; i++) {
// System.out.println(Math.abs(u(i * h)));//- u(i * h)
// }
//
// for (int i = 0; i < n - 1; i++) {
// System.out.println(resY[i]);
// }
// System.out.println();
// for (int i = 1; i < n; i++) {
// System.out.println(u(i * h));//-
// }
}
}

1 2 3 4 5 6 7 8 9 10