Python Scipy 练习

注：本文为国科大《Python 科学计算与数据处理》课程 Scipy 作业的代码，仅供学习参考，如有错误欢迎指出。

题目1

使用 scipy.optimize.fsolve 求解非线性方程组：

$$\{\begin{array}{l}\cos (a)=1-d^2/(2\cdot r^2)\\ L=a\cdot r\\ d=140\\ L=156\end{array}$$

from scipy.optimize import fsolve
from math import cos, sin

def f(x):
    d = 140
    l = 156
    a, r = x.tolist()
    return [cos(a)-1+d**2/(2*r**2), l-a*r]

# 不导入 jacobi 矩阵
ans1 = fsolve(f, [1, 1])
print(ans1)
print(f(ans1))

# 导入 jacobi 矩阵

def j(x):
    d = 140
    a, r = x.tolist()
    return [[-sin(a), -(d**2)/(r**3)], [-r, -a]]  # 雅可比矩阵

ans2 = fsolve(f, [1,1], fprime=j)
print(ans2)
print(f(ans2))

题目2

用 curve_fit 对高斯分布进行拟合，高斯分布函数为 $y=a\cdot \exp (-(x-b{)}^2/(2c^2)),x\in [0,10]$， 其中真实值 $a=1,b=5,c=2$。试对 $y$ 加入噪声之后进行拟合，并作图与真实数据进行比较。

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

def gauss(x, a, b, c):
    return a * np.exp(-(x-b)**2/(2*c**2))

x = np.linspace(0, 10, 1000)
y = gauss(x, 1, 5, 2) + np.random.normal(size=len(x)) * 0.1

popt, pcov = curve_fit(gauss, x, y)
y2 = gauss(x, popt[0], popt[1], popt[2])

plt.plot(x, y, 'bo')
plt.plot(x, y2, 'r')
plt.show()

题目3

对4个数据点x = [-1, 0, 2.0, 1.0]，y = [1.0, 0.3, -0.5, 0.8]进行Rbf插值，插值中使用三种插值方法分别是 multiquadric、gaussian、和 linear，在np.linspace(-3, 4, 100)上作图。

from scipy.interpolate import Rbf
import numpy as np
import matplotlib.pyplot as plt

x = [-1, 0, 2.0, 1.0]
y = [1.0, 0.3, -0.5, 0.8]
func1 = Rbf(x, y, function='multiquadric')
func2 = Rbf(x, y, function='gaussian')
func3 = Rbf(x, y, function='linear')

i = np.linspace(-3, 4, 100)

plt.plot(i, func1(i), 'r', label="multiquadric")
plt.plot(i, func2(i), 'g', label="gaussian")
plt.plot(i, func3(i), 'b', label="linear")
plt.plot(x, y, 'ko')
plt.legend()
plt.show()

题目4

分别用 optimize.fmin_bfgs、optimize.fminbound、optimize.brute 三种优化方法对函数 $x^2+10\sin (x)，x\in [-10,10]$ 求最小值，并作图验证。

from scipy.optimize import fmin_bfgs, fminbound, brute
import numpy as np
import matplotlib.pyplot as plt

def fun(x):
    return x ** 2 + 10 * np.sin(x)

x1 = fmin_bfgs(fun, np.array([0]))[0]  # 参数：初始猜测值
x2 = fminbound(fun, -10, 10)  # 参数：优化边界

ranges = (slice(-10, 10.1, 0.1),)
x3 = brute(fun, ranges)[0]  # ranges：元组类型，每个元素必须为 slice 或者 range

print('x1:', x1)
print('x2:', x1)
print('x3:', x1)

x_list = [x1, x2, x3]
text_list = ['fmin_bfgs', 'fminbound', 'brute']
color_list = ['ro', 'go', 'bo']

for i in range(3):
    x0 = np.linspace(-10, 10, 2000)
    plt.plot(x0, fun(x0), 'y')

    x, y = x_list[i], fun(x_list[i])
    plt.plot(x, y, color_list[i])
    plt.annotate(text_list[i], xy=(x, y), xytext=(x+1.5, y), arrowprops=dict(arrowstyle='->', facecolor='black'))
    plt.show()

参考：
Scipy Manual - scipy.optimize.fmin_bfgs
Scipy Manual - scipy.optimize.fminbound
Scipy Manual - scipy.optimize.brute

题目5

计算积分：a) ${\int }_0^3{\cos }^2(e^x)dx$，b) ${\int }_0^{0.5}dy{\int }_0^{\sqrt{1-4y^2}}16xydx$

