In addition to my previous post ODESolver Ι try to implement antoher different numerical schemes. For example, I implement Runge-Kutta 4 order scheme for solver ODE class, and Forward Euler scheme. Basic superclass ODE_Solver code posted below:
#ODS.py
import numpy as np
class ODE_Solver(object):
"""
Supercalss sover ODE-s
Attribute:
x: array of x coords
u: solution array y(x)
k: number of steps
f: right hand side ODE equation: du/dx = f(u, x)
"""
def __init__(self, f):
if not callable(f):
# check correct function f(u, x)
raise TypeError('f is %s, not function' % type(f))
self.f = lambda u, x: np.asarray(f(u, x), float)
def solver_st(self):
"""Implement numerical scheme"""
raise NotImplementedError
def set_initial_condition(self, u0):
if isinstance(u0, (float, int)): # ODE 1-th order
self.neq = 1
u0 = float(u0)
else: # ODE high order or system of ODE-s
u0 = np.asarray(u0) # (initial conds)
self.neq = u0.size
self.u0 = u0
# check correct lenght of vector f
try:
f0 = self.f(self.u0, 0)
except IndexError:
raise IndexError(
'index out of bounds f(u,x). right index %s' % (str(range(self.neq))))
if f0.size != self.neq:
raise ValueError('f(u,x) return %d elems, u has %d elems' % (f0.size, self.neq))
def solve(self, coord_points, terminate=None):
"""
Solve equations. Default False
"""
if terminate is None:
terminate = lambda u, x, step_no: False
if isinstance(coord_points, (float, int)):
raise TypeError('solve: array lists not iterable')
self.x = np.asarray(coord_points)
if self.x.size <= 1:
raise ValueError('ODESolver.solve requre coords x array')
n = self.x.size
if self.neq == 1: # ODE
self.u = np.zeros(n)
else:
self.u = np.zeros((n, self.neq))
# Assume self.x[0] corresponds to self.u0
self.u[0] = self.u0
# looping for x coords
for k in range(n - 1):
self.k = k
self.u[k + 1] = self.solver_st()
if terminate(self.u, self.x, self.k + 1):
break
return self.u[:k + 2], self.x[:k + 2]
Forwar Euler scheme implement posted below:
# ES.py
from abc import ABC
from ODS import ODE_Solver
class FE(ODE_Solver, ABC):
"""
Attribute:
x: array of x coords
u: solution array y(x)
k: number of steps
f: right hand side ODE equation: du/dx = f(u, x)
"""
def solver_st(self):
u, f, k, x = self.u, self.f, self.k, self.x
dx = x[k+1] - x[k]
u_new = u[k] + dx*f(u[k], x[k])
return u_new
If I call Forvard Euler scheme it is workin, but now I try to implement Adams-Bashfort-Moulton Scheme:
#ADS.py
from abc import ABC
from ODS import ODE_Solver
class ABM4(ODE_Solver, ABC):
"""
Attribute:
x: array of x coords
u: solution array y(x)
k: number of steps
f: right hand side ODE equation: du/dx = f(u, x)
"""
def solver_st(self, alp=None):
u, f, k, x = self.u, self.f, self.k, self.x
for k in range(0,4):
dx = x[k + 1] - x[k]
dx2 = dx / 2.0
K1 = dx * f(u[k], x[k])
K2 = dx * f(u[k] + 0.5 * K1, x[k] + dx2)
K3 = dx * f(u[k] + 0.5 * K2, x[k] + dx2)
K4 = dx * f(u[k] + K3, x[k] + dx)
u = u[k] + (1 / 6.0) * (K1 + 2 * K2 + 2 * K3 + K4)
for k in range(4, self.x.size):
dx = x[k + 1] - x[k]
alp[k] = dx / 24 * (55 * f(u[k - 1], x[k - 1]) - 59 * f(u[k - 2], x[k - 2]) + 37 * f(u[k - 3], x[k - 3])
- 9 * f(u[k - 4], x[k - 4])) + u[k - 1]
u_new = dx / 24 * (9 * f(u[k], alp[k]) + 19 * f(u[k - 1], x[k - 1]) - 5 * f(u[k - 2], x[k - 2])
+ f(u[k - 3], x[k - 3])) + u[k - 1]
return u_new
And if I try to call this solver (Adams-Bashfort-Moulton) method from another class in main file
import matplotlib.pyplot as plt
import numpy as np
from tabulate import tabulate
from ADS import ABM4
def exact(x):
return (np.exp(-3 * x) * (-129 * (x ** 4) - 16 * (x ** 3) + 54 * (x ** 2) + 36 * x)) / 12
def f(u, x):
return (u[1],
u[2],
u[3],
u[4],
- 15 * u[4] - 90 * u[3] -
270 * u[2] - 405 * u[1] - 243 * u[0])
y0 = [0, 3, -9, -8, 0]
solver = ABM4(f)
solver.set_initial_condition(y0)
x_points = np.linspace(0, 5, 50)
u, x = solver.solve(x_points)
y = u[:, 0]
plt.plot(x, y, 'bo', linewidth=3, markersize=3)
plt.plot(x, exact(x), 'r-')
plt.grid()
plt.show()
I get the same error:
Traceback (most recent call last):
File "C:\Fin_Proj_ODE\test1.py", line 29, in <module>
u, x = solver.solve(x_points)
File "C:\Fin_Proj_ODE\ODS.py", line 71, in solve
self.u[k + 1] = self.solver_st()
File "C:\Fin_Proj_ODE\ADS.py", line 22, in solver_st
K1 = dx * f(u[k], x[k])
File "C:\Fin_Proj_ODE\ODS.py", line 20, in <lambda>
self.f = lambda u, x: np.asarray(f(u, x), float)
File "C:\Fin_Proj_ODE\test1.py", line 17, in f
return (u[1],
IndexError: invalid index to scalar variable.
Process finished with exit code 1
How to improve code and fix this bug?
Probably, this bug tied with uncorrect loop implementation in #ADS.py, so the main loop iteration in ODE_Solver compute arrays of solution throughout the range coord_points. And so,now I don't understnd how to rewrite this loop in ABM4 class.
