Читать книгу Numerical Methods in Computational Finance - Daniel J. Duffy - Страница 61

We show hand-crafted code for the explicit Euler method as well as the schemes (3.22), (3.23) and (3.24). The main reason is to get hands-on experience with coding solvers for ODEs before using production solvers (for example, Python's or Boost C++ odeint libraries), especially for readers for whom numerical ODEs are new, for example readers with an economics or econometrics background. A good strategy is 1) get it working (write your own code and examples), then 2) get it right (use a library with the same examples) and then and only then get it optimised (use a library in a larger application).

The following code is for the methods:

 Explicit Euler

 Fourth-order Runge–Kutta (RK4)

 Second-order Ralston.

 Heun (improved Euler)

 Predictor-corrector method.

We now present the code for a scalar IVP for these methods:

# y' = f(x,y) Explicit Euler method, scalar def Euler(f, x, y, T, N, TOL): h = (T - x) / N for n in range(1,N+1): ynP1 = y + h*f(x,y) # Go to next level x += h y = ynP1 return ynP1 # y' = f(x,y) Explicit Euler method, system def Euler2(f, x, y0, T, N, TOL): h = (T - x) / N m = len(y0) ynP1 = y0 for n in range(N): for j in range(m): ynP1[j] = y0[j] + h*f[j](x,y0) # Go to next level x += h y0 = ynP1 return ynP1 def RK4(f, x, y, T, N, TOL): # Order 4 Runge Kutta method h = (T - x) / N k1 = 0; k2 = 0; k3 = 0; k4 = 0 ynP1 = 0 # y(n+1) for n in range(N): k1 = f(x, y); k2 = f(x + h / 2, y + h*k1 / 2) k3 = f(x + h / 2, y + h*k2 / 2); k4 = f(x + h, y + h*k3); ynP1 = y + h*(k1 + 2*k2 + 2*k3 + k4)/6 # Go to next level x += h; y = ynP1 return ynP1 # https://en.wikipedia.org/wiki/Heun%27s_method def Heun(f, x, y, T, N, TOL): # Modified Euler method (improved polygon method), Collatz 1960 h = (T - x) / N # Variables k1 = 0; k2 = 0 for n in range(N): ynP1 = y + h*f(x+h/2, y + h*f(x,y)/2) # Go to next level x += h; y = ynP1 return ynP1 def Ralston2(f, x, y, T, N, TOL): h = (T - x) / N # Variables k1 = 0; k2=0 ynP1=0 # y(n+1) for n in range(N): k1 = h*f(x, y); k2 = h*f(x + 2*h / 3, y + 2*k1 / 3) ynP1 = y + (k1 + 3 * k2 ) / 4 # Go to next level x += h; y = ynP1 return ynP1; def PredictorCorrector(f, x, y, T, N, TOL): # PC with fixed step size h = (T - x) / N for n in range(1,N+1): yCorr2 = y + h*f(x, y) # 1. Predictor (explicit Euler) # Produce y(n+1) at level x(n+1) from level n, Inner iteration diff = 2 * TOL while (diff > TOL): yCorr = y + 0.5*h*(f(x, y) + f(x + h, yCorr2)) diff = math.fabs(yCorr - yCorr2) / math.fabs(yCorr) yCorr2 = yCorr # Stuff has converged at level n+1, now update next level x += h y = yCorr return yCorr

A sample test program is:

# TestFDMSchemes.py # # A catalogue of hand-crafted numerical methods for solving ODEs # # (C) Datasim Education BV 2021 # # Testing import math import numpy as np import FDMSchemes as fdm def func(t,y): return y*(1.0 - y) # logistic ode t = 0; y = 0.5 #initial time and value T = 4.0 #end of integration interval N = 4000 #number of divisions of [0,T] TOL = 1.0e-5 #for some method, a measure of the desired accuracy # y' = f(x,y) Explicit Euler method value = fdm.Euler(func,t,y,T,N,TOL) #scalar print(value) value = fdm.PredictorCorrector(func,t,y,T,N,TOL) #scalar print(value) value = fdm.RK4(func,t,y,T,N,TOL) #scalar print(value) y = 0.5 value = fdm.Heun(func,t,y,T,N,TOL) #scalar print(value) value = fdm.Ralston2(func,t,y,T,N,TOL) #scalar print(value)

No doubt the code can be modified to make it more “pythonic” (whatever that means) and reduce the number of lines of code (at the possible expense of readability) if that is a requirement. We have also written these algorithms in C++11 as discussed in Duffy (2018).