欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识！
您现在的位置是： 首页

Python 科学计算学习二：SciPy-数值计算库

程序员文章站 2022-07-12 22:19:01
...

SciPy函数库在NumPy库的基础上增加了众多的数学、科学以及工程计算中常用的库函数。例如线性 代数、常微分方程数值求解、信号处理、图像处理、稀疏矩阵等等。
注：python版本不同，语法会略微有所差异。
我所使用的是Python 3.6.5 |Anaconda, Inc.|

1 最小二乘拟合

如果能找到一组 S(P)=∑i=1m[yi−f(xi,p)]2 中的p 值使得S函数最小，这种算法被称之为最小二乘拟合。
scipy中的子函数库optimize已经提供了实现最小二乘拟合算法的函数leastsq。下面是用leastsq进行 数据拟合的一个例子：
下面是用leastsq进行 数据拟合的一个例子：

# -*- coding: utf-8 -*-
import numpy as np
from scipy.optimize import leastsq
import pylab as pl
from pylab import mpl  
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False #使图表中负号正常输出
def func(x, p):
    A, k, theta = p
    return A*np.sin(2*np.pi*k*x+theta)
def residuals(p, y, x):
    return y - func(x, p)
x = np.linspace(0, -2*np.pi, 100)
A, k, theta = 10, 0.34, np.pi/6 # 真实数据的函数参数
y0 = func(x, [A, k, theta]) # 真实数据
y1 = y0 + 2 * np.random.randn(len(x)) # 加入噪声之后的实验数据
p0 = [7, 0.2, 0] # 第一次猜测的函数拟合参数
# 调用leastsq进行数据拟合
# residuals为计算误差的函数
# p0为拟合参数的初始值
# args为需要拟合的实验数据
plsq = leastsq(residuals, p0, args=(y1, x))
print (u"真实参数:", [A, k, theta])
print (u"拟合参数", plsq[0]) # 实验数据拟合后的参数
pl.plot(x, y0, label=u'真实数据')  #设置坐标轴
pl.plot(x, y1, label=u"带噪声的实验数据")
pl.plot(x, func(x, plsq[0]), label=u"拟合数据")
pl.legend()
pl.show()

Python 科学计算学习二：SciPy-数值计算库
实验结果：

>>> 真实参数: [10, 0.34000000000000002, 0.52359877559829882] 
>>> 拟合参数 [-9.84152775 0.33829767 -2.68899335]

Python 科学计算学习二：SciPy-数值计算库
我们看到拟合参数虽然和真实参数完全不同，但是由于正弦函数具有周期性，实际上拟合参数得到的 函数和真实参数对应的函数是一致的。

2 函数最小值

optimize库提供了几个求函数最小值的算法：fmin, fmin_powell, fmin_cg, fmin_bfgs。

3 非线性方程组求解

optimize库中的fsolve函数可以用来对非线性方程组进行求解。它的基本调用形式如下：

fsolve(func, x0)

func(x)是计算方程组误差的函数。 func返回将x代入方程组之后得到的误差；x0为未知数矢量的初始值。
实例：求解如下方程组的解
5∗x1+3=0
4∗x02−2∗sin(x1∗x2)=0
x1∗x2−1.5=0
代码：

# -*- coding: utf-8 -*-
from scipy.optimize import fsolve
from math import sin,cos
def f(x):
    x0 = float(x[0])
    x1 = float(x[1])
    x2 = float(x[2])
    return [
        5*x1+3,
        4*x0*x0 - 2*sin(x1*x2),
        x1*x2 - 1.5
    ]   
result = fsolve(f, [1,1,1])
print (result)
print (f(result))

实验结果：
Python 科学计算学习二：SciPy-数值计算库

4 B-Spline样条曲线

interpolate库提供了许多对数据进行插值运算的函数。下面是使用直线和B-Spline对正弦波上的点进行插值的例子。

# -*- coding: utf-8 -*-
import numpy as np
import pylab as pl
from scipy import interpolate
from pylab import mpl  
mpl.rcParams['font.sans-serif'] = ['SimHei'] #图表中的文字正常输出
mpl.rcParams['axes.unicode_minus'] = False #使图表中负号正常输出
x = np.linspace(0, 2*np.pi+np.pi/4, 10)
y = np.sin(x)
x_new = np.linspace(0, 2*np.pi+np.pi/4, 100)
f_linear = interpolate.interp1d(x, y)
tck = interpolate.splrep(x, y)
y_bspline = interpolate.splev(x_new, tck)

