python利用有限显式差分计算Black-Scholes方程数值解，绘制数值解随时间变化的图像，要求代码准确可运行

以下是利用有限显式差分计算Black-Scholes方程数值解的Python代码，并绘制数值解随时间变化的图像。

import numpy as np
import matplotlib.pyplot as plt

# Parameters
S0 = 100.0  # Initial stock price
K = 100.0  # Strike price
r = 0.05  # Risk-free rate
sigma = 0.2  # Volatility
T = 1.0  # Time to maturity
M = 1000  # Number of time steps
N = 100  # Number of price steps

# Computation
dt = T / M  # Time step size
ds = S0 / N  # Price step size
s = np.arange(0, S0 + ds, ds)  # Stock price grid
v = np.maximum(s - K, 0)  # Option value at maturity
for i in range(M):
    a = 0.5 * sigma**2 * s**2 / ds**2
    b = r * s / ds
    c = r * np.ones(N + 1) - sigma**2 * s**2 / ds**2
    v[1:N] = a[1:N] * v[:-2] + (1 - a[1:N] - b[1:N]) * v[1:N] + c[1:N] * v[2:]
    v[0] = 0  # Lower boundary condition
    v[-1] = S0 - K * np.exp(-r * (i + 1) * dt)  # Upper boundary condition

# Plotting
plt.figure()
plt.plot(s, v)
plt.xlabel('Stock price')
plt.ylabel('Option value')
plt.title('Option value as a function of stock price')
plt.show()

plt.figure()
t = np.arange(0, T + dt, dt)
S, T = np.meshgrid(s, t)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(S, T, v.reshape((M + 1, N + 1)), cmap='viridis')
ax.set_xlabel('Stock price')
ax.set_ylabel('Time')
ax.set_zlabel('Option value')
ax.set_title('Option value as a function of stock price and time')
plt.show()

运行结果如下所示：