Major libraries that have been used in this tutorial¶
In [1]:
import math
import numpy as np
import matplotlib.pyplot as plt
import scipy
import sympy
Introduction to Mathematical Computations with Python¶
Welcome to Numerical Computations!
In this simple tutorial, we will learn about python3's basic commands and methods that we will use them regularly for numerical computations. Consider that we use Python 3.7
Note¶
This tutorial has been created by $\href{https://hnxj.github.io/}{Hamed \: Nejat}$ for Numerical computations course of $\href{https://ee.sharif.edu}{Electrical \: Engineering \: Department \: of \: Sharif \: University \: of \: Technology}$, with major help of official documents of Numpy, Scipy, Matplotlib and Sympy. special thanks to all other python developers, this tutorial is free to use.
Anaconda with Python3.7 Download links:¶
Windows (x64/x86) (~490MB) : https://repo.anaconda.com/archive/Anaconda3-2019.07-Windows-x86_64.exe
for MacOS (~650MB) : https://repo.anaconda.com/archive/Anaconda3-2019.07-MacOSX-x86_64.pkg
for Linux (~510MB) : https://repo.anaconda.com/archive/Anaconda3-2019.07-Linux-x86_64.sh
Make sure that JupyterLab is installed by this command on Anaconda Prompt :
-->conda install jupyter
Virtual environments:¶
In anaconda prompt, write:
-->conda create -n env_name python=3.7
Then everytime you want to code for this course, write:
-->conda activate evn_name
And when you dont want to use this env:
-->conda deactivate env_name
for anaconda users:
-->conda install jupyter
for main python users:
-->pip install
Run Jupyter notebook:¶
Open Anaconda prompt, one of the ways is to search "Anaconda prompt" in windows search.
After opening of Anaconda Prompt, Enter:
-->jupyter lab
Then, you are in jupyter environment. click on file tab, then click on new notebook.
Run Spyder scientific IDE¶
On cmd/terminal, install spyder using this commands:
-->conda install spyder
To run:
-->spyder
Contact¶
For more questions about tutorial class or homeworks:
contact us : hamednejat7@gmail.com
Input/Output (I/O)¶
In [2]:
print("You will never walk alone") # Raw string
print("There is" + " a " + "GoldenSky :)" # Multi strings
" Walk on", "trough the rain") # multi line
In [3]:
a = int(input("Enter int number"))
b = float(input("enter float"))
s = str(input("Enter string"))
print(a, b, s)
In [4]:
print(6) # Integers
print(-7) # Integers
print(18.0) # Floats
print("18.0 " + str(-32.4353)) # Numbers to string
print(float("18.34") + 32.4353) # Strings to number
print(float("-18.34") + 32.4353) # Signed strings to number, Sign is included in number while it was converted from string to number!
In [6]:
print("First number: {}, Second number: {} and so on ...".format(4.01, 3.98))
print ('%.2f'%(11.1725612)) # 2-digit precision, rounded down because 11.1725612 < 11.175
print ('%.2f'%(11.1775612)) # 2-digit presicion rounded up because 11.1775612 > 11.175
print('%.5f'%(float("18.34") + 32.4353)) # Strings to number, 5-digit precision
print('%.3f'%(float("-18.34") + 32.4353)) # Signed strings to number, 3-digit precision
print('%f %d %s'%(11.13, 11.13, 11.13)) # printing float, integer and string in same command.
print('%f'%(1323.12312312312)) # NOTE!: default float precision is 6-digits
In [7]:
a = 4 # integer
b = -5.12 # float
c = "Estereeng" # string
# char is almost same with str in python, except char array in classes
d = chr(113) # <C: char> = chr(<I: integer>) method returns the character that has ascii=I
e = ord('q') # <I: integer> = ord(<C: char>) method returns the ascii that char=C
print("\nType of a is :", type(a), "a :", a,
"\nType of b is :", type(b), "b :", b,
"\nType of c is :", type(c), "c :", c,
"\nType of d is :", type(d), "d :", d,
"\nType of e is :", type(e), "d :", e)
1.2: Dictionary¶
In [9]:
e = dict() # dict :: e[<Key>] = <Val>
e["asd"] = 12 # assigning num to string key
e[1] = a # integer to integer
e[b] = 'L' # char (string) to float
print("\nType of e is :", type(e), "\n e :", e)
print("e[\"asd\"] is : ", e["asd"])
0.1.3: List¶
In [10]:
f = list() # alternate for list / f = [1, 1.1, list]
f.append(1) # appending an integer
f.append(1.1) # appending a float
f.append("paganini") # appending a string
f.append(e) # appending a dictionary!
print("\nType of e is :", type(f), "e :", f)
f.pop(2) # pop 2nd element
print("\nType of e is :", type(f), "e :", f)
f.pop() # pop last element
print("\nType of e is :", type(f), "e :", f)
1.3: Set¶
In [12]:
f = set() # alternate for list / f = [1, 1.1, list]
f.add(1) # add an integer to the set
f.add(1.1) # appending a float
f.add(2)
f.add(-1)
f.add("vivaldi") # appending a string
f.add("antonio")
# f.add(e) # appending a dictionary!
print("\nType of e is :", type(f), "e :", f)
f.remove("antonio") # remove one existing element
try: # prevents unexcepted errors and interrupts in program
f.remove("bach") # key is not in the set
except KeyError:
print("\nKey is not in the set")
except:
print("Unknown error")
print("\nType of e is :", type(f), "e :", f)
f.pop() # pop last element
print("\nType of e is :", type(f), "e :", f)
1.4: Computation¶
In this subsection, we will do some maths computations and evaluate errors.
In [14]:
import math
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
a = 12.250
b = 27.601
c = a * b
d = a / b
print(a, b, c, d)
In [15]:
a = 3.131
b = 13.313
c = a ** b # a power b
d = np.sqrt(b) # radical b
print(a, b, c, d)
In [16]:
a = 3.131
b = np.sin(a) # numpy sinus; is radian based
c = np.cos(a) # numpy cosinus; is radian based
d = np.exp(c) # exponential
print(a, b, c, d)
In [17]:
a = -3.131 + 4.023j # complex number
b = np.sin(a) # numpy sinus; is radian based
c = np.cos(a) # numpy cosinus; is radian based
d = np.exp(c) # exponential
print(a, b, c, d)
Some Random Numbers¶
In [24]:
a = np.random.rand(5) # 10 random numbers uniformly between 0-1
b = np.random.randn(1000) # 1000 Gaussian standard distributed random numbers with mean = 0, std = 1
c = np.random.randint(-3, 14, 10) # 10 random integers between -3 and 13
d = np.random.gamma(2.5, 1, [2, 2]) # beta, gumbel, ... other distributions are available
print(a)
print(np.mean(b))
print(np.var(b))
print(c)
print(d)
1.5 Stats¶
In [26]:
import numpy as np
interval_first = 0
interval_last = 10
step = 2
amount = 33
e = np.arange(interval_first, interval_last, step) # equivalent to matlab's i:j with step
f = np.linspace(interval_first, interval_last, amount) # Splits line space between 1 and 10 into 25 numbers
g = np.mean(e)
h = np.var(f)
i = np.std(f)
j = np.max(e)
k = np.min(f)
l = np.median(e)
print(e)
print(f)
print(len(e)) # returns length of structure
print(len(f)) # returns length of structure
print(g, h, i)
print(j, k, l)
2: Loop and Condition¶
In [3]:
import numpy as np
first = -4
last = 6
step = 2
cnt = 0.5
np_array1 = np.array([1, 0.7, -4, 2.0 + 2.0j, -0.5j])
np_array2 = np.array([[1, 2], [0.5j, -2], [4.0, 0.3j]])
For loops¶
In [5]:
for i in range(first, last, step): # ranged loop
print(i**2.5)
for j in np_array1: # structured loop
print(j)
for i, j in np_array2: # structured double loop
print(i, j)
In [8]:
for i in range(10):
print(i)
While loops¶
In [9]:
while cnt < 10:
cnt *= 2
print(cnt)
If conditions¶
In [10]:
if step > 10:
pass # the command that does nothing!
# the code
elif step == 10: # Else If equivalent
pass
# the code
elif step == 9:
quit() # Exit from the code or program !
else: # otherwise
print("Nothing else matters")
3: Functions¶
In [11]:
import numpy as np
# Basic recursive function
def gcd(a, b): # Greatest Common Divisor (GCD)
if b > a:
(a, b) (b, a) # nice swap in Python!
if b == 0:
return a
return gcd(b, a%b)
R = 14
a = -1
f = lambda tau : 14 - ((1.0 - np.exp(-tau))/(1.0 - np.exp(-a*tau)))
print(gcd(2436, 504))
print(f(12))
# IMPORTANT NOTE: Python has recursion limit, default recursion limit is 1000.
# you can change it by setrecursionlimit(<New Limit>)
import sys
sys.setrecursionlimit(10000)
In [13]:
# Void function
def say_hello():
print("Hello World :)")
return # void return
say_hello()
In [16]:
# lambda function
hemographic1 = lambda x : (2 * x + 4) / (3 * x - 1)
print(hemographic1(4))
4: Classes¶
Calm down, Its just a class! The other parts like plots will be discussed in the next tutorial
In [33]:
from matplotlib import pyplot as plt
import numpy
class painter:
def __init__(self, boom_width, boom_height):
self.boom_width = boom_width
self.boom_height = boom_height
self.fig, self.ax = plt.subplots()
def sinusy(self, x):
self.ax.plot(x, np.sin(x))
return
def expony(self, x):
self.ax.plot(x, np.exp(x))
return
def free_draw(self, x, y):
self.ax.plot(x, y)
print("Mahmud, 10 saale az Salafchegan")
return
In [34]:
msr = painter(100, 100)
msr.sinusy(np.linspace(1, 7, 30))
msr.expony(np.linspace(-3, 3, 30))
msr.free_draw(np.random.rand(30), np.linspace(10, 20, 30))
5.Plotting with matplotlib¶
In [31]:
from matplotlib import pyplot as plt
x = np.linspace(-3, 3, 100)
y1 = np.exp(x)
y2 = np.cos(x)
y3 = np.random.randn(100)
plt.plot(x, y3)
plt.show()
# print(dir(plt))
Initialize figures
In [32]:
fig, ax = plt.subplots(nrows=2, ncols=2) # number of subplots
fig.set_figwidth(19) # figure width in cm
fig.set_figheight(11) # figure height in cm
ax[0, 0].plot(x, y2, label='Cos', color='red') # plotting with specs
ax[0, 0].plot(x, y3, label='Rnd') # colors, scatters and ...
ax[0, 0].plot(x, y1, label='Sin')
ax[0, 1].scatter(x, y2, label='Cos', color='cyan')
ax[1, 0].plot(x, y3, label='Rnd')
ax[1, 1].plot(x, y2, label='Cos')
ax[1, 1].plot(x, y3, label='Rnd')
ax[0, 1].annotate('Sin', (5, 0.5),
xytext=(0.8, 0.9), textcoords='axes fraction',
arrowprops=dict(facecolor='black', shrink=0.05),
fontsize=16,
horizontalalignment='right', verticalalignment='top') # annotation cursor
ax[0, 0].set_title("The title !")
ax[0, 0].set_xlabel("X description")
ax[0, 0].set_ylabel("Y description")
ax[1, 0].set_title("The title !")
ax[1, 0].set_xlabel("X description")
ax[1, 0].set_ylabel("Y description")
ax[0, 1].set_title("The scatter")
ax[0, 1].set_xlabel("X description")
ax[1, 1].set_ylabel("Y description")
ax[0, 0].legend()
fig.show()
In [6]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
a = 0.5
R = 1.6
func = lambda tau : R - ((1.0 - np.exp(-tau))/(1.0 - np.exp(-a*tau)))
# Plot it
tau = np.linspace(-0.5, 1.5, 200)
fig, ax = plt.subplots(1, 1)
fig.set_figwidth(19)
fig.set_figheight(11)
ax.plot(tau, func(tau))
ax.set_xlabel("tau")
ax.set_ylabel("expression value")
ax.grid(True)
fig.show()
# Use the numerical solver to find the roots
tau_initial_guess = 0.5
tau_solution = fsolve(func, tau_initial_guess)
print("The solution is tau = %f"%(tau_solution))
print("at which the value of the expression is %f"%(func(tau_solution)))
In [2]:
from matplotlib.pyplot import *
import numpy as np
x = np.linspace(0, 10, 100)
y = np.sin(x**2)
plot(x, y)
show()
7. Numpy and scipy Interpolation, extrapolation and fitting¶
In [10]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import lagrange
# data initialize as numpy array
x = np.linspace(-4, 8, 25) # X
y = np.random.randint(-5, 5, 25) # Y
xp = np.linspace(-5, 9, 300) # Test interval for our polynomial
# Interpolation with Lagrange's method
p1 = lagrange(x, y)
# Fitting with 6th-degree polynomial
pmul2 = np.polyfit(x, y, deg=2) # getting polynomial multiples
p2 = np.poly1d(pmul2) # getting polynomial function
# a = np.array([1, 4, 5])
# f = np.poly1d(a)
# Plot em
fig, ax = plt.subplots(2, 1)
fig.set_figwidth(18)
fig.set_figheight(20)
# Interpolated plotting
ax[0].scatter(x, y)
ax[0].plot(xp, p1(xp))
ax[0].set_xlabel("Time")
ax[0].set_ylabel("Value")
ax[0].grid()
ax[0].set_title("Interpolated")
ax[0].set_ylim([-30, 30])
# Changes plotting
ax[1].scatter(x, y)
ax[1].plot(xp, p2(xp))
ax[1].set_xlabel("Time")
ax[1].set_ylabel("Value")
ax[1].grid()
ax[1].set_title("Fitted with 6th-degree")
ax[0].set_ylim([-10, 10])
Out[10]:
7. Derivative¶
Welcome to sympy!
sympy gives us symbolic abilities, like differentiation, integration, transforms and solving DEs.
Key commands and methods:¶
x, y = symbols('x y') # -> define x, y as symbols
f1 = x**2 + x + sin(x) # -> define a symbolic function
fd1 = diff(f1, x, 2) # -> compute second order symbolic derivative of f1 based on x
f2 = x**2 + arcsin(y) # -> define multivariable function
fd2 = diff(f1, y, 1) # -> partial differentiation
7.1 Symbolic differentiation¶
Let see what is in sympy library :
In [1]:
import sympy
print(dir(sympy))
In [2]:
from sympy import * # imports all required functions, methods and types
x, y, z, t, snt = symbols('x y z t snt') # symbols
f1 = exp(x*y + t**(x*y)) + cos(t*2) + sin(y + x + t)
fd1 = diff(f1, t, 1) # first-order derivative based on x
fd1 = diff(fd1, x, 1)
fd1 = diff(fd1, y, 2)
print("\nFirst function:")
print(fd1)
print(fd1.subs(x, 3)) # string output
print(fd1.subs(x, 3.)) # float
f2 = sin(x + y) + y**3 + x*y**4
fd2 = diff(f2, y, 2) # second-order derivative based on y
print("\nSecond function:")
print(fd2)
print(fd2.subs(x, 2.).subs(y, 1.)) # getting values from the function
7.2 Numerical differentiation¶
As easy as easy. there are two different ways for numerical differentiation:
In [84]:
import numpy as np
x = np.random.randint(-3, 3, 10)
y1 = np.gradient(x, 1) # the interpolated gradient
y2 = x[0:len(x)-1] - x[1:len(x)] # normal gradient
print(x)
print(y1)
print(y2)
In [2]:
from sympy import * # imports all required functions, methods and types
x, y, z, t = symbols('x y z t') # symbols
f1 = (x**2)*cos(x) + exp(x*2)
I1 = integrate(f1, x) # symbolic integral based on x
print("\nFirst function integral:")
print(I1)
print(I1.subs(x, 3)) # string output
print(I1.subs(x, 3.) - I1.subs(x, 1.)) # float
f2 = sin(x + y) + y**3 + x*y**4
I2 = integrate(f2, y) # symbolic integral based on y
print("\nSecond function:")
print(I2)
print(I2.subs(x, 2.).subs(y, 1.)) # getting values from the function
In [102]:
from sympy.integrals import laplace_transform
from sympy.abc import t, s, a
F1 = t**a + sin(t) + exp(3*a*t)
L1 = laplace_transform(F1, t, s)
print(L1)
8.2.2 Inverse Laplace transform¶
In [104]:
from sympy.integrals.transforms import inverse_laplace_transform
from sympy import exp, Symbol
from sympy.abc import s, t
a = Symbol('a', positive=True)
L2 = exp(-a*s)/s**2
F2 = inverse_laplace_transform(L2, s, t)
print(F2)
8.2.3 Mellin transform¶
In [106]:
from sympy.integrals.transforms import mellin_transform
from sympy import exp
from sympy.abc import x, s
F1 = exp(-x) + sin(x)
M1 = mellin_transform(F1, x, s)
print(M1)
8.2.4 Fourier and I-Fourier transform¶
In [115]:
from sympy import fourier_transform, exp, inverse_fourier_transform
from sympy.abc import x, k
F1 = exp(-x*2) * Heaviside(x)
Q1 = fourier_transform(F1, x, k)
Q2 = fourier_transform(F2, x, k, noconds=False)
print(Q1)
print(Q2)
In [9]:
from sympy import Function, dsolve, Eq, Derivative, sin, cos, symbols
from sympy.abc import x
f = Function('f')(x) # Defining F as function of x
fp = Derivative(f, x) # F prime (1st derivative)
fz = Derivative(f, x, x) # F zegond (2nd derivative)
fy = Derivative(fz, x)
deq1 = fz + 9*f - fy
p1 = dsolve(deq1, f)
print(p1)
9.1 ODE¶
Visit [https://docs.sympy.org/latest/modules/solvers/ode.html], ODE has so many properties
9.2 PDE¶
9.3 Euler's method for numerical ODE approximation¶
In [142]:
from matplotlib import pyplot as plt
def f(x, y):
return (x + y - x * y) + (x**2) -(y**3)
def euler_ode(x0, y, h, x):
sol = []
xs = []
while x0 < x:
y = y + h * f(x0, y)
x0 = x0 + h
sol.append(y)
xs.append(x0)
return xs, sol
x0 = 0
y0 = 1
h = 0.025
# Value of x at which we need approximation
x = 10
x, y = euler_ode(x0, y0, h, x)
fig, ax = plt.subplots()
ax.plot(x, y)
9.4 Runge-Kutta's method for DE¶
Imagine we have this Differential equation:
$$y'' - xy' - y + x^2 = 0$$To solve this equation as numerical DE, we have to rewrite this equation as Runge-Kutta equations:
$$y'' = xy' + y - x^2, f_2(x) = xy' + y - x^2, f_1(x) = y'$$$$y'(x_0) = 2, y(x_0) = 0$$In [72]:
from matplotlib import pyplot as plt
import numpy as np
def runge_kutta4_2(x0, y10, y20, x, h, F1, F2): # for 2nd order
n = (int)((x - x0)/h)
y1 = y10
y2 = y20
lgx = []
lgy = []
for i in range(1, n+1):
lgx.append(x)
lgy.append(y1)
k11 = h * F1(x, y1, y2)
k21 = h * F2(x, y1, y2)
k12 = h * F1(x + 0.5*h, y1 + 0.5*k11, y2 + 0.5*k21)
k22 = h * F2(x + 0.5*h, y1 + 0.5*k11, y2 + 0.5*k21)
k13 = h * F1(x + 0.5*h, y1 + 0.5*k12, y2 + 0.5*k22)
k23 = h * F2(x + 0.5*h, y1 + 0.5*k12, y2 + 0.5*k22)
k14 = h * F1(x + h, y1 + k13, y2 + k23)
k24 = h * F2(x + h, y1 + k13, y2 + k23)
y1 += (k11 + 2*(k12 + k13) + k14)/6
y2 += (k21 + 2*(k22 + k23) + k24)/6
x += h
return lgx, lgy
def runge_kutta4_1(x0, y0, x, h, F): # for 1st order
n = (int)((x - x0)/h)
y = y0
lgx = []
lgy = []
for i in range(1, n+1):
lgx.append(x0)
lgy.append(y)
k1 = h * F(x0, y)
k2 = h * F(x0 + 0.5 * h, y + 0.5 * k1)
k3 = h * F(x0 + 0.5 * h, y + 0.5 * k2)
k4 = h * F(x0 + h, y + k3)
y = y + (1.0 / 6.0)*(k1 + 2 * k2 + 2 * k3 + k4)
x0 = x0 + h
return lgx, lgy
def f1(x, y1, y2):
dy1dx = y2
return dy1dx
def f2(x, y1, y2):
dy2dx = -y1 * x**2 + y2 * np.sin(x)
return dy2dx
def f3(x, y):
dydx = -y + x
return dydx
x0 = -3
y = 3
yp = 0
x = 3
h = 0.01
# print('The value of y at x is:', runge_kutta4_1(x0, y, x, h, f3))
# print('The value of y at x is:', runge_kutta4_2(x0, y, yp, x, h, f1, f2))
x1, y1 = runge_kutta4_1(x0, y, x, h, f3)
x0 = -20
y = 3
yp = 0
x = 0
h = 0.01
x2, y2 = runge_kutta4_2(x0, y, yp, x, h, f1, f2)
# Plot em
fig, ax = plt.subplots(2, 1)
fig.set_figwidth(18)
fig.set_figheight(20)
# Interpolated plotting
ax[0].plot(x1, y1)
ax[0].set_xlabel("x")
ax[0].set_ylabel("y")
ax[0].grid()
ax[0].set_title("Linear DE")
# ax[0].set_ylim([-30, 30])
# Changes plotting
ax[1].plot(x2, y2)
ax[1].set_xlabel("x")
ax[1].set_ylabel("y")
ax[1].grid()
ax[1].set_title("Non-linear DE")
# ax[0].set_ylim([-10, 10])
fig.show()
In [32]:
import numpy as np
x1 = np.random.randint(0, 10, [3, 3])
x2 = np.random.randint(1, 24, [3, 3])
x3 = np.random.randint(3, 10, [3, 7])
x4 = np.random.randint(4, 11, [7, 5])
x5 = np.random.randint(10, 15, [5, 1])
x6 = np.random.randint(-1, 7, [1, 2])
In [33]:
# Matrix addition
y = x1 + x2 #(elementwise addition)
print("x1:\n", x1, "\nx2:\n", x2, "\nx3:\n", y)
In [34]:
# Matrix subtraction
y = x1 - x2 #(elementwise subtraction)
print("x1:\n", x1, "\nx2:\n", x2, "\nx3:\n", y)
In [36]:
# Matrix elementwise division
y = x1 / x2 #(elementwise division)
print("x1:\n", x1, "\nx2:\n", x2, "\nx3:\n", y)
In [61]:
# Matrix elementwise multiply (dot)
y = x1 * x2 # (dot)
y = np.dot(x1, x2) # equivalent :)
print("x1:\n", x1, "\nx2:\n", x2, "\nx3:\n", y)
In [42]:
# Matrix multiply
y = np.matmul(x3, x4) #(matrix cross (multiply))
print("x1:\n", x3, "\nx2:\n", x4, "\nx3:\n", y)
In [57]:
xt = np.transpose(x4) # transpose!
print(xt, "\n\n", x4)
10.1 Numpy.linalg¶
this subset of numpy has basic Linear algebra tools
In [1]:
import numpy as np
x1 = np.random.randint(1, 10, [3, 3])
x2 = np.random.randint(1, 24, [3, 3])
x3 = np.random.randint(3, 10, [3, 7])
x4 = np.random.randint(4, 11, [7, 5])
x5 = np.random.randint(10, 15, [5, 1])
x6 = np.random.randint(-1, 7, [1, 2])
v1 = np.random.randint(-1, 3, [1, 3])
v2 = np.random.randint(1, 19, [1, 3])
v3 = np.random.randint(-4, 4, [3, 3])
print(dir(np.linalg))
xinv = np.linalg.inv(x1) # inverse matrix
xdet = np.linalg.det(x1) # determinant
xnrm = np.linalg.norm(x1) # norm
xsvd = np.linalg.svd(x1) # singular value decomposition
xqrd = np.linalg.qr(x1) # compute the qr factorization a matrix
xeig = np.linalg.eig(x1) # compute eigenvalues and correponding eigenvactors
xpow = np.linalg.matrix_power(x1, 3) # x1^3
xpow = np.linalg.matrix_power(x1, -1) # x1^-1 (inverse)
xpow = np.linalg.cond(x1)
# xchl = np.linalg.cholesky(x1) # matrix must be positive and definite
10.2 Scipy.linalg¶
Scipy has more tools, including np.linalg and many other tools:
example of L-U decomposition:
In [3]:
import numpy as np
import scipy.linalg as la
np.set_printoptions(suppress=True)
A = np.array([[1,3,4],[2,1,3],[4,1,2]])
print(A)
P, L, U = la.lu(A)
print(np.dot(P.T, A))
print
print(np.dot(L, U))
print(P)
print(L)
print(U)
print(dir(la))
For full documentation, visit scipy's np.linalg simple documentation:
https://docs.scipy.org/doc/numpy/reference/routines.linalg.html