Решение дифференциальных уравнений в функции одной переменной на заданном отрезке методом Рунге-Кута в среде MATHCAD

Целью данной работы, является обучение студентов применению пакета MATHCAD для решения дифференциальных уравнений в функции одной действительной переменной на ПЭВМ.

Введение.

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

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

интегрального синуса решений дифференциального уравнения вида

х²у" + ху' + (х² — ν ²)у = 0,

называемых функциями Бесселя, и ряда других функций, называемых специальными функциями. При наличии таких таблиц исчезает принципиальная разница между вычислением функций sin(x), ln(х), и т.д. и специальных
функций. В том и другом случаях можно вычислять значения этих функций
при помощи таблиц, и те и другие функции можно вычислять, приближая
многочленами, рациональными дробями и т.д. Таким образом, в класс задач,
интегрируемых в явном виде, включились задачи, решения которых выражаются через специальные функции. Однако и этот, более широкий, класс
составляет относительно малую долю задач, предъявляемых к решению. Существенное расширение класса реально решаемых дифференциальных уравнений, а следовательно, и расширение сферы применения математики произошло с созданием численных методов и распространением применения
ЭВМ.

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

Разложение решения в ряд Тейлора.

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

Пусть требуется найти на отрезке [x₀, x₀ + Х] решение дифференциального уравнения

у' = f(x, у) (7-1)

 

при начальном условии у(х₀) = у₀, f(х, у) аналитична в точке
(х₀, y₀). Дифференцируя (1.1) по х, получим соотношения:

у " =f_x(x, y) + f_y(х, y)*y',

