calculus/calculus.go - master - mill/emp

calculus/calculus.go
83 строки3 KB
EdieLay
Улучшения в wave, реализация задачи теплицы
30 апр 2024, 21:21
30 апр 2024, 21:21acb6207
package calculus

import (
	"errors"
	"math"
)

// Функция, которая принимает квадратную матрицу коэффициентов системы линейных уравнений, столбец свободных членов и решает эту систему методом прогонки
func SolveTridiagonalSystem(matrix [][]float64, freeTerms []float64) ([]float64, error) {
	var n = len(matrix)
	if len(matrix[0]) != n {
		err := errors.New("matrix is not square")
		result := []float64{0}
		return result, err
	}
	if len(freeTerms) != n {
		err := errors.New("length of free terms doesn't match matrix size")
		result := []float64{0}
		return result, err
	}
	for i := 0; i < n; i++ {
		for j := 0; j < n; j++ {
			if math.Abs(float64(i-j)) > 1 && matrix[i][j] != 0 {
				err := errors.New("matrix is not tridiagonal")
				result := []float64{0}
				return result, err
			}
		}
	}

	// https://pro-prof.com/forums/topic/sweep-method-for-solving-systems-of-linear-algebraic-equations
	// прямой ход метода прогонки
	y, alpha, beta := make([]float64, n), make([]float64, n), make([]float64, n)
	y[0] = matrix[0][0]
	alpha[0] = -matrix[0][1] / y[0]
	beta[0] = freeTerms[0] / y[0]
	for i := 1; i < n-1; i++ {
		y[i] = matrix[i][i] + alpha[i-1]*matrix[i][i-1]
		alpha[i] = -matrix[i][i+1] / y[i]
		beta[i] = (freeTerms[i] - beta[i-1]*matrix[i][i-1]) / y[i]
	}
	y[n-1] = matrix[n-1][n-1] + alpha[n-2]*matrix[n-1][n-2]
	beta[n-1] = (freeTerms[n-1] - beta[n-2]*matrix[n-1][n-2]) / y[n-1]

	// обратный ход метода прогонки
	x := make([]float64, n)
	x[n-1] = beta[n-1]
	for i := n - 2; i >= 0; i-- {
		x[i] = beta[i] + alpha[i]*x[i+1]
	}

	return x, nil // возвращаем результат
}

func FillEquationSystemForEMPWithFDM(mu float64, equationsNumber int, leftOdd bool, rightOdd bool) [][]float64 {
	equationSystem := make([][]float64, equationsNumber) // матрица системы уравнений
	for i := 0; i < equationsNumber; i++ {
		equationSystem[i] = make([]float64, equationsNumber)
	}

	// заполнение матрицы коэффициентов
	for i := 0; i < equationsNumber; i++ {
		for j := 0; j < equationsNumber; j++ {
			if i == j {
				equationSystem[i][j] = 1 + 2*mu
			} else if i == j+1 || i == j-1 {
				equationSystem[i][j] = -mu
			} else {
				equationSystem[i][j] = 0
			}
		}
	}

	// Добавляем граничные условия
	if !leftOdd {
		equationSystem[0][0] -= mu
	}
	if !rightOdd {
		equationSystem[equationsNumber-1][equationsNumber-1] -= mu
	}

	return equationSystem
}