Гауссово исключение без поворота - PullRequest
Новая
       37

Гауссово исключение без поворота

0 голосов
/ 30 октября 2019

Я пытаюсь написать программу, которая может выполнять устранение по Гауссу без частичного поворота. Это то, что у меня есть, я знаю, что что-то напутал, но не могу понять, что. Он продолжает печатать ту же матрицу, с которой начал, а не 0 с 1 по диагонали. Любая помощь приветствуется! 

import numpy as np
import math

#A = np.array([[1,2,-1],[5,2,2],[-3,5,-1]])
#b = np.array([[2,9,1]])
A = np.array([[8,-2,-1,0,0],[-2,9,-4,-1,0],[-1,-3,7,-1,2],[0,-4,-2,12,-5],[0,0,-7,-3,15]])
b = np.array([[5],[2],[0],[1],[5]])

def forward_elim(A, b, n):
#calculates the forward part of Gaussian elimination.

    for row in range(0, n-1):
        for i in range(row+1, n):
            factor = A[i,row] / A[row,row]
            for j in range(row, n):
                A[i,j] = A[i,j] - factor * A[row,j]

            b[i] = b[i] - factor * b[row]

        print("A = \n%s and b = %s" % (A,b))
    return A, b

def back_sub(A, b, n):
#back substitution and returns result

    x = np.zeros((n,1))
    x[n-1] = b[n-1] / A[n-1, n-1]
    for row in range(n-2, -1, -1):
        sums = b[row]
        for j in range(row+1, n):
            sums = sums - A[row,j] * x[j]
        x[row] = sums / A[row,row]
    return x

def gauss_elim(A, b):
#performs gaussian elimination without pivoting

    n = A.shape[0]
    #checks for zero diagonal elements
    if any(np.diag(A)==0):
        raise ZeroDivisionError(("Can't divide by 0")) #raise used to raise exception

    A, b = forward_elim(A, b, n)
    return back_sub(A, b, n)

print(A)

1 Ответ

0 голосов
/ 30 октября 2019

Из приведенного кода не ясно, эффективно ли вы вызываете функцию gauss_elim(A, b).

Чтобы помочь выяснить, что происходит, попробуйте выполнить простой пример. Лучше всего с помощью отладчика или распечаткой промежуточных значений. Следуйте тому же примеру, записав шаги на листе бумаги.

Простой пример может быть: 

A = np.array([[2,4],[3,9]])
b = np.array([[2,30]])

Одна из проблем заключается в том, что вы объявляете b как 1xnматрица, но позже использовать его как просто вектор. Решением может быть объявление b как вектора: 

b = np.array([2,30])
...