у ^''' = f_xx(х, у) + 2f_xy(х, у) у ^' + f_yy(х, у) у^'2 + f_y(х, у) у ^'', ...

Подставляя в это выражение х = x₀ и у = у₀, получаем последовательно значения

у' (x₀), у" (x₀), у" '(х₀), ...

Таким образом, можно написать приближенное равенство

 

(7-2)

Если |х — x₀| больше радиуса сходимости ряда

,

то погрешность выражения (1.2) не стремится кнулю при n → ∞ . Поэтому иногда целесообразно поступить следующим образом.

Разобивают интервал [х₀, х₀+ Х] на отрезки

[x_j-1, x_j], j = 1, ..., N;

и последовательно получаем приближения у_j кзначениям решения

у(х_j), j = l, ..., N,

по следующему правилу: пусть значение у_j уже найдено, вычисляем значения в точке х_jпроизводных у_j⁽ⁱ⁾ решения исходного дифференциального уравнения, проходящего через точку (х_j, у_j);
на отрезке [х_j, х_j+1] полагаем

 

(7-3)

 

 

и, соответственно, полагаем

у_j+1 = z_j (х_j+1) (7-4)

Далее мы проведем подробную оценку погрешности одно-
шаговых способов интегрирования, к которым относится рассматриваемый способ. На настоящем этапе обсуждения можно
привести следующие соображения в пользу того, что погрешность этого метода может оказаться малой. Рассмотрим слу-
чай х_j+1 - х_j = h; если бы значение у_j совпадало со значением
точного значения у (х_j), то погрешность от замены у_j+1 на
z_j(х_j+1) имела бы порядок O(hⁿ⁺¹). Поскольку мы вносим по-
грешность на O(h^-1) отрезках, то можно ожидать, что при измельчении шага сетки будет выполняться соотношение

 

 

 

Это рассуждение является типичным примером правдоподобных рассуждений, приводящихся иногда в обоснование приме-
нения тех или иных численных методов; однако следует иметь
в виду, что во многих случаях встречается такая ситуация: в
обоснование сходимости метода такие рассуждения можно про-
вести, а на самом деле сходимости приближенного решения к
точному нет. Поэтому вопрос строгого обоснования сходимости

методов при “мельчении” шага, а также получения оценки погрешности является вопросом, имеющим не только теоретическое,
но и важнейшее практическое значение.

Описанный выше способ интегрирования дифференциальных уравнений
в случае n> 1 был исключен из широкой практики более столетия назад.
При использовании этого метода нужно вычислять значения функции f и всех ее производных f_x^j_y^m-j при m< n т.е. вычислить n(n+1)/2 различных функций.
Кроме других причин, при ручном счете, это неприятно еще и
тем, что рассеивается внимание расчетчика и возрастает вероятность ошибок.
При работе с ЭВМ применение этого метода требует написания большого
числа блоков вычисления производных, что противоречит основной тенденции
упрощения отношений между пользователем и машиной.

В настоящее время имеются системы программ, которые осуществляют
формульное дифференцирование функций; при заданной программе вычисления функции f(х, у) эта система программ расписывает программы вычисления производных функции f(x, у); таким образом, отпадает необходимость
написания программ вычисления этих производных. Однако по затратам машинного времени подобные алгоритмы, как правило, существенно усту-
пают рассматриваемым ниже методам типа Рунге — Кутта и конечно-разностным методам.

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

Методы Рунге — Кутта

В частном случае n = 1 расчетная формула (1.3) принимает вид

 

(7-5)

 

 

Этот метод называется методом Эйлера. Можно построить другой класс расчетных формул, ккоторому принадлежит метод
Эйлера. Укажем сначала простейшие методы этого класса, получаемые из наглядных соображений. Пусть известно значение
решения у(х) и требуется вычислить значение у(х+h). Рас-
смотрим равенство

 

(7-6)

 

При замене интеграла в правой части на величину hy'(х) по-
грешность имеет порядок O(h²),

т. е.

 

 

поскольку у'(х) = f(х, у(х)), то отсюда имеем

 

Отбрасывая член порядка O(h²) и обозначая х = х_j, х+h =
х_j+1, получим расчетную формулу Эйлера (7-5). Для получения более точной расчетной формулы нужно точнее аппроксимировать интеграл в правой части (7-6). Воспользовавшись формулой трапеций, получим

7-7

 

или
иначе,

 

Соответствующая расчетная формула:

 

 

Обычно это уравнение неразрешимо явно относительно у_j+1.
Поэтому произведем дальнейшее преобразование алгоритма.
Заменим у(х+h) в правой части (3) на некоторую величину

 

 

тогда правая часть изменится на величину

 

где у находится между у^* и у(х+h); вследствие предположе-
ния (5) эта величина имеет порядок O(h³). Таким образом, при
условии (5)

 

 

условию (5) удовлетворяет результат вычислений, по расчет-
ной формуле Эйлера

у^* = y(x) + hf(x, у (х)).

Последние соотношения определяют пару расчетных формул

 

 

Построим другую пару формул с погрешностью на шаге та-
кого же порядка. Интеграл в правой части (7-6) заменим по фор-
муле прямоугольников:

y(x+h) = y(x) + hy ^’ (x +h/2) + O(h³)

или, все равно, что
если

то, как и в предшествующем случае, имеем

у(х+ h) = у(х)+ hf (x+h/2, у^*)+ 0 (h³).

 

В качестве у* можно взять результат вычислений по формуле
Эйлера с шагом h/2:

у ^* = у(х)+h/2*f(х, у(х)).

Этим соотношениям соответствует пара расчетных формул

 

у_i+1/2 = у_i+ h/2 f ( x_i, y_j) (7-8)

 

y_i+1 = у_i+f (х_i+ h/2, у_j+1/2 )

Полученные методы относятся к семейству методов Рунге-
-Кутта, имеющих следующий вид. В процессе вычислений фиксированы некоторые числа

а₂, ..., а_q, p₁, ..., P_q (_if 0< j< i< q;

последовательно получаем

 

k₁(h) = hf(x, у),

 

k₂(h) = hf(х+a₂h), у+~₂₁k₁(h)),

k_q(h) = hf(х+ а_qh, у+ p_{q, 1}k₁(h)+... + p<, /г, (h))

и полагаем

y(x+h) ≈ z(h) = y(x) +

Рассмотрим вопрос овыборе параметров α _i, р_i, β _{i, j}. Обозна-
чим φ (h) = у(х+h) - z(h). Если f(x, у) - достаточно гладкая функция своих аргументов, то k₁(h), ..., k_q(h) и φ (h) -
гладкие функции параметра h. Предположим, что

φ (0)= φ ’(0))=... = φ ^(s)(h) =0

 

при произвольных достаточно гладких функциях f(x, у), а
φ ^(s+1)(0) ≠ 0 для некоторой гладкой функции f(x, у). Согласно формуле
Тейлора, справедливо равенство

 

(7-9)

 

где 0< Θ < 1. Величина φ (h) называется погрешностью метода на шаге, а s порядком погрешности метода.

Порядок выполнения работы.

I. Делим отрезок [а, b] на n частей:

a = x₀ < x₁ < x₂ < ... < x_i < x_i+1 < ... < x_n = b

Находим шаг вычисления

h = (а - b)/n

В Mathcad имеются следующие встроенные функции:

- Rkfixed (init, x1, x2, npoints, D) - используется для фиксированного шага метода Рунге-Кутта.

- Rkadapt (init, x1, x2, npoints, D) - используется для метода Рунге-Кутта с адаптивным размером шага.

- Bulstoer (init, x1, x2, npoints, D) - Bulirsch-Stoer метод более точный, но применяется только для гладко изменяющихся функций.

Каждая из этих функций возвращает матрицу, в которой первый столбец содержит “n points” – n точек между x1 и x2, отрезка на котором ищется решение Обыкновенного Дифференциального Уравнения (ОДУ). Последующие столбцы содержат соответствующие значения матрицы ОДУ с n-ым порядком.

В работе будем использовать встроенную Mathcad функцию - rkfixed.

Объяснения:

· Init – является, или вектором n реальных начальных значений, или одиночного скалярного начального значения, в случае одиночной ОДУ.

· х1 и x2 - реальные, скалярные концевые точки интервала, по которому ищется решение ОДУ. Начальные значения в init - значения функции ОДУ.

· N-points - целое число точек вне исходной точки, в которых ищется решение, при необходимости они должны быть аппроксимированы. Оно определяет число строк (n-points +1) в матрице решений.

· D(x, y) - n-элементный вектор решения функции независимой переменной x, и вектора функций y, содержит уравнения для первых производных всех неизвестных функций в системе ОДУ. Чтобы создать этот вектор, уравнение с первыми производными приводиться отдельно на левом - hand-side, без множителей, и никаких более высокого порядка производных в уравнении. Например, одиночная ОД функция y’(x), которая содержит вторую производную, должна быть написана как система уравнений в y0(x) и y1(x), где первая производная y0 - y1. Одна обыкновенная дифференциальная функция

перезаписывается для решающего устройства, используя векторные матрицы, как

Примечания:

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

Rkadapt - использует неравномерные размеры внутреннего шага, когда решается дифференциально-разностное уравнение, добавляется большее количество шагов в области большего изменения решения, но возвращает решение в равномерно распределенных точках.

Выполнение работы.

Ниже приведен пример решения задачи охлаждения нагретого тела до температуры окружающей среды методом Рунге-Кута. Решение ищется, в виде численного уравнения вида:

y_n+1= y_n+x/6(k₁(y_n)+1k₂(yn)+2k3(yn)+k4(yn))

Температура тела изменяется от начальной температуры заданной в переменной y₀ изменяющейся по экспоненциальному закону вида

r(at+t_s).

Выполняется блок - схема алгоритма на отдельном листе и строится график, как пример, рассмотрено решение дифференциального уравнения df(t, y) = -r(y₀ – T_s), прилагается распечатка решения этой задачи с таблицей данных и графиком. В пакете Mathcad соответствующие блоки будут выглядеть следующим образом:

1. Вводим начальные условия:

T_s = 22⁰С - температура окружающей среды;

y₀ = 90⁰С - температура нагретого тела;

r – коэффициент линеаризации;

Записываем начальные условия:

Отрезок 90 ÷ 22 делится на n =1000 частей

Решение выводится в массив Z. Все значения индексируются от 0 до

n-1 с шагом Δ x

Вид изменения температуры представлен на графике

Содержание отчета

Отчет к лабораторной работе должен включать следующие разделы:

· математическая постановка задачи;

· алгоритм задачи;

· блок-схему алгоритма;

· программа, составленная по блок-схеме;

· результат вычисления задачи по программе на ЭВМ;

· для задач, допускающих получение точного решения аналитическим методом, привести это решение и сравнить с полученным результатом. Привести абсолютную и относительную ошибки!

Контрольные вопросы.

1. Как решается задача о нахождении наибольшего (наименьшего) значения функции одной независимой переменной аналитически?

2. Почему не всегда можно решить задачу аналитически (учитывая использование численных методов для приближенного решения трансцендентных уравнений y^’(x) = 0)?

3. Что такое функция rkfixed и как она работает?

4. В чем достоинства и недостатки метода Рунге-Кута?

5. Как будет выглядеть программа реализации поиска минимума?

6. Если функция определена, но имеет разрыв в конечном числе точек отрезка [a, b], можно ли применить метод Рунге-Кута?

 

Варианты заданий.

Найти решение функции y = f (x) на отрезке [ a, b ] методом Рунге-Кута, в соответствии с вариантом заданным преподавателем. Сравнить результат с точным значением (решением методом Ньютона-Лейбница), оценить абсолютную и относительную погрешности вычислений.

 

Таблица вариантов 7-1

 

№ варианта	Функция y = f (x)	Отр-к инт-ния [a, b]	Точное значение
m = min f(x)[a, b]	M = max f(x)[a, b]
	y = (x-3)2*e|x|	[-1, 4]
	y = (x-3)3*e|x+1|	[-2, 4]
	y = e	[-2, 1]
	y=ln(1+ )	[-2, 1]
	y=arcctg	[-1, 2]
	y=	[-1, 2]
	y = arcctg	[-1, 2]
	y = cos(x + π /3signx) + sin(x+π /6)	[-π, π ]
	y=sin(x-π /4) - cos(x -2π /3*signx)	[-π, π ]
	y = sin(x-π /3) - cos(x-2π /3*signx)	[-π, π ]
	y = 2arcctgx + arcsin(2x/(x2 + 1)	[-π, π ]
	y = 4x + (9π 2)/x + sinx	[π, 2π ]
	y = cos2x + cos2(π /3 + x)-cosx*cos(π /3+x)	[- π, π ]
	y =2Sinx + Sin2x	[0, 3/2π ]
	y = x - 2lnx	[3/2, e]
	y =	[-2, 1]
	y=	[-1, 1]
	y =	[0, 1]
	y = x5 - 5x4 + 5x3 + 1	[-1, 2]
	y = x - 2	[0, 5]
	y = 0, 625x4 + 0, 5x2 - 1, 2502x	[0, 3]
	y = 0, 5x2 + x + -2, 5x	[0, 3]
	y = 0, 0393x3 - x2lnx + 0, 25x2	[1, 2]
	y = 1, 5x2 - 4xlnx + 4x	[2, 4]
	y = ex + e-x - 2x	[0, 1]
	y = 1, 5x2 - 14x + ex + e-x	[1, 3]
	y = 2x - 0, 5x2 – Cosx -(1+x)ln(1+x)	[0, 3]
	y = -	[-2, 1]
	Y = xlnx – ½ *x2 + 0, 8x	[-0, 5, 3]
	y = 1/3*x3 - (1+x)ln(1+x) + x -2	[-0, 5, 3]
	y = x/lnx	[-1, 1]
	y = 2x2 - lnx	[-1, e]
	y = x - 2sinx	[0, 2π ]
	y = 2sinc + cos2x	[0, 2 π ]
	y = x	[0, 2]
	y = 2x3 - 3x3	[-1, 1]
	y = x4 - 2x2 + 5	[-2, 2]
	y = x + 2	[0, 4]
	y = x5 - 5x2 + 5x3 + 1	[-1, 2]
	y = x3 - 3x2 + 6x - 2	[-1, 1]
	y =	[-6, 8]
	y = (1 – x + x2)/(1 + x - x2)	[0, 1]
	y = (x - 1)/(x + 1)	[0, 4]
	y = 1/x + 1/(1 - x)	[0, 1]
	y = sin2x - x	[-π /2, π /2]
	y = 2tgx – tg2x	[0, π /2]
	y = x4	[0, 1, +∞ ]
	y =	[0, 3]
	y = arctg	[0, 1]
	y = x*sinx + cosx – (¼ )*x2	[-π /2, π /2]

Назад

 

 

ЛАБОРАТ0РНАЯ РАБОТА № 8

⇐ Предыдущая234567891011Следующая ⇒

Поделиться: