SciPy科学计算

一、SciPy中的常数与特殊函数

1.SciPy的constants模块

from scipy import constants as C
print("圆周率：",C.pi)   #圆周率
print("黄金比例：",C.golden)   #黄金比例
print("真空中的光速：",C.c)   #真空中的光速
print("普朗克常数：",C.h)    #普朗克常数
print("一英里等于多少米：",C.mile)   #一英里等于多少米
print("一英寸等于多少米：",C.inch)   #一英寸等于多少米
print("一度等于多少弧度：",C.degree)   #一度等于多少弧度
print("一分钟等于多少秒：",C.minute)  #一分钟等于多少秒
print("标准重力加速度：",C.g)   #标准重力加速度

2.SciPy的special模块

from scipy import special as S
print("立方根：",S.cbrt(8)) #立方根
print("10**3：",S.exp10(3)) #10**3
print("正弦函数：",S.sindg(90)) #正弦函数，参数为角度
print("四舍五入函数3.1：",S.round(3.1)) #四舍五入函数
print("四舍五入函数3.5：",S.round(3.5))
print("四舍五入函数3.499：",S.round(3.499))
print("从5个中任选3个的组合数：",S.comb(5,3)) #从5个中任选3个的组合数
print("排列数：",S.perm(5,3)) #排列数
print("gamma函数：",S.gamma(4)) #gamma函数
print("beta函数：",S.beta(10,200)) #beta函数
print("sinc函数：",S.sinc(0)) #sinc函数

二、SciPy中的线性代数基本运算

1.基本的矩阵计算

#矩阵的创建及其简单运算
from scipy import linalg 
import numpy as np 
A = np.mat('[1,2;3,4]')
print('A矩阵为：\n',A)
print('A的转置矩阵为：\n',A.T)
print('A的逆矩阵为：\n',A.I)

#数组的创建及其简单运算
M = np.array([[1,2],[3,4]])
print('M矩阵为：\n',M)
print('M的转置矩阵为：\n',M.T)
print('M的逆矩阵为：\n',linalg.inv(M))

2.线性方程组求解

方程：x+3y+5z=10，2x+5y-z=6，2x+4y+7z=4

from scipy import linalg 
import numpy as np 
a = np.array([[1, 3, 5], [2, 5, -1], [2, 4, 7]])
b = np.array([10, 6, 4]) 
x = linalg.solve(a, b) 
print(x)

3.行列式的计算

M = np.array([[1,2],[3,4]])
linalg.det(M) 

结果为：-2.0

4.范数

矩阵的常用范数及其含义

矩阵范数	含义
Frobenius	矩阵所有元素平方和的平方根
L范数	矩阵每列元素和的最大值
-1范数	矩阵每列元素和的最小值
2范数	矩阵的最大奇异值
-2范数	矩阵的最小奇异值
正无穷范数	矩阵每行元素和的最大值
负无穷范数	矩阵每行元素和的最小值

#范数的运算
M = np.array([[1,2],[3,4]])
print('M矩阵为：\n',M)
print('M矩阵的L范数为：\n',linalg.norm(M,1))
print('M矩阵的2范数为：\n',linalg.norm(M,2))
print('M矩阵的正无穷范数为：\n',linalg.norm(M,np.inf))

5.特征值分解

#特征值分解
A = np.array([[1,2],[3,4]])
l,v = linalg.eig(A)
print(l)
print(v)

三、数据拟合

1.多项式拟合

import matplotlib.pyplot as plt
x = np.linspace(-5,5,20)
y = 3.5*x+2.1
y_noise = y+np.random.randn(20)*2
coeff = np.polyfit(x,y_noise,1)
plt.plot(x,y_noise,'x',x,coeff[0]*x+coeff[1])
plt.show()

#分别使用1阶、3阶和9阶多项式对数据进行拟合
x = np.linspace(0,np.pi*2)
y = np.sin(x)
f1 = np.poly1d(np.polyfit(x,y,1))
f3 = np.poly1d(np.polyfit(x,y,3))
f9 = np.poly1d(np.polyfit(x,y,9))
t = np.linspace(-3*np.pi,np.pi*3,80)
plt.plot(x,y,'x',t,f1(t),':',t,f3(t),'--',t,f9(t),'-.')
plt.axis([0,3*np.pi,-1.5,1.5])
plt.legend(['data',r'$n = 1$',r'$n = 3$',r'$n = 9$'])
plt.show()

2.最小二乘拟合

#最小二乘估计
from scipy import optimize
def myfunc(x,a,b,w,t):
    return a*np.exp(-b*np.sin(w*x+t))
x = np.linspace(0,2*np.pi)
par = [3,2,1.25,np.pi/4]
y = myfunc(x,*par)
y_noise = y+0.8*np.random.randn(len(y))
def errf(p,x,y):
    return y-myfunc(x,*p)
c,rv = optimize.leastsq(errf,[1,1,1,1],args = (x,y_noise))
#c返回找到的最小二乘估计
plt.plot(x,y_noise,'x',x,y,x,myfunc(x,*c),':')
plt.legend(['data','actual','leastsq'])
plt.show()

四、SciPy中的图像处理

1.图像平滑

import numpy as np
from scipy import ndimage
from scipy import misc
import matplotlib.pyplot as plt
%matplotlib inline
image = misc.face()
aa = plt.subplot(1,3,1)
plt.title("title")
plt.imshow(image)
plt.axis('off')
plt.subplot(1,3,2)
plt.title("medi_filter")
filter = ndimage.median_filter(image,size=10)
#使用SciPy的中值滤波处理图片
plt.imshow(filter)
plt.axis('off')
plt.subplot(1,3,3)
plt.title("gausfilter")
blurred_face = ndimage.gaussian_filter(image, sigma = 7)#高斯滤波
plt.imshow(blurred_face)
plt.axis('off') 

2.图像的旋转和锐化

image = misc.face() 
plt.subplot(131)
plt.title("title")
plt.imshow(image)
plt.axis('off')
plt.subplot(132)
rotate = ndimage.rotate(image,60)
plt.title("rotate")
plt.imshow(rotate)
plt.axis('off') #边缘检测
plt.subplot(133)
prewitt = ndimage.prewitt(image)
plt.title("prewitt")
plt.imshow(prewitt)
plt.axis('off')