Решение задачи Коши различными численными методами

[Alm99]

Добрый вечер, уважаемые форумчане. Сложно было решить, в каком разделе оставить данную тему, но я думаю все больше в этом разделе. Перейду непосредственно к самой задаче. Дано ОДУ, необходимо решить задачу Коши:

$\frac{d^{5}y}{dx^{5}} + 15 \frac{d^{4}y}{dx^{4}} + 90 \frac{d^{3}y}{dx^{3}} + 270 \frac{d^{2}y}{dx^{2}} + 405 \frac{dy}{dx} + 243y = 0$

$y\Bigg|_{x=0} = 0; \quad \frac{dy}{dx}\Bigg|_{x=0} = 3; \quad \frac{d^{2}y}{dx^{2}}\Bigg|_{x=0} = -9; \quad \frac{d^{3}y}{dx^{3}}\Bigg|_{x=0} = -8; \quad \frac{d^{4}y}{dx^{4}}\Bigg|_{x=0} = 0;$

Пробую написать так cказать программу решатель ОДУ 5 порядка, используя принципы ООП. Ноебходимо реализовать несколько численных схем решения, чтобы потом можно было решать ОДУ и сравнить результаты решения.

Каждая из численных схем реализована в виде отдельного класса, который наследует методы суперкласса ODE_Solver - по сути в нем осуществляется проверка, что передаваемая функция (правая часть уравнения является вектором или скаляром, также прописан сам цикл прохода по узлам, в которых отыскивается решение ОДУ. Каждя численная схема реализована в виде класса, при этом, содежит только метод solve_st. В качетсве примера рассмотрим реализацию прямой схемы Эйлера для решения ОДУ 5 степени с начальными условиями y0.

Ниже представлена реализация класса ODE_Solver

Используется синтаксис Python

#ODS.py
 
import numpy as np
 
 
class ODE_Solver(object):
    """
    Суперкласс для решения ОДУ или системы ОДУ нормированной и приведенной к виду: du/dx = f(u, x)
    Атрибуты класса:
    x: массив узлов точек координаты x
    u: массив решения ОДУ в точках узла x
    k: число шагов вычисления значения в узлах
    f: функция правой части ОДУ, реализованная в виде: f(u, x)
    """
 
 
    def __init__(self, f):
        if not callable(f):
            # проверка корректности типа передаваемой функции f(u, x)
            raise TypeError('f не %s, является функцией' % type(f))
        # результат работы вычисления функции решателя (тип float)
        self.f = lambda u, x: np.asarray(f(u, x), float)
 
    def solver_st(self):
        """Реализация численной схемы решателя"""
        raise NotImplementedError
 
    def set_initial_condition(self, u0):
        if isinstance(u0, (float, int)):  # ОДУ является одномерным
            self.neq = 1
            u0 = float(u0)
        else:  # система ОДУ, ОДУ порядка выше 1-го
            u0 = np.asarray(u0)  # (начальные условия вектор-функция - порядок 2+)
            self.neq = u0.size
        self.u0 = u0
 
        # Проверка, что функция возвращает вектор f корректной длины:
        try:
            f0 = self.f(self.u0, 0)
        except IndexError:
            raise IndexError(
                'Индекс вектора u выходит за границы f(u,x). Допустимые индексы %s' % (str(range(self.neq))))
        if f0.size != self.neq:
            raise ValueError('f(u,x) возвращено %d элементов, вектор u имеет %d элементов' % (f0.size, self.neq))
 
    def solve(self, coord_points, terminate=None):
        """
        Решение ОДУ и  запись полученных значений в виде массива или списка.
        По умолчанию функция возвращает False, если нет элементов.
        """
        if terminate is None:
            terminate = lambda u, x, step_no: False
 
        if isinstance(coord_points, (float, int)):
            raise TypeError('solve: массив точек x не является итерируемым объектом')
        self.x = np.asarray(coord_points)
        if self.x.size <= 1:
            raise ValueError('ODESolver.solve требует на вход массив координат'
                             ' точек поиска решения. Число точек меньше 2!')
 
        n = self.x.size
        if self.neq == 1:  # ОДУ
            self.u = np.zeros(n)
        else:  # ОДУ порядка 2+ или система ОДУ
            self.u = np.zeros((n, self.neq))
 
        # Присвоить self.x[0] corresponds to self.u0
        self.u[0] = self.u0
 
        # Проход по сетке координат x (использование векторизации)
        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  # прервать цикл при последнем k
        return self.u[:k + 2], self.x[:k + 2] 

Непосредственно реализация класса численной схемы Эйлера FE:

Используется синтаксис Python

 
#ES.py
 
from abc import ABC
from ODS import ODE_Solver
 
 
class FE(ODE_Solver, ABC):
    """
    Реализация классической прямой схемы Эйлера первого порядка точности.
    Возвращает массив решения ОДУ или системы ОДУ.
    Атрибуты класса:
    x: массив узлов точек координаты x
    u: массив решения ОДУ в точках узла x
    k: число шагов вычисления значения в узлах
    f: функция правой части ОДУ, реализованная в виде: f(u, x)
    """
 
 
    @classmethod
    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
 

Непосрественно пытаюсь вызвать теперь схему решателя и передать в нее ОДУ пятого порядка def f(u,x):

Используется синтаксис Python

#test.py
 
import matplotlib.pyplot as plt
import numpy as np
from ES import FE
 
 
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 = FE(f)
solver.set_initial_condition(y0)
x_points = np.linspace(0, 5, 150)
u, x = solver.solve(x_points)
plt.plot(x, u)
 

При вызове требуемой схемы решателя я получаю странную ошибку следующего типа:

1. File "C:\Fin_Proj_ODE\test1.py", line 19, in <module> 
2. u, x = solver.solve(x_points) 
3. File "C:\Fin_Proj_ODE\ODS.py", line 71, in solve 
4. self.u[k + 1] = self.solver_st() 
5. File "C:\Fin_Proj_ODE\ES.py", line 18, in solver_st 
6. u, f, k, x = self.u, self.f, self.k, self.x 
7. AttributeError: type object 'FE' has no attribute 'u' 

Уважаемые форумчане, мне очень нужна Ваша помощь и возникло несколько вопросов:

1) Каким образом решить подобную ошибку, потому что я же наследовал класс численной схемы от класса ODE_Solver, тогад почему при вызове конкретного решателя возникает такая проблема?

