Action disabled: source
num-ansys:start
Table of Contents
Numerical Analysis
Solving Linear Equatuion
The following algorithm solves the equation Ax=b.
import math def factorize(A): n = int(math.sqrt(len(A))) diagA = True L = [] # check if no diagonal entries are zero for i in range(0, len(A), n + 1): if A[i] == 0: diagA = False # factorize with no row exchanges if diagA: for i in range(len(A)): r, c = i2rc(i, n) if r == c: # diagonal entries L.append(1) continue elif r > c: # lower triangular entries ml = A[i] / A[n * (c - 1) + c - 1] L.append(ml) r0 = A[n * (r - 1):n * r] r1 = A[n * (c - 1):n * c] A[n * (r - 1):n * r] = rewrite(r0, r1, ml) continue elif r < c: # upper triangular entries L.append(0) continue U = A return (L, U) def rewrite(r0, r1, ml): for i in range(len(r0)): r0[i] = r0[i] - ml * r1[i] return r0 def rc2i(r, c, n): i = n * (r - 1) + c - 1 return i def i2rc(i, n): r = i // n + 1 c = i % n + 1 return (r, c) def solv_lw_tri(L, b): n = len(b) c = [] for i in range(n): c.append(b[i]) for j in range(i): c[i] = c[i] - L[rc2i(i + 1, j + 1, n)] * c[j] c[i] = c[i] / L[rc2i(i + 1, i + 1, n)] return c def solv_up_tri(U, c): n = len(c) x = [0] * n for i in range(n - 1, -1, -1): x[i] = c[i] for j in range(i + 1, n): x[i] = x[i] - U[rc2i(i + 1, j + 1, n)] * x[j] x[i] = x[i] / U[rc2i(i + 1, i + 1, n)] return x class Test(object): """Test methods""" A = [2, 1, 3, 4, 3, 10, -2, 1, 7] b = [11, 28, 3] L = [1, 0, 0, 2, 1, 0, -1, 2, 1] U = [2, 1, 3, 0, 1, 4, 0, 0, 2] c = [11, 6, 2] x = [3, 2, 1] def test_factorize(self, A): assert factorize(A) == (self.L, self.U) def test_solv_lw_tri(self, L, b): assert solv_lw_tri(L, b) == self.c def test_solv_up_tri(self, U, c): assert solv_up_tri(U, c) == self.x if __name__ == "__main__": test = Test() test.test_factorize(test.A) test.test_solv_lw_tri(test.L, test.b) test.test_solv_up_tri(test.U, test.c)
Polynomial curve fitting
import numpy as np def polyfit(x, fx, w, n): """Calculate the coeffients of the polynomial that best approximate (x, y). :param x (list): indipendent variable :param fx (list): dipendent variable :param w (list): weight of each node :param n (int): degree of the polynomial :return a (list): coeffientes of the polynomial """ A = [] for i in range(len(x)): # rows r = [] for j in range(n, 0, -1): # columns r.append(w[i]**.5 * x[i]**(n-j)) A.append(r) b = [] for i in range(len(x)): # columns b.append([w[i]**.5 * fx[i]]) A = np.array(A) b = np.array(b) a = np.linalg.lstsq(np.dot(A.T, A), np.dot(A.T, b)) # (A' * A) \ (A' * b) a = a[0].T.tolist()[0] return a
num-ansys/start.txt · Last modified: 2023/05/28 16:07 by 127.0.0.1