2014 dxdy logo

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

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




Начать новую тему Ответить на тему
 
 Решение задачи Коши различными численными методами
Сообщение20.05.2022, 22:48 


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

$\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 


21/05/16
4292
Аделаида
Уберите @classmethod. Иначе solver_st принимает не инстанс класса, а сам класс.

 Профиль  
                  
 
 Re: Решение задачи Коши различными численными методами
Сообщение21.05.2022, 11:34 


17/03/20
183
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 


21/05/16
4292
Аделаида
Alm99 в сообщении #1555031 писал(а):
Например, при реализации схем при вызове AM4 возникает ошибка:

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

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

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

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 4 ] 

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group