2) Каким образом можно реализовать класс Problem с исходной задачей и передавать требуемые аргументы, а также вернуть массив погрешности работы решателя (погрешность - разность численного расчета и аналитического решения)? Про погрешность метода в данном случае речи не идет и получения ее оценки - она известна.

3) Насколько правильно в данном случае реализованы остальные методы решателя? Здесь также хотел бы услышать ценные замечания а нужны ли данные методы в рамках решения представленной задачи, или простой схемы Эйлера, Рунге-Кутты 4 порядка достаточно?

Используется синтаксис Python

#RKS.py
 
from abc import ABC
from ODS import ODE_Solver
 
 
class RK4(ODE_Solver, ABC):
    """
    Реализация классической схемы Рунге-Кутта 4 порядка. Является наиболее стабильной
    неадаптивной многошаговой схемой 4 порядка точности. Явная реализация схемы.
    Возвращает массив решения ОДУ или системы ОДУ.
    Атрибуты класса:
    x: массив узлов точек координаты x
    u: массив решения ОДУ в точках узла x
    k: число шагов вычисления значения в узлах
    f: функция правой части ОДУ, реализованная в виде: 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]
        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_new = u[k] + (1 / 6.0) * (K1 + 2 * K2 + 2 * K3 + K4)
        return u_new
 

Используется синтаксис Python

#RKFS.py
 
from abc import ABC
from ODS import ODE_Solver
import numpy as np
 
 
class RKF(ODE_Solver, ABC):
    """
    Реализация  Рунге-Кутта-Фельберга 4 порядка точности (модификация метода Рунге-Кутта с адаптивным шагом)
    Возвращает массив решения ОДУ или системы ОДУ.
    Атрибуты класса:
    x: массив узлов точек координаты x
    u: массив решения ОДУ в точках узла x
    k: число шагов вычисления значения в узлах
    f: функция правой части ОДУ, реализованная в виде: f(u, x)
    """
 
    def __init__(self, f):
        super().__init__(f)
        self.tol = None
        self.coeff_err = None
        self.beta = None
        self.coeff = None
        self.alpha = None
        self.coeff_star = None
        self.coeff_err = None
        self.alpha = None
 
    def solver_st(self):
        u, f, k, x, alpha, beta, coeff, coeff_star, coeff_err, tol = self.u, self.f, self.k, self.alpha, self.x, \
                                                                     self.beta, self.coeff, self.coeff_star, \
                                                                     self.tol, self.coeff_err
 
        # Таблица Бутчера коэффициентов схемы Рунге-Кутта 4 порядка и 5 порядка точности.
        # alpha = np.array([0.0, 1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0])
        beta = np.array([[0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [1.0 / 4.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [3.0 / 32.0, 9.0 / 32.0, 0.0, 0.0, 0.0, 0.0],
                         [1932.0 / 2197.0, (-7200.0) / 2197.0, 7296.0 / 2197.0, 0.0, 0.0, 0.0],
                         [439.0 / 216.0, -8.0, 3680.0 / 513.0, (-845.0) / 4104.0, 0.0, 0.0],
                         [(-8.0) / 27.0, 2.0, (-3544.0) / 2565.0, 1859.0 / 4104.0, (-11.0) / 40.0, 0.0]])
        coeff = np.array(
            [25.0 / 216.0, 0.0, 1408.0 / 2565.0, 2197.0 / 4104.0, (-1.0) / 5.0, 0.0])
        # коэффициенты 4 порядка точности
        coeff_star = np.array([16.0 / 135.0, 0.0, 6656.0 / 12825.0, 28561.0 / 56430.0, (-9.0) / 50.0,
                               2.0 / 55.0])  # коэффициенты 5 порядка точности
 
        coeff_err = coeff - coeff_star  # разность коэффициентов методов 4 и 5 порядка
        dx = x[k + 1] - x[k]  # начальное приближение адаптивного шага решателя
 
        K1 = dx * f(u[k], x[k])
        K2 = f(u[k] + dx * beta[1, 0] * K1)
        K3 = f(u[k] + dx * (beta[2, 0] * K1 + beta[2, 1] * K2))
        K4 = f(u[k] + dx * (beta[3, 0] * K1 + beta[3, 1] * K2 + beta[3, 2] * K3))
        K5 = f(u[k] + dx * (beta[4, 0] * K1 + beta[4, 1] * K2 + beta[4, 2] * K3 + beta[4, 3] * K4))
        K6 = f(u[k] + dx * (beta[5, 0] * K1 + beta[5, 1] * K2 + beta[5, 2] * K3 + beta[5, 3] * K4 + beta[5, 4] * K5))
        errorfield = dx * (coeff_err[0] * K1 + coeff_err[1] * K2 + coeff_err[2] * K3 + coeff_err[3] * K4 +
                           coeff_err[4] * K5 + coeff_err[5] * K6)
        max_error = np.absolute(errorfield).max()
 
        if max_error <= tol:
            u[k + 1] = u[k] + dx * (coeff_star[0] * K1 + coeff_star[1] * K2 + coeff_star[2] * K3 +
                                    coeff_star[3] * K4 + coeff_star[4] * K5 + coeff_star[5] * K6)
            dxx = dx * (tol / max_error) ** 0.2
            dx = dxx
            x[k + 1] = x[k] + dx
            eps = max_error
        else:
            dx = dx * (tol / max_error) ** 0.25
        u_new = u[k]
        return u_new, eps
 

Используется синтаксис Python

#ADS.py
 
from abc import ABC
from ODS import ODE_Solver
 
 
class ABF4(ODE_Solver, ABC):
    """
    Реализация классической 4-х шаговой схемы Адамса-Башфорта. Явная реализация схемы.
    Возвращает массив решения ОДУ или системы ОДУ.
    Атрибуты класса:
    x: массив узлов точек координаты x
    u: массив решения ОДУ в точках узла x
    k: число шагов вычисления значения в узлах
    f: функция правой части ОДУ, реализованная в виде: 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]
        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)
 
        u_new = (u[-1] + (dx / 24) * (55 * f(x[-1], u[-1]) - 59 * f(x[-2], u[-2])
                                      + 37 * f(x[-3], u[-3]) - 9 * f(x[-4], u[-4])))
        # x = (x[-1] + dx)
        return u_new
 

