2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




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

$\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
  1. #ODS.py
  2.  
  3. import numpy as np
  4.  
  5.  
  6. class ODE_Solver(object):
  7.     """
  8.    Суперкласс для решения ОДУ или системы ОДУ нормированной и приведенной к виду: du/dx = f(u, x)
  9.    Атрибуты класса:
  10.    x: массив узлов точек координаты x
  11.    u: массив решения ОДУ в точках узла x
  12.    k: число шагов вычисления значения в узлах
  13.    f: функция правой части ОДУ, реализованная в виде: f(u, x)
  14.    """
  15.  
  16.  
  17.     def __init__(self, f):
  18.         if not callable(f):
  19.             # проверка корректности типа передаваемой функции f(u, x)
  20.             raise TypeError('f не %s, является функцией' % type(f))
  21.         # результат работы вычисления функции решателя (тип float)
  22.         self.f = lambda u, x: np.asarray(f(u, x), float)
  23.  
  24.     def solver_st(self):
  25.         """Реализация численной схемы решателя"""
  26.         raise NotImplementedError
  27.  
  28.     def set_initial_condition(self, u0):
  29.         if isinstance(u0, (float, int)):  # ОДУ является одномерным
  30.             self.neq = 1
  31.             u0 = float(u0)
  32.         else:  # система ОДУ, ОДУ порядка выше 1-го
  33.             u0 = np.asarray(u0)  # (начальные условия вектор-функция - порядок 2+)
  34.             self.neq = u0.size
  35.         self.u0 = u0
  36.  
  37.         # Проверка, что функция возвращает вектор f корректной длины:
  38.         try:
  39.             f0 = self.f(self.u0, 0)
  40.         except IndexError:
  41.             raise IndexError(
  42.                 'Индекс вектора u выходит за границы f(u,x). Допустимые индексы %s' % (str(range(self.neq))))
  43.         if f0.size != self.neq:
  44.             raise ValueError('f(u,x) возвращено %d элементов, вектор u имеет %d элементов' % (f0.size, self.neq))
  45.  
  46.     def solve(self, coord_points, terminate=None):
  47.         """
  48.        Решение ОДУ и  запись полученных значений в виде массива или списка.
  49.        По умолчанию функция возвращает False, если нет элементов.
  50.        """
  51.         if terminate is None:
  52.             terminate = lambda u, x, step_no: False
  53.  
  54.         if isinstance(coord_points, (float, int)):
  55.             raise TypeError('solve: массив точек x не является итерируемым объектом')
  56.         self.x = np.asarray(coord_points)
  57.         if self.x.size <= 1:
  58.             raise ValueError('ODESolver.solve требует на вход массив координат'
  59.                              ' точек поиска решения. Число точек меньше 2!')
  60.  
  61.         n = self.x.size
  62.         if self.neq == 1:  # ОДУ
  63.             self.u = np.zeros(n)
  64.         else:  # ОДУ порядка 2+ или система ОДУ
  65.             self.u = np.zeros((n, self.neq))
  66.  
  67.         # Присвоить self.x[0] corresponds to self.u0
  68.         self.u[0] = self.u0
  69.  
  70.         # Проход по сетке координат x (использование векторизации)
  71.         for k in range(n - 1):
  72.             self.k = k
  73.             self.u[k + 1] = self.solver_st()
  74.             if terminate(self.u, self.x, self.k + 1):
  75.                 break  # прервать цикл при последнем k
  76.         return self.u[:k + 2], self.x[:k + 2]


Непосредственно реализация класса численной схемы Эйлера FE:
код: [ скачать ] [ спрятать ]
Используется синтаксис Python
  1.  
  2. #ES.py
  3.  
  4. from abc import ABC
  5. from ODS import ODE_Solver
  6.  
  7.  
  8. class FE(ODE_Solver, ABC):
  9.     """
  10.    Реализация классической прямой схемы Эйлера первого порядка точности.
  11.    Возвращает массив решения ОДУ или системы ОДУ.
  12.    Атрибуты класса:
  13.    x: массив узлов точек координаты x
  14.    u: массив решения ОДУ в точках узла x
  15.    k: число шагов вычисления значения в узлах
  16.    f: функция правой части ОДУ, реализованная в виде: f(u, x)
  17.    """
  18.  
  19.  
  20.     @classmethod
  21.     def solver_st(self):
  22.         u, f, k, x = self.u, self.f, self.k, self.x
  23.         dx = x[k + 1] - x[k]
  24.         u_new = u[k] + dx * f(u[k], x[k])
  25.         return u_new
  26.  


Непосрественно пытаюсь вызвать теперь схему решателя и передать в нее ОДУ пятого порядка def f(u,x):
код: [ скачать ] [ спрятать ]
Используется синтаксис Python
  1. #test.py
  2.  
  3. import matplotlib.pyplot as plt
  4. import numpy as np
  5. from ES import FE
  6.  
  7.  
  8. def f(u, x):
  9.     return (u[1],
  10.             u[2],
  11.             u[3],
  12.             u[4],
  13.             - 15 * u[4] - 90 * u[3] -
  14.             270 * u[2] - 405 * u[1] - 243 * u[0])
  15.  
  16.  
  17. y0 = [0, 3, -9, -8, 0]
  18. solver = FE(f)
  19. solver.set_initial_condition(y0)
  20. x_points = np.linspace(0, 5, 150)
  21. u, x = solver.solve(x_points)
  22. plt.plot(x, u)
  23.  


При вызове требуемой схемы решателя я получаю странную ошибку следующего типа:
  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
  1. #RKS.py
  2.  
  3. from abc import ABC
  4. from ODS import ODE_Solver
  5.  
  6.  
  7. class RK4(ODE_Solver, ABC):
  8.     """
  9.    Реализация классической схемы Рунге-Кутта 4 порядка. Является наиболее стабильной
  10.    неадаптивной многошаговой схемой 4 порядка точности. Явная реализация схемы.
  11.    Возвращает массив решения ОДУ или системы ОДУ.
  12.    Атрибуты класса:
  13.    x: массив узлов точек координаты x
  14.    u: массив решения ОДУ в точках узла x
  15.    k: число шагов вычисления значения в узлах
  16.    f: функция правой части ОДУ, реализованная в виде: f(u, x)
  17.    """
  18.  
  19.     def solver_st(self):
  20.         u, f, k, x = self.u, self.f, self.k, self.x
  21.         dx = x[k + 1] - x[k]
  22.         dx2 = dx / 2.0
  23.         K1 = dx * f(u[k], x[k])
  24.         K2 = dx * f(u[k] + 0.5 * K1, x[k] + dx2)
  25.         K3 = dx * f(u[k] + 0.5 * K2, x[k] + dx2)
  26.         K4 = dx * f(u[k] + K3, x[k] + dx)
  27.         u_new = u[k] + (1 / 6.0) * (K1 + 2 * K2 + 2 * K3 + K4)
  28.         return u_new
  29.  


код: [ скачать ] [ спрятать ]
Используется синтаксис Python
  1. #RKFS.py
  2.  
  3. from abc import ABC
  4. from ODS import ODE_Solver
  5. import numpy as np
  6.  
  7.  
  8. class RKF(ODE_Solver, ABC):
  9.     """
  10.    Реализация  Рунге-Кутта-Фельберга 4 порядка точности (модификация метода Рунге-Кутта с адаптивным шагом)
  11.    Возвращает массив решения ОДУ или системы ОДУ.
  12.    Атрибуты класса:
  13.    x: массив узлов точек координаты x
  14.    u: массив решения ОДУ в точках узла x
  15.    k: число шагов вычисления значения в узлах
  16.    f: функция правой части ОДУ, реализованная в виде: f(u, x)
  17.    """
  18.  
  19.     def __init__(self, f):
  20.         super().__init__(f)
  21.         self.tol = None
  22.         self.coeff_err = None
  23.         self.beta = None
  24.         self.coeff = None
  25.         self.alpha = None
  26.         self.coeff_star = None
  27.         self.coeff_err = None
  28.         self.alpha = None
  29.  
  30.     def solver_st(self):
  31.         u, f, k, x, alpha, beta, coeff, coeff_star, coeff_err, tol = self.u, self.f, self.k, self.alpha, self.x, \
  32.                                                                      self.beta, self.coeff, self.coeff_star, \
  33.                                                                      self.tol, self.coeff_err
  34.  
  35.         # Таблица Бутчера коэффициентов схемы Рунге-Кутта 4 порядка и 5 порядка точности.
  36.         # alpha = np.array([0.0, 1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0])
  37.         beta = np.array([[0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
  38.                          [1.0 / 4.0, 0.0, 0.0, 0.0, 0.0, 0.0],
  39.                          [3.0 / 32.0, 9.0 / 32.0, 0.0, 0.0, 0.0, 0.0],
  40.                          [1932.0 / 2197.0, (-7200.0) / 2197.0, 7296.0 / 2197.0, 0.0, 0.0, 0.0],
  41.                          [439.0 / 216.0, -8.0, 3680.0 / 513.0, (-845.0) / 4104.0, 0.0, 0.0],
  42.                          [(-8.0) / 27.0, 2.0, (-3544.0) / 2565.0, 1859.0 / 4104.0, (-11.0) / 40.0, 0.0]])
  43.         coeff = np.array(
  44.             [25.0 / 216.0, 0.0, 1408.0 / 2565.0, 2197.0 / 4104.0, (-1.0) / 5.0, 0.0])
  45.         # коэффициенты 4 порядка точности
  46.         coeff_star = np.array([16.0 / 135.0, 0.0, 6656.0 / 12825.0, 28561.0 / 56430.0, (-9.0) / 50.0,
  47.                                2.0 / 55.0])  # коэффициенты 5 порядка точности
  48.  
  49.         coeff_err = coeff - coeff_star  # разность коэффициентов методов 4 и 5 порядка
  50.         dx = x[k + 1] - x[k]  # начальное приближение адаптивного шага решателя
  51.  
  52.         K1 = dx * f(u[k], x[k])
  53.         K2 = f(u[k] + dx * beta[1, 0] * K1)
  54.         K3 = f(u[k] + dx * (beta[2, 0] * K1 + beta[2, 1] * K2))
  55.         K4 = f(u[k] + dx * (beta[3, 0] * K1 + beta[3, 1] * K2 + beta[3, 2] * K3))
  56.         K5 = f(u[k] + dx * (beta[4, 0] * K1 + beta[4, 1] * K2 + beta[4, 2] * K3 + beta[4, 3] * K4))
  57.         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))
  58.         errorfield = dx * (coeff_err[0] * K1 + coeff_err[1] * K2 + coeff_err[2] * K3 + coeff_err[3] * K4 +
  59.                            coeff_err[4] * K5 + coeff_err[5] * K6)
  60.         max_error = np.absolute(errorfield).max()
  61.  
  62.         if max_error <= tol:
  63.             u[k + 1] = u[k] + dx * (coeff_star[0] * K1 + coeff_star[1] * K2 + coeff_star[2] * K3 +
  64.                                     coeff_star[3] * K4 + coeff_star[4] * K5 + coeff_star[5] * K6)
  65.             dxx = dx * (tol / max_error) ** 0.2
  66.             dx = dxx
  67.             x[k + 1] = x[k] + dx
  68.             eps = max_error
  69.         else:
  70.             dx = dx * (tol / max_error) ** 0.25
  71.         u_new = u[k]
  72.         return u_new, eps
  73.  