# 计算第一个积分

from scipy import integrate
import numpy as np

ans = integrate.quad(lambda x: np.cos(np.exp(x)) ** 2, 0, 3)
print(ans[0])

# 计算第二个积分

from scipy import integrate
import numpy as np

# 方法一：integrate.dblquad 计算二重积分
f = lambda x, y: 16 * x * y
g = lambda x: 0
h = lambda y: np.sqrt(1 - 4 * y ** 2)
ans = integrate.dblquad(f, 0, 0.5, g, h)
print(ans[0])

# 方法二：integrate.nquad 计算多重积分
f = lambda x, y: 16 * x * y
bounds_x = lambda y: [0, np.sqrt(1 - 4 * y ** 2)]
bounds_y = [0, 0.5]
ans = integrate.nquad(f, [bounds_x, bounds_y])
print(ans[0])

参考：
Scipy Manual - scipy.integrate.quad
Scipy Manual - scipy.integrate.dblquad
Scipy Manual - scipy.integrate.nquad

题目6

弹簧系统每隔1ms周期的系统状态 ，试用 scipy.integrate.odeint 对该系统进行求解并作图，
其中参数(M, k, b, F) = (1.0, 0.5, 0.2, 1.0)，初始状态init_status = [-1, 0.0]，t = np.arange(0, 50, 0.02)。

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def fun(init_state, t, M, k, b, F):
    theta, omega = init_state
    dydt = [omega, (F - b * omega - k * theta) / M]
    return dydt

y0 = [-1, 0]  # 初值(x, x')
t = np.arange(0, 50, 0.02)
args = (1, 0.5, 0.2, 1)  # 参数(M, k, b, F)

x = odeint(fun, y0, t, args)

plt.plot(t, x[:, 0], 'r', label="x")
plt.plot(t, x[:, 1], 'g', label="x\'")
plt.legend()
plt.xlabel('t')
plt.grid()
plt.show()

参考：Scipy Manual - scipy.integrate.odeint

题目7

以参数为1的伽马分布生成1000个随机数，绘制这些样点的直方图，在其上绘制此伽马分布的 PDF（概率密度函数）

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma

fig, ax = plt.subplots(1, 1)

a = 1.0
x = np.linspace(gamma.ppf(0.01, a), gamma.ppf(0.99, a), 100)
ax.plot(x, gamma.pdf(x, a), 'r', label='gamma pdf')

r = gamma.rvs(a, size=1000)
ax.hist(r, bins=50, density=True, histtype='stepfilled', alpha=0.2)
ax.legend(loc='best', frameon=False)
plt.show()

参考：Scipy Manual - scipy.stats.gamma

题目8

scipy.sparse 提供了多种表示稀疏矩阵的格式，试用 dok_matrix、lil_matrix 表示矩阵[[3 0 8 0] [0 2 0 0] [0 0 0 0] [0 0 0 1]]，并与 coo_matrix 表示法进行比较。

import numpy as np
from scipy.sparse import coo_matrix, dok_matrix, lil_matrix

a = np.matrix([[3, 0, 8, 0], [0, 2, 0, 0], [0, 0, 0, 0], [0, 0, 0, 1]])

doc = dok_matrix(a)
print('doc:')
print(doc.items())

lil = lil_matrix(a)
print('\nlil:')
print('rows:', lil.rows)
print('data:', lil.data)

coo = coo_matrix(a)
print('\ncoo:')
print('col:', coo.col)
print('row:', coo.row)
print('data:', coo.data)

doc_matrix： 采用字典存储矩阵中的非0元素，字典的 key 是元素位置信息的元组，value 是元素值。

lil_matrix： 使用两个列表存储非0元素，列表 rows 保存非零元素所在的列，列表 data 存储元素值。

coo_matrix： 采用三元组 (row, col, data) 的形式存储非0元素，数组 row、col 和 data 分别保存非0元素的行下标、列下标与值。

参考：【SciPy】Sparse稀疏矩阵主要存储格式总结