数值分析 之高斯消元法和它的改进

# -*- coding: utf-8 -*-
"""
Created on Sun Nov 19 16:40:16 2017

@author: www
"""
import numpy as np


def Gaussian_elimination(X,y):
    
    row,col = X.shape
    theta = np.zeros(row)
    y = y.reshape(row,1)
    
    X_y = np.hstack((X,y))
    
    for i in range(row):
        for j in range(i+1,row):
            X_y[j] = X_y[j] - X_y[i] * (X_y[j][i] / X_y[i][i])
            
    for i in range(row-1, -1, -1):
        sum = X_y[i][col]

        for j in range(i+1,row):
            sum = sum -X_y[i][j] * theta[j]

        theta[i] = sum / X_y[i][i]
    
    theta.resize(col, 1)
    return theta
    

def Gauss_better(X,y):
    row,col = X.shape
    theta = np.zeros(row)
    y = y.reshape(row,1)
    
    X_y = np.hstack((X,y))
    X_y_max = np.max(X_y, axis=1, keepdims=True)
    X_y_stand = np.zeros_like(X_y)
    
    for i in range(row):
        X_y_stand[i] = X_y[i] / X_y_max[i]

    for i in range(row):
        a = np.argmax(X_y_stand[i:, i]) + i
        if a != i:
            theta = np.copy(X_y[i])
            X_y[i] = X_y[a]
            X_y[a] = theta

    X = X_y[:, :-1]
    y = X_y[:, -1]

    Gaussian_elimination(X,y)            
        
    
            



def main():
    X = np.array([[2.,-1.,3.],
                [4., 2.,5.],
                [1., 2.,0.]])
    y = np.array([1., 4., 7.])
    theta = Gaussian_elimination(X,y)
    print(theta)
    
    theta = Gauss_better(X,y)
    print(theta)
    
main()