код: [ скачать ] [ спрятать ]
Используется синтаксис Python
  1. #ADS.py
  2.  
  3. from abc import ABC
  4. from ODS import ODE_Solver
  5.  
  6.  
  7. class ABF4(ODE_Solver, ABC):
  8.     """
  9.    Реализация классической 4-х шаговой схемы Адамса-Башфорта. Явная реализация схемы.
  10.    Возвращает массив решения ОДУ или системы ОДУ.
  11.    Атрибуты класса:
  12.    x: массив узлов точек координаты x
  13.    u: массив решения ОДУ в точках узла x
  14.    k: число шагов вычисления значения в узлах
  15.    f: функция правой части ОДУ, реализованная в виде: f(u, x)
  16.    """
  17.  
  18.     def solver_st(self):
  19.         u, f, k, x = self.u, self.f, self.k, self.x
  20.         dx = x[k + 1] - x[k]
  21.         dx2 = dx / 2.0
  22.         K1 = dx * f(u[k], x[k])
  23.         K2 = dx * f(u[k] + 0.5 * K1, x[k] + dx2)
  24.         K3 = dx * f(u[k] + 0.5 * K2, x[k] + dx2)
  25.         K4 = dx * f(u[k] + K3, x[k] + dx)
  26.         u = u[k] + (1 / 6.0) * (K1 + 2 * K2 + 2 * K3 + K4)
  27.  
  28.         u_new = (u[-1] + (dx / 24) * (55 * f(x[-1], u[-1]) - 59 * f(x[-2], u[-2])
  29.                                       + 37 * f(x[-3], u[-3]) - 9 * f(x[-4], u[-4])))
  30.         # x = (x[-1] + dx)
  31.         return u_new
  32.  


код: [ скачать ] [ спрятать ]
Используется синтаксис Python
  1. #AMS.py
  2.  
  3. from abc import ABC
  4. from ODS import ODE_Solver
  5.  
  6.  
  7. class AM4(ODE_Solver, ABC):
  8.     """
  9.    Реализация классической 4-х шаговой схемы Адамса-Моултона. Явная реализация схемы.
  10.    Возвращает массив решения ОДУ или системы ОДУ.
  11.    Атрибуты класса:
  12.    x: массив узлов точек координаты x
  13.    u: массив решения ОДУ в точках узла x
  14.    k: число шагов вычисления значения в узлах
  15.    f: функция правой части ОДУ, реализованная в виде: f(u, x)
  16.    """
  17.  
  18.     def solver_st(self):
  19.         u, f, k, x = self.u, self.f, self.k, self.x
  20.         dx = x[k + 1] - x[k]
  21.         dx2 = dx / 2.0
  22.         K1 = dx * f(u[k], x[k])
  23.         K2 = dx * f(u[k] + 0.5 * K1, x[k] + dx2)
  24.         K3 = dx * f(u[k] + 0.5 * K2, x[k] + dx2)
  25.         K4 = dx * f(u[k] + K3, x[k] + dx)
  26.         u = u[k] + (1 / 6.0) * (K1 + 2 * K2 + 2 * K3 + K4)
  27.  
  28.         u_new = (u[-1] + (dx / 720) * (251 * f(x[-1] + dx) + 646 * f(x[-1], u[-1])
  29.                                        - 264 * f(x[-2], u[-2]) + 106 * f(x[-3], u[-3])
  30.                                        - 19 * f(x[-4], u[-4])))
  31.         # x = (x[-1] + dx)
  32.         return u_new
  33.  



Аналитическое решение данной задачи Коши:
$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}$


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

 
 
 
 Re: Решение задачи Коши различными численными методами
Сообщение21.05.2022, 05:43 
Уберите @classmethod. Иначе solver_st принимает не инстанс класса, а сам класс.

 
 
 
 Re: Решение задачи Коши различными численными методами
Сообщение21.05.2022, 11:34 
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, не хочется его делать атрибутом всех классов, хотя это будет верно и правильно с точки зрения вычислений.

 
 
 
 Re: Решение задачи Коши различными численными методами
Сообщение21.05.2022, 13:08 
Alm99 в сообщении #1555031 писал(а):
Например, при реализации схем при вызове AM4 возникает ошибка:

Потому что ODE_Solver.f принимает два аргумента, а не один.
Alm99 в сообщении #1555031 писал(а):
а при вызове ABF4

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

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

 
 
 [ Сообщений: 4 ] 


Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group