Индивидуальное задание (Игорь Фатеев)
Решение задачи Коши при помощи методов Рунге-Кутты
Рассматривается задача Коши вида
{u′=f(x,u),u(x0)=u0,
где u - вектор неизвестных функций, x - переменная, f - набор вектор-функций нескольких переменных.
Методы Рунге-Кутты задаются формулами:u(x+h)=u(x)+q∑i=1pikiгде введено обозначение
k1=hf(x,u(x))k2=hf(x+α2h,u(x)+β21k1).....................................kq=hf(x+αqh,u(x)+βq1k1+...+βq,q−1kq−1)=hf(x+αqh,u(x)+q−1∑i=1βqiki)
Конкретный метод определяется числом q и коэффициентами αi и βij. Эти коэффициенты часто упорядочивают в таблицу (называемую таблицей Батчера):
0α2β21α3β31β32⋮⋮⋮⋱αqβq1βq2…βq,q−1p1p2…pq−1pq
Вложенные формулы Рунге-Кутты
Используются следующие формулы
u(x+h)=u(x)+q∑i=1piki
˜u(x+h)=u(x)+q∑i=1˜piki
Формула для u(x+h) имеет погрешность порядка p, Формула для ˜u(x+h) имеет погрешность s=p±1. Для практической оценки погрешности получается выражение
r(x+h)=u(x+h)−˜u(x+h)≈q∑i=1(pi−˜pi)ki
Таблица Батчера приобретает вид:
0α2β21α3β31β32⋮⋮⋮⋱αqβq1βq2…βq,q−1p1p2…pq−1pq˜p1˜p2…˜pq−1˜pq
Решение задачи Коши средствами языка C#
- Для решения задачи создаём новый проект консольного приложения под названием RungeKutta и добавляем в проект новый обобщенный интерфейс (в отдельном файле).
/// <summary> /// Интерфейс для решения дифференциальных уравнений /// </summary> /// <typeparam name="T">тип переменной</typeparam> /// <typeparam name="U">тип значения функции</typeparam> interface IOdEsolver<T, U> { /// <summary> /// правые части дифференциальных уравнений /// </summary> List<Func<T, List<U>, U>> Equations { get; } // List<T>,
Func<
in
T1,
in
T2,
out
TResult>
/// <summary> /// погрешность вычислений /// </summary> double Epsilon { get; } /// <summary> /// метод для решения задачи Коши /// </summary> /// <param name="a">точка, в которой заданы начальные условия</param> /// <param name="b">точка, в которой ищется решение</param> /// <param name="initialConditions">вектор начальных условий</param> /// <returns>решение задачи Коши в точке x=b</returns> List<T> Solve(T a, T b, List<T> initialConditions); } - Добавляем в проект статический класс ButcherTableau, содержащий таблицы Батчера для нескольких методов из семейства методов Рунге-Кутты: метода Хойна, метода Богацки-Шампайна, метода Кеша-Карпа, метода Рунге-Кутты-Фельдберга и двух методов Дорманда-Принса разных порядков точности.
- Добавляем в проект класс, реализующий интерфейс
IOdEsolver
public class OdeSolver : IOdEsolver<double, double>
- Добавляем в новый класс конструктор
public OdeSolver(List<Func<double, List<double>, double>> funcs, double[][] butcherTableau, double epsilon = 0.1e-6)
- Реализуем метод Solve. Алгоритм его работы следующий
- Вычисляем новую переменную h = b - a
- В цикле, пока a<b (используем условие Math.Abs(b - a) > double.Epsilon)
- Вычисляем u(a+h) и ˜u(a+h).
- Вычисляем погрешность r(a+h)
- Если модуль погрешность превышает Epsilon, делим h пополам и снова вычисляем u(a+h) и ˜u(a+h).
- a += h;
- Вектор начальных условий приравниваем ˜u(a+h)
- Возвращаем ˜u(a+h).
- Создаём модульные тесты. Код тестирующего метода должен быть примерно таким
[TestMethod()] public void SolveTest() { var equation = new List<Func<double, List<double>, double>> { (x, y) => y[0] }; OdeSolver odeSolver = new OdeSolver(equation, ButcherTableau.RungeKuttaFeldberg, 0.1e-6); var solution = odeSolver.Solve(0, 1, new List<double> { 1 }); Assert.IsTrue(Math.Abs(solution[0] - Math.E) < 0.1e-5); var equations = new List<Func<double, List<double>, double>> {(x, y) => y[1], (x, y) => -y[0]}; OdeSolver solver = new OdeSolver(equations, ButcherTableau.RungeKuttaFeldberg, 0.1e-6); solution = solver.Solve(0, 1, new List<double> {0, 1}); Assert.IsTrue(Math.Abs(solution[0] - Math.Sin(1)) < 0.1e-5); Assert.IsTrue(Math.Abs(solution[1] - Math.Cos(1)) < 0.1e-5); // добавьте свои тесты для каждого из методов }
- Продемонстрируйте использование класса OdeSolver в основной программе
9 июня 2019, 16:55 |