Используется синтаксис Python

#AMS.py
 
from abc import ABC
from ODS import ODE_Solver
 
 
class AM4(ODE_Solver, ABC):
    """
    Реализация классической 4-х шаговой схемы Адамса-Моултона. Явная реализация схемы.
    Возвращает массив решения ОДУ или системы ОДУ.
    Атрибуты класса:
    x: массив узлов точек координаты x
    u: массив решения ОДУ в точках узла x
    k: число шагов вычисления значения в узлах
    f: функция правой части ОДУ, реализованная в виде: 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]
        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)
 
        u_new = (u[-1] + (dx / 720) * (251 * f(x[-1] + dx) + 646 * f(x[-1], u[-1])
                                       - 264 * f(x[-2], u[-2]) + 106 * f(x[-3], u[-3])
                                       - 19 * f(x[-4], u[-4])))
        # x = (x[-1] + dx)
        return u_new
 

Аналитическое решение данной задачи Коши:

$y(t)\ = \frac{\mathrm{e}^{-3\thinsp t}\thinsp\left(-129\thinsp t^4-16\thinsp t^3+54\thinsp t^2+36\thinsp t\right)}{12}$

Может часть методов и вовсе не даст существенных преимуществ решения достаточно простой задачи Коши? Заранее, огромное спасибо за внимание и терпение!

[kotenok gav]

Уберите @classmethod. Иначе solver_st принимает не инстанс класса, а сам класс.

[Alm99]

kotenok gav
Это действительно разрешило проблему, спасибо огромное! Да, попутно возник вопрос отноистельно реализации численных схем. Все же имеет ли смысл реализовать численные схемы Адамса-Башфорта, Адамса-Моултона и Рунге-Кутта модифицированную? Например, при реализации схем при вызове AM4 возникает ошибка:

Используется синтаксис Python

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 = AM4(f)  # call AM4
# solver = ABF4(f)  # call ABF4
solver.set_initial_condition(y0)
x_points = np.linspace(0, 5, 100)
u, x = solver.solve(x_points)
y = u[:, 0]
 

1. File "C:\Fin_Proj_ODE\test1.py", line 23, in <module> 
2. u, x = solver.solve(x_points) 
3. File "C:\Fin_Proj_ODE\ODS.py", line 71, in solve 
4. self.u[k + 1] = self.solver_st() 
5. File "C:\Fin_Proj_ODE\AMS.py", line 26, in solver_st 
6. u_new = (u[k-1] + (dx / 720) * (251 * f(x[k-1] + dx) + 646 * f(x[k-1], u[k-1]) 
7. TypeError: ODE_Solver.__init__.<locals>.<lambda>() missing 1 required positional argument: 'x' 

а при вызове ABF4

1. Traceback (most recent call last): 
2. File "C:\Fin_Proj_ODE\test1.py", line 23, in <module> 
3. u, x = solver.solve(x_points) 
4. File "C:\Fin_Proj_ODE\ODS.py", line 71, in solve 
5. self.u[k + 1] = self.solver_st() 
6. File "C:\Fin_Proj_ODE\ADS.py", line 26, in solver_st 
7. u_new = (u[k-1] + (dx / 24) * (55 * f(x[k-1], u[k-1]) - 59 * f(x[k-2], u[k-2]) 
8. File "C:\Fin_Proj_ODE\ODS.py", line 20, in <lambda> 
9. self.f = lambda u, x: np.asarray(f(u, x), float) 
10. File "C:\Fin_Proj_ODE\test1.py", line 11, in f 
11. return (u[1], 
12. IndexError: invalid index to scalar variable. 
13.  

Их реализация почти та же, что и у RK, который работает правильно, почему тогда при обновлении уточненного значения u_new не удается вернуть массив самого решения и производных? Согласен, что схема RKF реализована неверно, потому что требуется передавать дополнительно параметр tol, не хочется его делать атрибутом всех классов, хотя это будет верно и правильно с точки зрения вычислений.

[kotenok gav]

Alm99 в сообщении #1555031 писал(а):

Например, при реализации схем при вызове AM4 возникает ошибка:

Потому что ODE_Solver.f принимает два аргумента, а не один.

Alm99 в сообщении #1555031 писал(а):

а при вызове ABF4

Потому что вы передаёте ODE_Solver.f один численный аргумент, а их должно быть два - и первый из них должен быть вектором (насчёт второго не понял).

И неужели так сложно находить эти ошибки вам самим? Достаточно просто внимательно читать, что выдаёт вам Питон, и смотреть сам ваш код.