python解zuobiaoxi方程_Python数值计算——Scipy库的应用

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

1. 拟合与优化一optimize

非线性方程组求解

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

from 

最小二乘拟合

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

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

线性函数的拟合其实比较简单，我们数据分析中数据可能更加符合其他类型的函数，比如指数函数，正弦函数等，能否用最小二乘法进行拟合呢，答案当然是可以的。下面展示一个指数函数的拟合。希望大家仔细品味这寥寥几行代码。这里同样给出两种方法。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

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