CoCalc -- Collaborative Calculation in the Cloud
SharedWork6.ipynbOpen in CoCalc

Jupyter notebook Work6.ipynb

# -*- coding: utf-8 -*-
from __future__ import division, absolute_import, print_function
import numpy as np

def form_a(n, m):
    a = np.zeros((n, n))
    for i in range(1, n+1):
        for j in range(1, n+1):
            if i == j:
                a[i-1][j-1] = n + m**2 + j/m + i/n
            else:
                a[i-1][j-1] = (i + j)/(m + n)
    return a

def form_b(n):
    b = np.zeros(n)
    for i in range(1,n+1):
        b[i-1] = 200 + 50*i
    return b

def p_print(x):
    return str(x.tolist()).replace(']','').replace('[', '').replace(',','')

def cholesky(A):
    n = len(A)
    
    L = np.zeros((n, n))
    
    for i in range(n):
        for j in range(i + 1):
            t_sum = 0
            for k in range(j):
                t_sum += L[i][k] * L[j][k]
            
            if i == j:
                L[i][j] = np.sqrt(A[i][i] - t_sum)
            else:
                L[i][j] = 1.0 / L[j][j] * (A[i][j] - t_sum)
            
    return L

def solve_ls(L, b):
    n = len(L)
    
    y = np.matrix([0.0 for i in range(n)]).T
    for i in range(n):
        y[i, 0] = (b[i] - L[i, :i] * y[:i, 0]) / L[i, i]
    
    L = L.T
    x = np.matrix([0.0 for i in range(n)]).T
    for i in range(n - 1, -1, -1):
        x[i, 0] = (y[i] - L[i, i:] * x[i:, 0]) / L[i, i]
    return p_print(x)

if __name__ == "__main__":
    n = 20
    m = 8
    
    print(solve_ls(cholesky(form_a(n, m)), form_b(n)))
2.132160949534369 2.663065103035645 3.1926648789835226 3.720965078564629 4.247970479431361 4.773685835845899 5.298115878823181 5.821265316272813 6.343138833139963 6.863741091545193 7.3830767309233005 7.901150368161144 8.417966597734434 8.93352999184357 9.447845100548482 9.960916451902452 10.47274855208504 10.983345885533966 11.492712915076117 12.000854082057543