Приближенное решение задачи Коши на C#

Решением задачи Коши для дифференциального уравнения порядка m называется такая функция y(x), которая удовлетворяет m начальным условиям. Исходя из названия поста и вышесказанного определения, можно догадаться, что сегодня мы будем писать программу, а точнее класс, на C#, в котором будут описаны функции вычисления приближенного значения решения задачи Коши для дифференциального уравнения первого порядка.

Небольшая подготовка и составление плана

В нашем классе, который назовем NumDifferentiation, должны быть описаны методы, вычисляющие приближенное значение решения задачи Коши с заданной точностью в заданной точке x. Поэтому первым шагом я предлагаю написать класс CauchyConditions, который будет содержать три свойства:

  • функцию первой производной — y' = f(x, y)
  • начальную точку x0
  • точку y0 — значение решения задачи Коши в точке x0

После чего можно приступать к основному классу. Как я сказал ранее, он будет содержать не один метод, а целых четыре, каждый из которых будет находить значения, используя различные численные методы:

  • метод Эйлера
  • метод Эйлера-Коши
  • модифицированный метод Эйлера
  • метод Рунге-Кутты

Подробнее познакомиться с каждым из этих методов Вы можете в этой статье, а здесь я просто выложу формулы.

Приближенное решение задачи Коши на C# (Выч. мат, С#)

Класс CauchyConditions

public class CauchyConditions
{
public double X0 { get; private set; }
public double Y0 { get; private set; }
public Func<double, double, double> FirstDeritative { get; private set; }
public CauchyConditions(double x0, double y0, Func<double, double, double> firstDeritative)
{
this.X0 = x0;
this.Y0 = y0;
this.FirstDeritative = firstDeritative;
}
}

Все три свойства делаем автореализуемыми с доступом "только для чтения" и описываем один конструктор, в котором инициализируем каждое из свойств.

Класс NumDifferentiation

Этот класс получился довольно большим, поэтому будем изучать его по частям, а в конце соберем это все воедино.

Немного о концепции

Хотя публичных методов будет всего четыре, приватных будет намного больше. Все это сделано с целью избавиться от повторяющегося кода. Главный метод — метод под названием RungeMethod. Он вызывается при вычислении каждым численным методом. Его назначение — добиться указанной (переданной в качестве параметра) точности. Этот метод в свою очередь вызывает метод calculation, который собственно и вычисляет нужное значение тем численным методом, который передадим в качестве параметра.

Думаю, так сразу понять довольно сложно, поэтому плавно перейдем к коду.

Публичные методы

// Метод Эйлера
public static double calcByEulerMethod(double x, double eps, CauchyConditions conditions)
{
return RungeMethod(x, eps, conditions, __calcByEulerMethod, 1);
}
// Метод Эйлера-Коши
public static double calcByEulerCauchyMethod(double x, double eps, CauchyConditions conditions)
{
return RungeMethod(x, eps, conditions, __calcByEulerCauchyMethod, 2);
}
// Модифицированный метод Эйлера
public static double calcByModifiedEulerMethod(double x, double eps, CauchyConditions conditions)
{
return RungeMethod(x, eps, conditions, __calcByModifiedEulerMethod, 2);
}
// Метод Рунге-Кутты
public static double calcByRungeKuttaMethod(double x, double eps, CauchyConditions conditions)
{
return RungeMethod(x, eps, conditions, __calcByRungeKuttaMethod, 4);
}

Как и говорилось ранее, в качестве параметров передается точка x, точность и условие Коши, а возвращается результат метода RungeMethod. Вновь вызываемый метод принимает 5 параметров:

  • точку x
  • точность eps
  • условия Коши conditions
  • способ вычисления (функция)
  • порядок точности численного метода (используется для оценки погрешности по принципу Рунге)

Метод RungeMethod

private static double RungeMethod(double x, double eps, CauchyConditions conditions,
Func<double, CauchyConditions, double> wayCalc, int power)
{
double step = (x - conditions.X0) / 2;
double resultFullStep = calculation(x, step, conditions, wayCalc);
double resultHalfStep = calculation(x, step / 2, conditions, wayCalc);
double error = Math.Abs((resultHalfStep - resultFullStep) / (Math.Pow(2, power) - 1));
int i = 0;
// i нужен, чтобы избежать зацикливания,
// когда задана очень маленькая точность
while (error > eps && i < 20)
{
step /= 2;
resultFullStep = calculation(x, step, conditions, wayCalc);
resultHalfStep = calculation(x, step / 2, conditions, wayCalc);
error = Math.Abs(resultHalfStep - resultFullStep);
i++;
}
return resultHalfStep + error;
}

Здесь просто реализована оценка погрешности по принципу Рунге. Значение будет перевычисляться до тех пор, пока заданная точность не будет достигнута, или же количество заходов в цикл не станет равным двадцати.

В метод calculation передаются параметры:

  • та же точка x
  • шаг (используется в формулах численных методов)
  • те же условия Коши
  • способ вычисления (wayCalc)

Метод calculation

private static double calculation(double x, double step, CauchyConditions conditions,
Func<double, CauchyConditions, double> wayCalc)
{
double previousX = conditions.X0;
double result = conditions.Y0;
while (previousX < x)
{
result += wayCalc(step, new CauchyConditions(previousX, result, conditions.FirstDeritative));
previousX += step;
}
return result;
}

wayCalc — любой из нижеперечисленных методов

