Численное решение дифференциальных уравнений

Численное решение дифференциальных уравнений

Методы решения дифференциальных уравнений очень сложны и во многом определяются характером этих уравнений. Во многих случаях подготовка программы для решения этих уравнений проводится строго индивидуально. Однако, есть некоторые достаточно общие методы решения таких уравнений и, что самое главное, есть определенная проверенная практикой методика получения результатов решения в виде графика.

Мы рассмотрим решение двумя методами:

·  методом Эйлера;

·  методом Рунге-Кутта.

Эти методы имеют наиболее широкое распространение на практике

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

Задание состоит в следующем.

1.  Подготовить алгоритм решения заданного уравнения или системы уравнений указанным в задании методом.

2.  Подготовить программу для алгоритма п. 1 в виде полностью готового программного продукта.

3.  Провести отладку программы и убедиться, что она отражает выбранный алгоритм (например, путем решения некоторого уравнения, имеющего точное решение, и последующей логической проверки полученного решения).

4.  Выполнить решение заданного уравнения при указанных в задании условиях.

5.  В процессе выполнения каждое приближение выводить на экран и определить сколько итераций необходимо для получения нужной точности.

Рассмотрим методы решения дифференциальных уравнений.

1. Метод Эйлера

Дифференциальное уравнение
y¢ = f(x, y) (1)
определяет на плоскости поле направлений, т.е. определяет в каждой точке плоскости, в которой существует функция f(x, y), направление интегральной кривой уравнения, проходящей через эту точку.

Пусть необходимо решить задачу Коши для (1), удовлетворяющую начальному условию y(x₀) = y₀.

Разделим отрезок [x₀, X] на n равных частей (X – x₀)/n = h, где h – шаг изменения аргумента. Тогда значения решения в следующей после начальной точке будет y₁ » y₀ + h × f(x₀, y₀), если участок достаточно мал для того, чтобы считать на нем функцию постоянной.

Общее решение примет вид: y_k+1 » y_k + h × f(x_k, y_k).

Для определения сходимости шаг решения делится пополам h¹ = h/2 и все решение повторяется. Для каждой точки, имеющей координату x_k = x₀ + k × h (т.е. в первоначально заданных точках поиска решения) определяется разность
a_k = abs(y_k^i-1 – y_kⁱ), где i – номер итерации. Решение прекращается, если max(a_k) < e, где e – точность вычисления.

Аналогично для системы уравнений

y₁¢ = f₁(x, y₁, ..., y_n)

...

y_n¢ = f_n(x, y₁, ..., y_n)

получим

y₁(x_k + h) » y₁(x_k) + h × f₁(y₁(x_k), ..., y_n(x_k))

...

y_n(x_k + h) » y_n(x_k) + h × f₁(y₁(x_k), ..., y_n(x_k))

Далее все аналогично решению одного уравнения.

2. Метод Рунге-Кутта

Рассмотрим уравнение y¢ = f(x, y) с начальным условием y(x₀) = y₀.

Для вычисления y в каждой следующей после начальной точке определяются 4 числа:

k₁ = h × f(x, y)

k₂ = h × f(x + h/2, y + k₁/2)

k₃ = h × f(x + h/2, y + k₂/2)

k₄ = h × f(x + h, y + k₃)

Тогда приращение y определяется по формуле: Dy = (k₁ + 2 × k₂ + 2 × k₃ + k₄)/6.

Расчет повторяется во всех узловых точках, а узловые точки определяются так же, как и для метода Эйлера.

Для систем уравнений все выполняется аналогично, но количество вычислений возрастает. Например, для системы 2-х уравнений

y¢ = f(x, y, z)

z¢ = g(x, y, z)

с заданными начальными приближениями x₀, y₀, z₀ значения y₁ и z₁ в точке x₁ определяются следующим образом:

k₁ = h ×f(x₀, y₀, z₀)                                           l₁ = h × g(x₀, y₀, z₀)

k₂ = h ×f(x₀ + h/2, y₀ + k₁/2, z₀ + l₁/2)                        l₂ = h × g(x₀ + h/2, y₀ + k₁/2, z₀ + l₁/2)

k₃ = h ×f(x₀ + h/2, y₀ + k₂/2, z₀ + l₂/2)                        l₃ = h × g(x₀ + h/2, y₀ + k₂/2, z₀ + l₂/2)

k₄ = h ×f(x₀ + h, y₀ + k₃, z₀ + l₃)                      l₄ = h × g(x₀ + h, y₀ + k₃, z₀ + l₃)

Dy = (k₁ + 2 × k₂ + 2 × k₃ + k₄)/6                     Dz = (l₁ + 2 × l₂ + 2 × l₃ + l₄)/6

Задание

№	Уравнение	Отрезок	Нач. приближ.	Метод
1		0 £ x £ 1 hнач = 0.1	v(0) = 1 w(0) = 2	Метод Эйлера
2		0 £ x £ 1 hнач = 0.1	y(0) = 0 y¢(0) = 1	Метод Эйлера
3		0 £ t £ 1 hнач = 0.1 a = 0.1; 1.0; 10.0	x(0) = 0 y(0) = 1	Метод Эйлера
4		0 £ x £ 1 hнач = 0.1	y(0) = 0	Метод Эйлера
5		0 £ x £ 1 hнач = 0.1	y(0) = 1	Метод Эйлера
6		1 £ t £ 2 hнач = 0.1	x(1) = 1 y(1) = 1	Метод Эйлера
7		0 £ x £ 1 hнач = 0.1	y(0) = 1	Метод Эйлера
8		2 £ x £ 6 hнач = 0.1	y(2) = 4	Метод Эйлера
9		0 £ x £ 2 hнач = 0.1	y(0) = 0	Метод Эйлера
10		1 £ x £ 2 hнач = 0.1	y(1) = 0 Точное решение	Метод Эйлера
11		0 £ x £ 1 hнач = 0.1	y(0) = 1	Метод Эйлера
12		0 £ x £ 1 hнач = 0.1	y(0) = 1	Метод Эйлера

Варианты с 13 по 24 соответствуют вариантам 1-12, но вместо метода Эйлера следует использовать метод Рунге-Кутта.

На практической работе особое внимание следует уделить процессу получения решения последовательными итерациями в каждом из предложенных методов. Определите количество итераций для получения решения каждым из указанных методов.