fork(7) download
copy 
  1. import numpy as np
  2.  
  3. def forward_sub(L, b):
  4.     """Given a lower triangular matrix L and right-side vector b,
  5.     compute the solution vector y solving Ly = b."""
  6.  
  7.     y = np.zeros_like(b)
  8.     for i in range(len(b)):
  9.       s = np.dot(L[i, :i], y[:i])
  10.       y[i] = b[i] - s
  11.     return y
  12.  
  13. def backward_sub(U, y):
  14.     """Given a lower triangular matrix U and right-side vector y,
  15.     compute the solution vector x solving Ux = y."""
  16.  
  17.     x = np.zeros_like(y)
  18.  
  19.     for i in range(len(x), 0, -1):
  20.       x[i-1] = (y[i-1] - np.dot(U[i-1, i:], x[i:])) / U[i-1, i-1]
  21.  
  22.     return x
  23.  
  24. def lu_factor(A):
  25.  
  26.   #LU decompostion using Doolittles method
  27.  
  28.   L = np.zeros_like(A)
  29.   U = np.zeros_like(A)
  30.  
  31.   N = np.size(A,0)
  32.  
  33.   """for i in range(N):
  34.     for k in range(i, N):
  35.       U[i,k] = A[i,k] - np.dot(L[i, :i], U[:i, k])
  36.     L[i,i] = 1
  37.     for k in range(i+1, N):
  38.       L[k,i] = (A[k,i] - np.dot(L[k, :i], U[:i, i])) / U[i,i]"""
  39.   for k in range(N):
  40.     L[k, k] = 1
  41.     U[k, k] = (A[k, k] - np.dot(L[k, :k], U[:k, k]))
  42.     for j in range(k+1, N):
  43.       U[k, j] = (A[k, j] - np.dot(L[k, :k], U[:k, j]))
  44.     for i in range(k+1, N):
  45.       L[i, k] = (A[i, k] - np.dot(L[i, :k], U[:k, k])) / U[k, k]
  46.   return (L, U)
  47.  
  48.  
  49. def lu_solve(L, U, b):
  50.     # Step 1: Solve Uy = b using forward substitution
  51.  
  52.     # Step 2: Solve Lx = y using backward substitution
  53.  
  54.     y = forward_sub(L,b)
  55.     x = backward_sub(U,y)
  56.  
  57.     return x
  58.  
  59.  
  60. def linear_solve(A, b):
  61.     # ...
  62.  
  63.     L, U = lu_factor(A)
  64.     print(L)
  65.     print(U)
  66.     print(L @ U)
  67.     x = lu_solve(L,U,b)
  68.     return x
  69.  
  70.  
  71. b = np.array([6.0,-4,27])
  72. A = np.matrix([[1.,1,1],
  73.                [0,2,5],
  74.                [2,5,-1]])
  75.  
  76. print(linear_solve(A,b))
  77.  
Success #stdin #stdout 0.19s 27444KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
[[1.  0.  0. ]
 [0.  1.  0. ]
 [2.  1.5 1. ]]
[[  1.    1.    1. ]
 [  0.    2.    5. ]
 [  0.    0.  -10.5]]
[[ 1.  1.  1.]
 [ 0.  2.  5.]
 [ 2.  5. -1.]]
[ 5.  3. -2.]
https://ideone.com/2mhhf7
language:
Python 3 (python  3.12)
created:
6 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!