private static double __calcByEulerMethod(double step, CauchyConditions conditions)
{
return step * conditions.FirstDeritative(conditions.X0, conditions.Y0);
}
private static double __calcByEulerCauchyMethod(double step, CauchyConditions conditions)
{
double x = conditions.X0;
double y = conditions.Y0;
double firstValue = conditions.FirstDeritative(x, y);
double secondValue = conditions.FirstDeritative(x + step, y + step * firstValue);
return step / 2 * (firstValue + secondValue);
}
private static double __calcByModifiedEulerMethod(double step, CauchyConditions conditions)
{
double x = conditions.X0;
double y = conditions.Y0;
double halfStep = step / 2;
double firstValue = conditions.FirstDeritative(x, y);
double secondValue = conditions.FirstDeritative(x + halfStep, y + halfStep * firstValue);
return step * secondValue;
}
private static double __calcByRungeKuttaMethod(double step, CauchyConditions conditions)
{
double x = conditions.X0;
double y = conditions.Y0;
double halfStep = step / 2;
double firstValue = conditions.FirstDeritative(x, y);
double secondValue = conditions.FirstDeritative(x + halfStep, y + halfStep * firstValue);
double thirdValue = conditions.FirstDeritative(x + halfStep, y + halfStep * secondValue);
double fourthValue = conditions.FirstDeritative(x + halfStep, y + step * thirdValue);
return step / 6 * (firstValue + 2 * secondValue + 2 * thirdValue + fourthValue);
}

Эти методы возвращают значение второго слагаемого формулы соответствующего метода, т.е. то, что идет после yi.

#image.jpg

Теперь код всего класса NumDifferentiation

	public static class NumDifferentiation
{
private static double RungeMethod(double x, double eps, CauchyConditions conditions,
Func<double, CauchyConditions, double> wayCalc, int power)
{
double step = (x - conditions.X0) / 2;
double resultFullStep = calculation(x, step, conditions, wayCalc);
double resultHalfStep = calculation(x, step / 2, conditions, wayCalc);
double error = Math.Abs((resultHalfStep - resultFullStep) / (Math.Pow(2, power) - 1));
int i = 0;
while (error > eps && i < 20)
{
step /= 2;
resultFullStep = calculation(x, step, conditions, wayCalc);
resultHalfStep = calculation(x, step / 2, conditions, wayCalc);
error = Math.Abs(resultHalfStep - resultFullStep);
i++;
}
return resultHalfStep + error;
}
private static double calculation(double x, double step, CauchyConditions conditions,
Func<double, CauchyConditions, double> waycalc)
{
double previousX = conditions.X0;
double result = conditions.Y0;
while (previousX < x)
{
result += waycalc(step, new CauchyConditions(previousX, result, conditions.FirstDeritative));
previousX += step;
}
return result;
}
// Формула:
private static double __calcByEulerMethod(double step, CauchyConditions conditions)
{
return step * conditions.FirstDeritative(conditions.X0, conditions.Y0);
}
private static double __calcByEulerCauchyMethod(double step, CauchyConditions conditions)
{
double x = conditions.X0;
double y = conditions.Y0;
double firstValue = conditions.FirstDeritative(x, y);
double secondValue = conditions.FirstDeritative(x + step, y + step * firstValue);
return step / 2 * (firstValue + secondValue);
}
private static double __calcByModifiedEulerMethod(double step, CauchyConditions conditions)
{
double x = conditions.X0;
double y = conditions.Y0;
double halfStep = step / 2;
double firstValue = conditions.FirstDeritative(x, y);
double secondValue = conditions.FirstDeritative(x + halfStep, y + halfStep * firstValue);
return step * secondValue;
}
private static double __calcByRungeKuttaMethod(double step, CauchyConditions conditions)
{
double x = conditions.X0;
double y = conditions.Y0;
double halfStep = step / 2;
double firstValue = conditions.FirstDeritative(x, y);
double secondValue = conditions.FirstDeritative(x + halfStep, y + halfStep * firstValue);
double thirdValue = conditions.FirstDeritative(x + halfStep, y + halfStep * secondValue);
double fourthValue = conditions.FirstDeritative(x + halfStep, y + step * thirdValue);
return step / 6 * (firstValue + 2 * secondValue + 2 * thirdValue + fourthValue);
}
// Метод Эйлера
public static double calcByEulerMethod(double x, double eps, CauchyConditions conditions)
{
return RungeMethod(x, eps, conditions, __calcByEulerMethod, 1);
}
// Метод Эйлера-Коши
public static double calcByEulerCauchyMethod(double x, double eps, CauchyConditions conditions)
{
return RungeMethod(x, eps, conditions, __calcByEulerCauchyMethod, 2);
}
// Модифицированный метод Эйлера
public static double calcByModifiedEulerMethod(double x, double eps, CauchyConditions conditions)
{
return RungeMethod(x, eps, conditions, __calcByModifiedEulerMethod, 2);
}
// Метод Рунге-Кутты
public static double calcByRungeKuttaMethod(double x, double eps, CauchyConditions conditions)
{
return RungeMethod(x, eps, conditions, __calcByRungeKuttaMethod, 4);
}
}

2 отзыва на “Приближенное решение задачи Коши на C#

  1. Сергей

    Спасибо! Очень полезная статья! Вы не могли бы описать процедуру решения модифицированным методом Эйлера для ОДУ 2-го порядка?

  2. Mike_Device

    Благодарю за отзыв! Возьму на заметку предложенную тему, однако сейчас не так много времени, поэтому не могу сказать, когда появится подобная статья.

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *