Sistemas de Ecuaciones Lineales Especiales#
import numpy as np
import numpy.linalg as la
help(range)
Help on class range in module builtins:
class range(object)
| range(stop) -> range object
| range(start, stop[, step]) -> range object
|
| Return an object that produces a sequence of integers from start (inclusive)
| to stop (exclusive) by step. range(i, j) produces i, i+1, i+2, ..., j-1.
| start defaults to 0, and stop is omitted! range(4) produces 0, 1, 2, 3.
| These are exactly the valid indices for a list of 4 elements.
| When step is given, it specifies the increment (or decrement).
|
| Methods defined here:
|
| __bool__(self, /)
| True if self else False
|
| __contains__(self, key, /)
| Return bool(key in self).
|
| __eq__(self, value, /)
| Return self==value.
|
| __ge__(self, value, /)
| Return self>=value.
|
| __getattribute__(self, name, /)
| Return getattr(self, name).
|
| __getitem__(self, key, /)
| Return self[key].
|
| __gt__(self, value, /)
| Return self>value.
|
| __hash__(self, /)
| Return hash(self).
|
| __iter__(self, /)
| Implement iter(self).
|
| __le__(self, value, /)
| Return self<=value.
|
| __len__(self, /)
| Return len(self).
|
| __lt__(self, value, /)
| Return self<value.
|
| __ne__(self, value, /)
| Return self!=value.
|
| __reduce__(...)
| Helper for pickle.
|
| __repr__(self, /)
| Return repr(self).
|
| __reversed__(...)
| Return a reverse iterator.
|
| count(...)
| rangeobject.count(value) -> integer -- return number of occurrences of value
|
| index(...)
| rangeobject.index(value) -> integer -- return index of value.
| Raise ValueError if the value is not present.
|
| ----------------------------------------------------------------------
| Static methods defined here:
|
| __new__(*args, **kwargs)
| Create and return a new object. See help(type) for accurate signature.
|
| ----------------------------------------------------------------------
| Data descriptors defined here:
|
| start
|
| step
|
| stop
Sustitución hacia adelante#
def forsub(L, bs):
"""
n: size de x
"""
n = bs.size
xs = np.zeros(n)
for i in range(n):
# i recorre desde 0 hasta n-1
xs[i] = (bs[i] - L[i, :i] @ xs[:i]) / L[i, i]
return xs
help(np.tril)
Help on _ArrayFunctionDispatcher in module numpy:
tril(m, k=0)
Lower triangle of an array.
Return a copy of an array with elements above the `k`-th diagonal zeroed.
For arrays with ``ndim`` exceeding 2, `tril` will apply to the final two
axes.
Parameters
----------
m : array_like, shape (..., M, N)
Input array.
k : int, optional
Diagonal above which to zero elements. `k = 0` (the default) is the
main diagonal, `k < 0` is below it and `k > 0` is above.
Returns
-------
tril : ndarray, shape (..., M, N)
Lower triangle of `m`, of same shape and data-type as `m`.
See Also
--------
triu : same thing, only for the upper triangle
Examples
--------
>>> np.tril([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], -1)
array([[ 0, 0, 0],
[ 4, 0, 0],
[ 7, 8, 0],
[10, 11, 12]])
>>> np.tril(np.arange(3*4*5).reshape(3, 4, 5))
array([[[ 0, 0, 0, 0, 0],
[ 5, 6, 0, 0, 0],
[10, 11, 12, 0, 0],
[15, 16, 17, 18, 0]],
[[20, 0, 0, 0, 0],
[25, 26, 0, 0, 0],
[30, 31, 32, 0, 0],
[35, 36, 37, 38, 0]],
[[40, 0, 0, 0, 0],
[45, 46, 0, 0, 0],
[50, 51, 52, 0, 0],
[55, 56, 57, 58, 0]]])
L = np.array([[1, 3, 8], [-2, np.pi, 0], [2, -1, 1]])
L
array([[ 1. , 3. , 8. ],
[-2. , 3.14159265, 0. ],
[ 2. , -1. , 1. ]])
L = np.tril(L)
L
array([[ 1. , 0. , 0. ],
[-2. , 3.14159265, 0. ],
[ 2. , -1. , 1. ]])
i = 2
L[i, :i] # desde 0 hasta i - 1
array([ 2., -1.])
b = np.array([1, 2, 3])
b
array([1, 2, 3])
b[:2]
array([1, 2])
x = np.zeros(3)
i = 1
x[i] = (b[i] - L[i, :i] @ x[:i]) / L[i, i]
L[i, :i]
array([-2.])
x[:1]
array([0.])
x
array([0. , 0.63661977, 0. ])
for i in range(0, 3):
print(L[i, i])
1.0
3.141592653589793
1.0
B = np.tril(np.random.rand(3, 3))
B
array([[0.65162791, 0. , 0. ],
[0.35777475, 0.43116038, 0. ],
[0.00869951, 0.83185453, 0.37251968]])
Sustitución hacia atrás#
def backsub(U, bs):
n = bs.size
xs = np.zeros(n)
for i in reversed(range(n)):
xs[i] = (bs[i] - U[i, i + 1 :] @ xs[i + 1 :]) / U[i, i]
return xs
for i in reversed(range(5)):
print(i)
4
3
2
1
0
help(np.triu)
Help on _ArrayFunctionDispatcher in module numpy:
triu(m, k=0)
Upper triangle of an array.
Return a copy of an array with the elements below the `k`-th diagonal
zeroed. For arrays with ``ndim`` exceeding 2, `triu` will apply to the
final two axes.
Please refer to the documentation for `tril` for further details.
See Also
--------
tril : lower triangle of an array
Examples
--------
>>> np.triu([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], -1)
array([[ 1, 2, 3],
[ 4, 5, 6],
[ 0, 8, 9],
[ 0, 0, 12]])
>>> np.triu(np.arange(3*4*5).reshape(3, 4, 5))
array([[[ 0, 1, 2, 3, 4],
[ 0, 6, 7, 8, 9],
[ 0, 0, 12, 13, 14],
[ 0, 0, 0, 18, 19]],
[[20, 21, 22, 23, 24],
[ 0, 26, 27, 28, 29],
[ 0, 0, 32, 33, 34],
[ 0, 0, 0, 38, 39]],
[[40, 41, 42, 43, 44],
[ 0, 46, 47, 48, 49],
[ 0, 0, 52, 53, 54],
[ 0, 0, 0, 58, 59]]])
Eliminación Gaussiana#
def gauelim(inA, inbs):
A = np.copy(inA)
bs = np.copy(inbs)
n = bs.size
for j in range(n - 1):
for i in range(j + 1, n):
coeff = A[i, j] / A[j, j]
A[i, j:] -= coeff * A[j, j:]
bs[i] -= coeff * bs[j]
xs = backsub(A, bs)
return xs
A = np.array([[1, 3, 8], [-2, np.pi, 0], [2, -1, 1]])
A
array([[ 1. , 3. , 8. ],
[-2. , 3.14159265, 0. ],
[ 2. , -1. , 1. ]])
b = np.array([1, 2, 3])
b
array([1, 2, 3])
x_solution = gauelim(A, b)
x_solution
array([ 3.6887277 , 2.98493676, -1.45544225])
A @ x_solution
array([1. , 2. , 2.93707639])
b
array([1, 2, 3])
x_numpy_solution = la.solve(A, b)
x_numpy_solution
array([ 3.75167348, 3.02500929, -1.47833767])
Descomposición LU#
def ludec(A):
n = A.shape[0]
U = np.copy(A)
L = np.identity(n)
for j in range(n - 1):
for i in range(j + 1, n):
coeff = U[i, j] / U[j, j]
U[i, j:] -= coeff * U[j, j:]
L[i, j] = coeff
return L, U
Resolución de un sistema lineal vía LU#
def lusolve(A, bs):
L, U = ludec(A)
ys = forsub(L, bs)
xs = backsub(U, ys)
return xs
x_lu = lusolve(A, b)
x_lu
array([ 3.75167348, 3.02500929, -1.47833767])
%time
y = [i for i in range(1000000)]
CPU times: user 2 μs, sys: 0 ns, total: 2 μs
Wall time: 5.48 μs