pl.plot(x, y, "o", label=u"原始数据")  #设置坐标轴
pl.plot(x_new, f_linear(x_new), label=u"线性插值")
pl.plot(x_new, y_bspline, label=u"B-spline插值")
pl.legend()
pl.show()

实验结果：
Python 科学计算学习二：SciPy-数值计算库
在这段程序中，通过interp1d函数直接得到一个新的线性插值函数。而B-Spline插值运算需要先使用 splrep函数计算出B-Spline曲线的参数，然后将参数传递给splev函数计算出各个取样点的插值结果。

5 数值积分

数值积分是对定积分的数值求解，多重定积分的求值可以通过多次调用quad函数实现，为了调用方便，integrate库提供了dblquad函数进行二重定积分，tplquad函数进行三重定积分。
下面以计算单位半球体积为例说明dblquad函数的用法。
单位半球上的点(x,y,z)符合如下方程：

x2+y2+z2=1

可以如下定义通过(x,y)坐标计算球面上点的z值的函数：

def half_sphere(x, y): 
    return (1-x**2-y**2)**0.5

X-Y轴平面与此球体的交线为一个单位圆，因此积分区间为此单位圆，可以考虑为X轴坐标从-1到1进 行积分，而Y轴从 -half_circle(x) 到 half_circle(x) 进行积分，于是可以调用dblquad函数：

>>> integrate.dblquad(half_sphere, -1, 1, 
    lambda x:-half_circle(x), 
    lambda x:half_circle(x)) 
>>> (2.0943951023931988, 2.3252456653390915e-14) >>> np.pi*4/3/2 # 通过球体体积公式计算的半球体积
2.0943951023931953

dblquad函数的调用方式为：

dblquad(func2d, a, b, gfun, hfun)

对于func2d(x,y)函数进行二重积分，其中a,b为变量x的积分区间，而gfun(x)到hfun(x)为变量y的积分区间。

6 解常微分方程组

scipy.integrate库提供了数值积分和常微分方程组求解算法odeint。下面让我们来看看如何用odeint 计算洛仑兹吸引子的轨迹。洛仑兹吸引子由下面的三个微分方程定义：

dxdt=σ(y−x)

dydt=x(ρ−z)−y

dzdt=xy−βz

这三个方程定义了三维空间中各个坐标点上的速度矢量。从某个坐标开始沿着速度矢量进行积分，就 可以计算出无质量点在此空间中的运动轨迹。其中 σ, ρ, β 为三个常数，不同的参数可以计算出不同的 运动轨迹： x(t), y(t), z(t)。 当参数为某些值时，轨迹出现馄饨现象：即微小的初值差别也会显著地影 响运动轨迹。下面是洛仑兹吸引子的轨迹计算和绘制程序：

# -*- coding: utf-8 -*-
from scipy.integrate import odeint
import numpy as np

def lorenz(w, t, p, r, b):
    #给出位置矢量w，和三个参数p, r, b计算出 dx/dt, dy/dt, dz/dt的值
    x, y, z = w
    # 直接与lorenz的计算公式对应
    return np.array([p*(y-x), x*(r-z)-y, x*y-b*z])

t = np.arange(0, 30, 0.01) # 创建时间点
#调用ode对lorenz进行求解, 用两个不同的初始值
track1 = odeint(lorenz, (0.0, 1.00, 0.0), t, args=(10.0, 28.0, 3.0))
track2 = odeint(lorenz, (0.0, 1.01, 0.0), t, args=(10.0, 28.0, 3.0))
#绘图
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

fig = plt.figure()
ax = Axes3D(fig)
ax.plot(track1[:,0], track1[:,1], track1[:,2])
ax.plot(track2[:,0], track2[:,1], track2[:,2])
plt.show()

实验结果：
Python 科学计算学习二：SciPy-数值计算库
我们看到即使初始值只相差0.01，两条运动轨迹也是完全不同的。

7 滤波器设计

scipy.signal库提供了许多信号处理方面的函数。下面利用signal库设计简单滤波器，查看滤波器的频率响应，以及如何使用滤波器对信号进行滤波。
导入signal库:

import scipy.signal as signal

下面的程序设计一个带通IIR滤波器：

# -*- coding: utf-8 -*-
import numpy as np
import scipy.signal as signal
import pylab as pl

b, a = signal.iirdesign([0.2, 0.5], [0.1, 0.6], 2, 40)#滤波器的通带为0.2*f0到0.5*f0，阻带为小于0.1*f0和大于0.6*f0，其中f0为1/2的信号取样频率
#iirdesgin返回的两个数组b和a， 它们分别是IIR滤波器的分子和分母部分的系数。其中a[0]恒等于1
w, h = signal.freqz(b, a)#freqz返回两个数组,其中w是圆频率数组,h是w中的对应频率点的响应，它是一个复数数组，其幅值为滤波器的增益，相角为滤波器的相位特性

power = 20*np.log10(np.clip(np.abs(h), 1e-8, 1e100))
pl.plot(w/np.pi*4000, power) #绘制滤波器的增益特性图，假设取样频率为8kHz
t = np.arange(0, 2, 1/8000.0) #产生2秒钟取样频率为8kHz的取样时间数组
sweep = signal.chirp(t, f0=0, t1 = 2, f1=4000.0)#调用chirp得到2秒钟的频率扫描波形的数据：
out = signal.lfilter(b, a, sweep)#调用lfilter函数计算sweep波形经过带通滤波器之后的结果
out = 20*np.log10(np.abs(out))#将输出波形数据转换为能 量值
index = np.where(np.logical_and(out[1:-1] > out[:-2], out[1:-1] > out[2:]))[0] + 1
pl.plot(t[index]/2.0*4000, out[index] )#将时间转换为对应的频率，绘制所有局部最大点的能量值

下图显示freqz计算的频谱和频率扫描波得到的频率特性（蓝色为freqz计算的滤波器频谱，橙色为频率扫描波测量的滤波器频谱），我们看到其结果是一致的。
Python 科学计算学习二：SciPy-数值计算库

8 用Weave嵌入C语言

Python作为动态语言其功能虽然强大，但是在数值计算方面有一个最大的缺点：速度不够快。在Python级别的循环和计算的速度只有C语言程序的百分之一。SciPy提供了快速调用C++语言程序的方法– Weave。但 Weave 仅适用于 Python 2.x 。
下面是对NumPy的数组求和的例子：

# -*- coding: utf-8 -*-
import scipy.weave as weave
import numpy as np
import time

def my_sum(a):
    n=int(len(a))
    code="""
    int i;
    double counter;
    counter =0;
    for(i=0;i<n;i++){
        counter=counter+a(i);
    }
    return_val=counter;
    """
    err=weave.inline(
        code,['a','n'],
        type_converters=weave.converters.blitz,
        compiler="gcc"
    )
    return err

a = np.arange(0, 10000000, 1.0)
# 先调用一次my_sum，weave会自动对C语言进行编译，此后直接运行编译之后的代码
my_sum(a)

start = time.clock()
for i in xrange(100):
    my_sum(a)  #直接运行编译之后的代码
print ("my_sum:"),(time.clock() - start) / 100.0

start = time.clock()
for i in xrange(100):
    np.sum(a) #numpy中的sum，其实现也是C语言级别
print("np.sum:"),(time.clock() - start) / 100.0

start = time.clock()
sum(a)  #Python内部函数sum通过数组a的迭代接口访问其每个元素，因此速度很慢
print("sum:"), time.clock() - start

weave.inline函数的第一个参数为需要执行的C++语言代码，第二个参数是一个列表，它告诉weave 要把Python中的两个变量a和n传递给C++程序，注意我们用字符串表示变量名。converters.blitz是 一个类型转换器，将numpy的数组类型转换为C++的blitz类。C++程序中的变量a不是一个数组，而 是blitz类的实例，因此它使用a(i)获得其各个元素的值，而不是用a[i]。最后我们通过compiler参数告 诉weave要采用gcc为C++编译器。
运行结果：

my_sum: 0.0294527349146 
np.sum: 0.0527649547638 
sum: 9.11022322669
相关标签： Python SciPy

上一篇： Mybatis 学习标注

下一篇： gulp与webpack对比

推荐阅读