title: SciPy-数值计算库 create: 2016.12.7 modified: 2016.12.7 tags: python 最小二乘拟合 数值积分

  8

[TOC]

SciPy函数库在NumPy库的基础上增加了众多的数学、科学以及工程计算中常用的库函数。例如线性代数、常微分方程数值求解、信号处理、图像处理、稀疏矩阵等等。由于其涉及的领域众多、本书没有能力对其一一的进行介绍。

1 最小二乘拟合

假设有一组实验数据(x[i], y[i]),我们知道它们之间的函数关系:y = f(x),通过这些已知信息,需要确定函数中的一些参数项。例如,如果f是一个线型函数f(x) = k*x+b,那么参数k和b就是我们需要确定的值。 scipy中的子函数库optimize已经提供了实现最小二乘拟合(Least-square fitting)算法的函数leastsq。下面是用leastsq进行数据拟合的一个例子:


In [3]:
# -*- 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 residuals(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进行数据拟合
# 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()


真实参数: [10, 0.34, 0.5235987755982988]
拟合参数 [ 10.09912902   0.3412268   -5.70118097]

这个例子中我们要拟合的函数是一个正弦波函数,它有三个参数 A, k, theta ,分别对应振幅、频率、相角。假设我们的实验数据是一组包含噪声的数据 x, y1,其中y1是在真实数据y0的基础上加入噪声的到了。
通过leastsq函数对带噪声的实验数据x, y1进行数据拟合,可以找到x和真实数据y0之间的正弦关系的三个参数: A, k, theta。上面是程序的输出:
我们看到拟合参数虽然和真实参数完全不同,但是由于正弦函数具有周期性,实际上拟合参数得到的函数和真实参数对应的函数是一致的。

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可以如下定义:


In [ ]:
def func(x):     
    u1,u2,u3 = x   
    return [f1(u1,u2,u3), f2(u1,u2,u3), f3(u1,u2,u3)]

下面是一个实际的例子,求解如下方程组的解:

5*x1 + 3 = 0
4*x0*x0 - 2*sin(x1*x2) = 0
x1*x2 - 1.5 = 0
程序如下:


In [10]:
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)


[-0.70622057 -0.6        -2.5       ]
[0.0, -9.126033262418787e-14, 5.329070518200751e-15]

由于fsolve函数在调用函数f时,传递的参数为数组,因此如果直接使用数组中的元素计算的话,计算速度将会有所降低,因此这里先用float函数将数组中的元素转换为Python中的标准浮点数,然后调用标准math库中的函数进行运算。

3 数值积分

数值积分是对定积分的数值求解,例如可以利用数值积分计算某个形状的面积。下面让我们来考虑一下如何计算半径为1的半圆的面积,根据圆的面积公式,其面积应该等于PI/2。单位半圆曲线可以用下面的函数表示:


In [13]:
def half_circle(x):
    return (1-x**2)**0.5

调用scipy.integrate库中的quad函数的话,将会得到非常精确的结果:


In [15]:
from scipy import integrate
pi_half, err = integrate.quad(half_circle, -1, 1)
pi_half*2


Out[15]:
3.141592653589797

4 解常微分方程组

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

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


In [17]:
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()

In [ ]: