python求取Black-Scholes方程数值解，并画出不同波动率下期权价格的曲线

Black-Scholes方程是一个偏微分方程，求解它的数值解需要使用数值方法。这里我们使用有限差分法，将偏微分方程离散化为一个差分方程，然后用迭代法求解。

首先，我们需要导入一些必要的库：

import numpy as np
import matplotlib.pyplot as plt

然后，我们定义一些模型参数：

S0 = 100.0  # 股票初始价格
K = 100.0  # 行权价格
r = 0.05  # 无风险利率
sigma = 0.2  # 波动率
T = 1.0  # 到期时间

接下来，我们定义离散化参数：

N = 1000  # 离散化网格数
M = 1000  # 时间网格数
dt = T / M  # 时间步长
ds = S0 / N  # 价格步长

然后，我们创建一个数组来存储股票价格和时间的网格：

s = np.arange(0, S0 * 2, ds)
t = np.arange(0, T + dt, dt)

接下来，我们创建一个二维数组来存储期权价格的数值解。数组的第一维表示时间，第二维表示股票价格：

v = np.zeros((M + 1, N + 1))

我们初始化期权价格的边界条件。在到期时间，期权价格等于max(S-K, 0)，而在股票价格为0时，期权价格为0：

v[-1, :] = np.maximum(s - K, 0)
v[:, 0] = 0

然后，我们使用有限差分法来迭代求解期权价格。有限差分法的基本思想是将偏微分方程中的导数用差商近似。具体来说，我们使用向前差分法来近似时间导数，使用中心差分法来近似空间导数。

for i in range(M - 1, -1, -1):
    for j in range(1, N):
        alpha = 0.5 * sigma**2 * j**2
        beta = r * j
        gamma = r + alpha - beta
        v[i, j] = v[i+1, j] + dt * (alpha * v[i+1, j+1] + gamma * v[i+1, j] + beta * v[i+1, j-1]) / j**2

最后，我们画出不同波动率下期权价格的曲线：

plt.plot(s, v[0, :])
plt.xlabel('Stock Price')
plt.ylabel('Option Price')
plt.title('Option Price vs. Stock Price')
plt.show()

完整代码如下：