طرق الفروق المنتهية

المشتقات على شبكة

المشتقة نهاية، لكن الحاسوب يعمل بخطوات منتهية. تستبدل الفروق المنتهية f'(x) بالميل على خطوة صغيرة h: أمامية (f(x + h) - f(x)) / h أو مركزية (f(x + h) - f(x - h)) / (2h) وهي أدق لأن خطأي الجانبين يلغيان بعضهما (ينقص خطؤها مثل h**2 بينما الأمامية مثل h فقط). وتعطي الفكرة نفسها المشتقة الثانية، أي التحدب: (f(x + h) - 2 f(x) + f(x - h)) / h**2. وبالنسبة إلى سعر خيار كدالة في السعر الفوري فهذه هي غاما.

import math

def forward(f, x, h):
    return (f(x + h) - f(x)) / h

def central(f, x, h):
    return (f(x + h) - f(x - h)) / (2 * h)

exact = math.cos(1.0)                               # the derivative of sin at 1
for h in (0.1, 0.01, 0.001):
    print(h, abs(forward(math.sin, 1.0, h) - exact), abs(central(math.sin, 1.0, h) - exact))

# Gamma of a call from the second difference of its Black-Scholes price
import math

def norm_cdf(x):
    return 0.5 * (1 + math.erf(x / math.sqrt(2)))

def bs_call(spot, strike, rate, sigma, t):
    d1 = (math.log(spot / strike) + (rate + 0.5 * sigma ** 2) * t) / (sigma * math.sqrt(t))
    return spot * norm_cdf(d1) - strike * math.exp(-rate * t) * norm_cdf(d1 - sigma * math.sqrt(t))

def second_difference(f, x, h):
    return (f(x + h) - 2 * f(x) + f(x - h)) / h ** 2

price = lambda s: bs_call(s, 100, 0.05, 0.2, 1)
print(round(second_difference(price, 100, 0.5), 6))     # 0.018761
print(round(second_difference(price, 100, 0.01), 6))    # 0.018762

من مشتقة إلى معادلة

هدف الطريقة حل المعادلات التفاضلية الجزئية. وبلاك-شولز واحدة منها: يتغير سعر الخيار مع الزمن بحسب تحدبه بالنسبة إلى السعر الفوري. وبعد تغيير المتغيرات تصبح معادلة الحرارة du/dt = D * d2u/dx2: الحرارة تنتشر، وكلما كان منحنى الحرارة أشد تحدبًا تحرك أسرع. نضع القيم على شبكة ونستبدل المشتقتين بفروق، فتصير المعادلة قاعدة لتحديث كل نقطة انطلاقًا من جاراتها.

# The explicit scheme: a spike of heat spreads to its neighbours, step by step
import numpy as np

def explicit_step(values, alpha):
    new = values.copy()
    new[1:-1] = values[1:-1] + alpha * (values[2:] - 2 * values[1:-1] + values[:-2])
    return new

grid = np.zeros(9)
grid[4] = 1.0
for step in range(4):
    print(step, np.round(grid, 3))
    grid = explicit_step(grid, 0.25)

في الكود يربط alpha = D * dt / dx**2 بين الخطوة الزمنية dt وتباعد الشبكة dx. وتحديث واحد، new[1:-1] = values[1:-1] + alpha * (values[2:] - 2 * values[1:-1] + values[:-2])، يحرك كل النقاط الداخلية دفعة واحدة بالشرائح، والطرفان هما الشروط الحدية وقد ثُبّتا هنا. المخطط بسيط، والمأزق كذلك.

الاستقرار: المأزق

المخطط الصريح مستقر فقط إذا كان alpha <= 0.5، أي D * dt / dx**2 <= 0.5. وما بعد ذلك تتضخم أخطاء التقريب في كل خطوة وتنفجر الأرقام بنمط سنّ المنشار. والشرط قاسٍ: فتنصيف تباعد الشبكة يفرض خطوة زمنية أصغر بأربع مرات، فتتزايد الكلفة بسرعة. ولهذا يستخدم الكود الاحترافي غالبًا مخططات ضمنية تحل نظامًا خطيًا صغيرًا في كل خطوة وتستقر مع أي dt؛ أما الصريح فأبسط في الكتابة والفهم.

# The same spike with alpha = 0.4 (stable) and alpha = 0.6 (unstable) after 60 steps
import numpy as np

def explicit_step(values, alpha):
    new = values.copy()
    new[1:-1] = values[1:-1] + alpha * (values[2:] - 2 * values[1:-1] + values[:-2])
    return new

for alpha in (0.4, 0.6):
    grid = np.zeros(21)
    grid[10] = 1.0
    for _ in range(60):
        grid = explicit_step(grid, alpha)
    print(alpha, "largest value", float(np.abs(grid).max()))       # 0.4: 0.0556, 0.6: about 3.1e7

• تستبدل الفروق المنتهية المشتقة بقيم على شبكة: المركزية (f(x + h) - f(x - h)) / (2h) خطؤها نحو h**2؛ والفرق الثاني (f(x + h) - 2 f(x) + f(x - h)) / h**2 يعطي التحدب (غاما).
• تصبح المعادلات التفاضلية الجزئية قواعد تحديث على شبكة: الخطوة الصريحة لمعادلة الحرارة هي values[i] + alpha * (values[i + 1] - 2 * values[i] + values[i - 1]) مكتوبة بالشرائح.
• يتطلب المخطط الصريح D * dt / dx**2 <= 0.5 ليبقى مستقرًا؛ وترفع المخططات الضمنية هذا الحد مقابل حل نظام خطي في كل خطوة.

التالي: مصفوفتا التغاير والارتباط، الخطوة الأولى في تحسين المحافظ.

المختبر: ما الذي يُقيَّم

• قدّر التحدب باستخدام second_difference25 نقطة
• حرّك الشبكة خطوة واحدة باستخدام explicit_step25 نقطة
• تحقق من الخطوة الزمنية باستخدام is_stable25 نقطة
• تقدّم عبر الزمن باستخدام solve_heat25 نقطة