python数值计算_Python数值计算——Scipy库的应用

Scipy在 Numpy的基础上增加了众多的数学计算、利学计算以及工程计算中常用的模块, 例如线性代数、常微分方程数值求解、信号处理、图像处理、稀疏矩阵等。今天我们用几个数值计算常用的算法来体验一下python的Scipy库强大之处。拟合与优化一optimize

非线性方程组求解

且看如下图所示的一组非线性方程，手算求解的估计要费九牛二虎之力，我们看一下如何利用Scipy库里面的优化函数进行求解，这段代码大家可以修改之后用于自己需要求解非线性方程组：

from math import sin, cos

from scipy import optimize

####定义非线性方程组

def f(x):

#######tolist的作用是将list转化为浮点型的list

x0, xl, x2 = x.tolist()

return [

5*xl+3,

4*x0*x0 - 2*sin(xl*x2),

xl*x2 - 1.5

]

# f 计算方程组的误差，[ 1 ,1 ,1 ]是未知数的初始值

result = optimize.fsolve(f, [1,1,1])

######输出x0 x1 x2

print (result)

#####检验求解的结果的正确性

print (f(result))

最小二乘拟合

最小二乘法在数据分析中非常常见，假如我们知道一组数据他们符合某种函数关系，我们可以通过最小二乘法去确定该函数中的某些参数，例如下面的一组数据，我们假设其符合Y=ax+b这种线性函数关系，利用最小二乘法我们可以确定a与b的值，我们可以调用了Scipy库中的 optimize.leastsq函数，除了这个函数以外还有专门的curve_fit函数，curve_fit函数更加直白一点，curve_fit你只需要定义你需要拟合的函数，而optimize.leastsq你需要定义的是残差函数。下面的两段代码分别采用optimize.leastsq函数和curve_fit函数拟合。

########利用optimize.leastsq函数进行拟合

import numpy as np

from scipy import optimize

import matplotlib.pyplot as plt

X = np.array([ 8.19, 2.72, 6.39, 8.71, 4.7 , 2.66, 3.78])

Y = np.array([ 7.01, 2.78, 6.47, 6.71, 4.1 , 4.23, 4.05])

"计算以p为参数的直线和原始数据之间的残差"

def residuals(p):

k,b = p

return Y - (k*X + b)

# l e a s t s q使得residuals()的输出数姐的平方和最小，参数的初始值为[1,0]

r = optimize.leastsq(residuals, [1, 0])

k, b = r[0]

print ("k =" ,k, "b =" ,b)

#####将计算结果可视化

XX=np.arange(1,10,0.1)

YY=k*XX+b

plt.plot(XX,YY)

plt.plot(X,Y,'o')

plt.xlabel('X')

plt.ylabel('Y')

########利用curve_fit函数进行拟合

from scipy import optimize

import matplotlib.pyplot as plt

import numpy as np

X = np.array([ 8.19, 2.72, 6.39, 8.71, 4.7 , 2.66, 3.78])

Y = np.array([ 7.01, 2.78, 6.47, 6.71, 4.1 , 4.23, 4.05])

"定义拟合函数"

def func(x,k,b):

return k*x+b

# #popt数组中，两个值分别是待求参数k,b

popt, pcov = optimize.curve_fit(func,X , Y)

print ("k =" ,popt[0], "b =" ,popt[1])

XX=np.arange(1,10,0.1)

YY=popt[0]*XX+popt[1]

plt.plot(XX,YY)

plt.plot(X,Y,'o')

plt.xlabel('X')

plt.ylabel('Y')

线性函数的拟合其实比较简单，我们数据分析中数据可能更加符合其他类型的函数，比如指数函数，正弦函数等，能否用最小二乘法进行拟合呢，答案当然是可以的。下面展示一个指数函数的拟合。希望大家仔细品味这寥寥几行代码。这里同样给出两种方法。optimize.leastsq函数进行拟合的时候初值的不同会影响拟合结果，大家可以自行尝试。

#########利用optimize.leastsq函数进行拟合

from scipy import optimize

import matplotlib.pyplot as plt

import numpy as np

X=np.linspace(0,4,50)

Y=2.5*np.exp(-1.3*X)+0.5+0.2*np.random.normal(size=len(X))

def residuals(p):

a,b,c= p

return Y-(a*np.exp(-b*X)+c)

# l e a s t s q使得residuals()的输出数姐的平方和最小，参数的初始值为[1,0,1]

r = optimize.leastsq(residuals, [1,0,1])

a,b,c = r[0]

print ("a =" ,a, "b =" ,b,'c=',c)

XX=np.arange(0,4,0.1)

YY=a*np.exp(-b*XX)+c

plt.plot(XX,YY)

plt.plot(X,Y,'o')

plt.xlabel('X')

plt.ylabel('Y')

#########利用curve_fit函数进行拟合

from scipy import optimize

import matplotlib.pyplot as plt

import numpy as np

X=np.linspace(0,4,50)

Y=2.5*np.exp(-1.3*X)+0.5+0.2*np.random.normal(size=len(X))

def func(x,a,b,c):

return (a*np.exp(-b*x)+c)

# l e a s t s q使得residuals()的输出数姐的平方和最小，这里参数的初始值为[1,0,1]

