SciPy是一种使用NumPy来做高等数学、信号处理、优化、统计和许多其它科学任务的语言扩展,SciPy函数库在NumPy库的基础上增加了众多的数学、科学以及工程计算中常用的库函数。例如线性代数、常微分方程数值求解、信号处理、图像处理、稀疏矩阵等等。使用时要首先安装合适版本的scipy模块,win7_32位下载后一步步安装即可,直到安装完成。简单的例子:
from scipy import *
import numpy as np
from pylab import *
a = zeros(1000)
a[:10]=1
b=fft(a)
plot(abs(b))
show()
(1)最小二乘法拟合数据
scipy中的子函数库optimize已经提供了实现最小二乘拟合算法的函数leastsq。函数原型为:
leastsq(func, x0, args=(), Dfun=None, full_output=0, col_deriv=0, ftol=1.49012e-08, xtol=1.49012e-08, gtol=0.0, maxfev=0, epsfcn=0.0, factor=100, diag=None, warning=True)
一般我们只要指定前三个参数就可以了。
func 是我们自己定义的一个计算误差的函数,
x0 是计算的初始参数值
args 是指定func的其他参数
下面用一个示例程序来解释他的用法,同样也是使用上面的求一次函数的拟合参数
# -*- coding: utf-8 -*-
import numpy as np
from scipy.optimize import leastsq
import pylab as pl
def func(x,p):
"""数据拟合所用的函数: A*sin(2*pi*k*x + theta)"""
A, k, theta = p
return A*np.sin(2*np.pi*k*x+theta)
def errors(p,y,x):
"""实验数据x, y和拟合函数之间的差,p为拟合需要找到的系数"""
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进行数据拟合,errors为计算误差的函数,p0为拟合参数的初始值
#args为需要拟合的实验数据
plsq =leastsq(errors, p0, args=(y1, x))
print("real pares:", [A, k, theta])
print("fitting paras", plsq[0]) # 实验数据拟合后的参数
pl.plot(x, y0, label="real data")
pl.plot(x, y1, label="noise data")
pl.plot(x, func(x,plsq[0]),label="fitting data")
pl.legend()
pl.show()
F5运行结果如下:
>>>
('real pares:', [10, 0.34, 0.5235987755982988])
('fitting paras', array([-10.51740809, 0.34191249, 3.66939538]))
注:由于正弦函数具有周期性,实际上拟合参数得到的函数和真实参数对应的函数基本是一致的。
(2)非线性方程组求解
optimize库中的fsolve函数可以用来对非线性方程组进行求解。它的基本调用形式如下:
fsolve(func, x0)
func(x)是计算方程组误差的函数,它的参数x是一个矢量,表示方程组的各个未知数的一组可能解,func返回将x代入方程组之后得到的误差;x0为未知数矢量的初始值。如果要对如下方程组进行求解的话:
f1(u1,u2,u3) = 0
f2(u1,u2,u3) = 0
f3(u1,u2,u3) = 0
那么func可以如下定义:
def func(x):
u1,u2,u3 = x
return [f1(u1,u2,u3), f2(u1,u2,u3), f3(u1,u2,u3)]
例如,求解如下方程组的解:
x1-2 = 0
x0*x0-*sin(x1*x2) = 0
x1*x2-2 = 0
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 [
x1-2,
x0*x0-sin(x1*x2),
x1*x2-2
]
result = fsolve(f, [1,1,1]) #设置初始值为[1 1 1]
print (result)
print (f(result))
运行结果如下:
>>>
[ 0.95357088 2. 1. ]
[0.0, -6.661338147750939e-16, -4.440892098500626e-16]
参考: