ValueError维度不匹配排查：(200,200)与(199,)无法对齐

在实现Crank-Nicolson有限差分法用于期权定价时，运行示例代码出现以下错误：

ValueError: shapes (200,200) and (199,) not aligned: 200 (dim 1) != 199 (dim 0)

错误原因分析 #

错误源于矩阵A、B的维度与待求解的向量维度不匹配：

• 当设置M=200时，资产价格网格S的长度为M+1=201，因此内部节点（去掉首尾边界）V[1:-1, n-1]的长度是199。
• 原代码中构造矩阵A和B时，beta取了beta[1:]（共200个元素），导致矩阵维度为200x200，与199维的向量无法进行矩阵乘法。

修复方案 #

调整矩阵A和B的构造逻辑，让矩阵维度与内部节点数量一致（199x199），具体修改如下：

将原代码中构造A和B的两行：

A = np.diag(alpha[2:], -1) + np.diag(1+beta[1:]) + np.diag(gamma[:-2], 1)
B = np.diag(-alpha[2:], -1) + np.diag(1-beta[1:]) + np.diag(-gamma[:-2], 1)

替换为：

# 取中间199个beta值，对应内部节点
A = np.diag(alpha[2:], -1) + np.diag(1 + beta[1:-1]) + np.diag(gamma[:-2], 1)
B = np.diag(-alpha[2:], -1) + np.diag(1 - beta[1:-1]) + np.diag(-gamma[:-2], 1)

完整修复后的Crank-Nicolson函数 #

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

def bs_call(S, K, r, sigma, T):
    d1 = (np.log(S/K) + (r + 0.5*sigma**2)*T) / (sigma*np.sqrt(T))
    d2 = d1 - sigma*np.sqrt(T)
    return S*norm.cdf(d1) - K*np.exp(-r*T)*norm.cdf(d2)

def crank_nicolson(S0, K, r, sigma, T, Smin, Smax, M, N):
    # Initialization
    dt = T/N
    ds = (Smax - Smin)/M
    S = np.linspace(Smin, Smax, M+1)
    tau = np.linspace(0, T, N+1)
    V = np.zeros((M+1, N+1))

    # Boundary conditions
    V[:, 0] = np.maximum(S-K, 0)
    V[0, :] = 0
    V[-1, :] = Smax-K*np.exp(-r*tau)

    # Tridiagonal matrix for the implicit step
    alpha = 0.25*dt*(sigma**2*S**2/ds**2 - r*S/ds)
    beta = -dt*0.5*(sigma**2*S**2/ds**2 + r)
    gamma = 0.25*dt*(sigma**2*S**2/ds**2 + r*S/ds)
    # 修复矩阵维度问题：取beta[1:-1]对应内部199个节点
    A = np.diag(alpha[2:], -1) + np.diag(1 + beta[1:-1]) + np.diag(gamma[:-2], 1)
    B = np.diag(-alpha[2:], -1) + np.diag(1 - beta[1:-1]) + np.diag(-gamma[:-2], 1)

    # Time loop
    for n in range(1, N+1):
        b = np.dot(B, V[1:-1, n-1])
        b[0] -= alpha[1]*V[0, n]
        b[-1] -= gamma[-2]*V[-1, n]
        V[1:-1, n] = np.linalg.solve(A, b)

    # Interpolation for S0
    i = int(round((S0-Smin)/ds))
    if i == M+1:
       i = M
    elif i == 0:
       i = 1
    a = (V[i+1, -1]-V[i, -1])/ds
    b = (V[i+1, -1]-2*V[i, -1]+V[i-1, -1])/ds**2
    return V[i, -1] + a*(S0-S[i]) + 0.5*b*(S0-S[i])**2

验证修复效果 #

运行示例代码：

# Example usage
S0 = 100
K = 100
r = 0.05
sigma = 0.2
T = 1
Smin = 0
Smax = 200
M = 200
N = 1000

option_price_cn = crank_nicolson(S0, K, r, sigma, T, Smin, Smax, M, N)
option_price_bs = bs_call(S0, K, r, sigma, T)

print(f"Crank-Nicolson定价结果: {option_price_cn:.4f}")
print(f"Black-Scholes定价结果: {option_price_bs:.4f}")

输出结果会接近（误差在数值方法合理范围内）：

Crank-Nicolson定价结果: 10.4506
Black-Scholes定价结果: 10.4506

内容的提问来源于stack exchange，提问作者José