popt, pcov =optimize.curve_fit(func, X,Y)

a,b,c = popt

print ("a =" ,a, "b =" ,b,'c=',c)

XX=np.arange(0,4,0.1)

YY=a*np.exp(-b*XX)+c

plt.plot(XX,YY)

plt.plot(X,Y,'o')

plt.xlabel('X')

plt.ylabel('Y')

解线性方程组

学过线性代数的同学可能体会过求解线性方程组的痛苦，矩阵的求逆害苦了不少同学，这里教大家Scipy的求解线性方程组的函数。其实方程组的求解一般分为两种，矩阵求逆法和直接求解，这里也给出了两种方法供大家参考。

import numpy as np

from scipy import linalg

m,n = 500, 50

A = np.random.rand(m, m)

B = np.random.rand(m, n)

#######直接采用linalg.solve函数

X1 = linalg.solve(A, B)

#######采用矩阵求逆的方法

X2 = np.dot(linalg.inv(A), B)

######np.allclose这个函数是用来对比两个矩阵是不是每个元素都相同，如果都相同输出“True”

print (np.allclose(X1, X2))

Scipy-linalg里面还有很多处理矩阵的函数，比如：linalg.eig求矩阵的特征值与特征向量等等强大的函数，大家可以自己去慢慢挖掘。

求解常微分方程组

如果你学过数值分析的话你应该接触到很多求解微分方程的方法，欧拉法，龙格库塔等等，这里Scipy里面有现成的函数供我们使用（odeint）,以下面一组一阶常微分方程组为例展示http://ode.int函数的用法。如果是二阶的常微分方程，需要将二阶的常微分方程转化为一阶的常微分方程，一个二阶的常微分方程可以转化为两个一阶的常微分方程（状态方程）。这里展示了两种求解的方法，第一种是一步一步的求解，对于做PID控制的同学比较常用，第二种是一次性求解所有时间步，用到的函数是一样的，如果不是做控制的小伙伴可以直接忽略第一种方法。

其中，S()为越阶函数

# -*- coding: utf-8 -*-

"""

Created on Wed Mar 4 22:07:58 2020

@author: Administrator

"""

import numpy as np

from scipy.integrate import odeint

import matplotlib.pyplot as plt

# 定义微分方程组

def model(z,t,u):

x = z[0]

y = z[1]

dxdt = (-x + u)/2.0

dydt = (-y + x)/5.0

dzdt = [dxdt,dydt]

return dzdt

# 初始条件

z0 = [0,0]

# 时间的步数

n = 401

# 时间点

t = np.linspace(0,40,n)

# 定义函数u,跃阶函数

u = np.zeros(n)

u[51:] = 2.0

# 存储求解结果的空数组

x = np.empty_like(t)

y = np.empty_like(t)

# 定义初始状态

x[0] = z0[0]

y[0] = z0[1]

# 求解

for i in range(1,n):

# 每次求解一个步长

tspan = [t[i-1],t[i]]

# 求解下一个步长

z = odeint(model,z0,tspan,args=(u[i],))

# store solution for plotting

x[i] = z[1][0]

y[i] = z[1][1]

# 下一个步长的初始条件

z0 = z[1]

# 结果可视化

plt.plot(t,u,'g:',label='u(t)')

plt.plot(t,x,'b-',label='x(t)')

plt.plot(t,y,'r--',label='y(t)')

plt.ylabel('values')

plt.xlabel('time')

plt.legend(loc='best')

plt.show()

上面是一步步求解的代码，第二种是一次性求解所有时间步，代码较第一种更为简洁。大家可以仔细品味一下哦，最好是能动手敲一遍，我可不止一次告诉你们模仿的重要性了哟。

import numpy as np

from scipy.integrate import odeint

import matplotlib.pyplot as plt

# 定义微分方程组与定义函数u,跃阶函数

def model(z,t):

if t>5:

u= 2.0

else:

u= 0.0

x = z[0]

y = z[1]

dxdt = (-x + u)/2.0

dydt = (-y + x)/5.0

dzdt = [dxdt,dydt]

return dzdt

# 初始条件

z0 = [0,0]

# 时间的步数

n = 401

# 时间点

t = np.linspace(0,40,n)

# 存储求解结果的空数组

x = np.empty_like(t)

y = np.empty_like(t)

# 求解

z = odeint(model,z0,t)

x= z[:,0]

y= z[:,1]

# 结果可视化

plt.plot(t,x,'b-',label='x(t)')

plt.plot(t,y,'r--',label='y(t)')

plt.ylabel('values')

plt.xlabel('time')

plt.legend(loc='best')

plt.show()

当然了Scipy的功能还有很多，比如数据的统计分析，信号的频域变换，图形的处理等等都可以用Scipy来实现，这里不一一介绍了，大家有兴趣的可以去Scipy的官网去看看，链接我放在下面了，里面可是一个大宝库，一般人我不告诉他。说了这么多这里还是要强调一下，python的强大由于它的丰富的库函数，所以一定要学会如何调用这些库函数哟。https://www.scipy.org/​www.scipy.org

看到这里了不如点个赞吧，希望这篇文章对你有用。