Целью данной работы, есть обучение студентов применению пакета MATHCAD для ответа дифференциальных уравнений в функции одной настоящей переменной на ПЭВМ.
Введение.
Численные способы ответа задачи Коши.Задача ответа обычных дифференциальных уравнений сложнее задачи вычисления однократных интегралов, а
часть задач, интегрируемых в явном виде, тут значительно
меньше.
В рассуждениях об интегрируемости в явном виде в большинстве случаев имеется в
виду, что ответ возможно вычислено при помощи конечного числа «элементарных» операций: сложения, вычитания, умножения, деления, возведения
в степень, логарифмирования, потенцирования, косинуса и вычисления синуса. Еще во время, предшествующий появлению ЭВМ, понятие «элементарной» функции и «элементарной» операции претерпело изменение. Решения некоторых частных задач так довольно часто виделись в приложениях, что было нужно составить таблицы их значений, к примеру, таблицы интегралов Френеля
интегрального синуса ответов дифференциального уравнения вида
х2у + ху’ + (х2 — ?2)у = 0,
именуемых функциями Бесселя, и последовательности вторых функций, именуемых особыми функциями. При наличии таких таблиц исчезает принципиальная отличие между вычислением функций sin(x), ln(х), и т.д. и особых
функций. В том и втором случаях возможно вычислять значения этих функций
при помощи таблиц, и те и другие функции возможно вычислять, приближая
многочленами, рациональными дробями и т.д. Так, в класс задач,
интегрируемых в явном виде, включились задачи, решения которых выражаются через особые функции. Но и данный, более широкий, класс
образовывает довольно малую долю задач, предъявляемых к ответу. Значительное расширение класса реально решаемых дифференциальных уравнений, а следовательно, и расширение сферы применения математики случилось с созданием численных способов и распространением применения
ЭВМ.
На данный момент объем работы по подготовке ответа на ЭВМ
дифференциального уравнения не так уж значительно превышает объем работы по написанию этого уравнения; при жажде возможно конкретно
измашины взять график ответа либо изображение ответа на экране.
В следствии этого отпала необходимость в изучении интегрирования уравнений и большого числа в явном виде.
Разложение ответа в ряд Тейлора.
Одним из старейших способов ответа дифференциальных
уравнений есть способ разложения в ряд Тейлора.
Пускай требуется обнаружить отрезке [x0, x0 + Х] ответ дифференциального уравнения
у’ = f(x,у) (7-1)
при начальном условии у(х0) = у0, f(х, у) аналитична в точке
(х0, y0). Дифференцируя (1.1) по х, возьмём соотношения:
у=fx(x,y) + fy(х,y)*y’,
у»’ = fxx(х,у) + 2fxy(х,у) у’+ fyy(х, у) у’2 + fy(х, у) у», …
Подставляя в это выражение х=x0 и у = у0, приобретаем последовательно значения
у’ (x0), у (x0), у'(х0), …
Так, возможно написать приближенное равенство
(7-2)
В случае, если |х — x0| больше радиуса сходимости последовательности
,
то погрешность выражения (1.2) не пытается к нулю при n ?. Исходя из этого время от времени целесообразно поступить следующим образом.
Разобивают промежуток [х0, х0 + Х] на отрезки
[xj-1,xj], j = 1, …, N;
и последовательно приобретаем приближения уj к значениям ответа
у(хj), j = l, …, N,
по следующему правилу: пускай значение уj уже отыскано, вычисляем значения в точке хj производных уj(i) решения исходного дифференциального уравнения, проходящего через точку (хj, уj);
на отрезке [хj, хj+1] полагаем
(7-3)
и полагаем
уj+1 = zj (хj+1) (7-4)
Потом мы совершим подробную оценку погрешности одно-
шаговых способов интегрирования, к каким относится разглядываемый метод. На настоящем этапе дискуссии возможно
привести следующие мысли в пользу того, что погрешность этого способа может оказаться малой. Разглядим слу-
чай хj+1 — хj = h; если бы значение уj совпадало со значением
правильного значения у(хj), то погрешность от замены уj+1 на
zj(хj+1) имела бы порядок O(hn+1). Потому, что мы вносим по-
грешность на O(h-1) отрезках, то возможно ожидать, что при измельчении шага сетки будет выполняться соотношение
Это рассуждение есть обычным примером правдоподобных рассуждений, приводящихся время от времени в обоснование приме-
нения тех либо иных численных способов; но направляться иметь
в виду, что во многих случаях видится такая обстановка: в
обоснование сходимости способа такие рассуждения возможно про-
вести, а в действительности сходимости приближенного ответа к
правильному нет. Исходя из этого вопрос строгого обоснования сходимости
способов при “мельчении” шага, и получения оценки погрешности есть вопросом, имеющим не только теоретическое,
но и наиболее значимое практическое значение.
Обрисованный выше метод интегрирования дифференциальных уравнений
при n1 был исключен из широкой практики более столетия назад.
При применении этого способа необходимо вычислять значения функции f и всех ее производных fxjym-j при m Кроме других обстоятельств, при ручном счете, это не очень приятно еще и
тем, что рассеивается внимание расчетчика и возрастает возможность неточностей.
При работе с ЭВМ использование этого способа требует написания солидного
числа блоков вычисления производных, что противоречит главной тенденции
упрощения взаимоотношений между машиной и пользователем.
На данный момент имеются совокупности программ, каковые реализовывают
формульное дифференцирование функций; при заданной программе вычисления функции f(х,у) эта совокупность программ расписывает программы вычисления производных функции f(x,у); так, отпадает необходимость
написания программ вычисления этих производных. Но по затратам машинного времени подобные методы, в большинстве случаев, значительно усту-
пают разглядываемым ниже способам типа Рунге — Кутта и конечно-разностным способам.
Отечественные рассуждения не нужно осознавать, как окончательные советы воздерживаться от применения способов, применяющих разложение ответа в ряд Тейлора. Во многих случаях ЭВМ употребляются для ответа дифференциальных уравнений строго фиксированного типа. К примеру, при расчетах траекторий перемещения небесных тел приходится многократно интегрировать совокупности дифференциальных уравнений в полной мере определенного вида при
разных различных значениях и начальных условиях параметров в правых
частях. То событие, что все время решается одинаковая совокупность дифференциальных уравнений, может дать следующее преимущество: конкретные формулы для производных правых частей совокупности в большинстве случаев имеют довольно много
неспециализированных частей; в следствии для отдельных классов совокупностей может произойти,
что одновременное вычисление всех этих производных потребует, относительно малого числа арифметических операций и разглядываемый способ
окажется действеннее вторых способов численного интегрирования. В будущем мы, по большей части, будем изучать способы численного интегрирования,
требующие вычисления лишь значений правой части; на базе этих способов
в большинстве случаев и строятся стандартные программы ответа широких классов типовых задач.
Способы Рунге — Кутта
В частном случае n = 1 расчетная формула (1.3) принимает вид
(7-5)
Данный способ именуется способом Эйлера. Возможно выстроить второй класс расчетных формул, к у которого в собствености способ
Эйлера. Укажем сперва несложные способы этого класса, приобретаемые из наглядных мыслей. Пускай известно значение
ответа у(х) и требуется вычислить значение у(х+h). Рас-
наблюдаем равенство
(7-6)
При замене интеграла в правой части на величину hy'(х) по-
грешность имеет порядок O(h2),
т. е.
потому, что у'(х) = f(х, у(х)), то из этого имеем
Отбрасывая член порядка O(h2) и обозначая х = хj, х+h =
хj+1, возьмём расчетную формулу Эйлера (7-5). Для получения более правильной расчетной формулы необходимо правильнее аппроксимировать интеграл в правой части (7-6). Воспользовавшись формулой трапеций, возьмём
7-7
либо
в противном случае,
Соответствующая расчетная формула:
В большинстве случаев это уравнение неразрешимо очевидно довольно уj+1.
Исходя из этого произведем предстоящее преобразование метода.
Заменим у(х+h) в правой части (3) на некую величину
тогда правая часть изменится на величину
где у находится между у* и у(х+h); благодаря предположе-
ния (5) эта величина имеет порядок O(h3). Так, при
условии (5)
условию (5) удовлетворяет итог вычислений, по расчет-
ной формуле Эйлера
у* = y(x) + hf(x, у (х)).
Последние соотношения определяют несколько расчетных формул
Выстроим другую несколько формул с погрешностью на шаге та-
кого же порядка. Интеграл в правой части (7-6) заменим по фор-
муле прямоугольников:
y(x+h) = y(x) + hy’(x +h/2) + O(h3)
либо, все равно, что
в случае, если
то, как и в предшествующем случае, имеем
у(х+ h) = у(х)+ hf (x+h/2, у*)+ 0 (h3).
В качестве у* возможно забрать итог вычислений по формуле
Эйлера с шагом h/2:
у*= у(х)+h/2*f(х, у(х)).
Этим соотношениям соответствует пара расчетных формул
уi+1/2 = уi+ h/2 f ( xi,yj) (7-8)
yi+1 = уi+f (хi+ h/2, уj+1/2 )
Полученные способы относятся к семейству способов Рунге-
-Кутта, имеющих следующий вид. В ходе вычислений фиксированы кое-какие числа
а2, …, аq, p1, …, Pq (if 0
последовательно приобретаем
k1(h) = hf(x,у),
k2(h) = hf(х+a2h), у+~21k1(h)),
kq(h) = hf(х+ аqh, у+ pq,1k1(h)+ … + p
и полагаем
y(x+h) ? z(h) = y(x) +
Разглядим вопрос о выборе параметров ?i, рi, ?i,j. Обозна-
чим ?(h) = у(х+h) — z(h). В случае, если f(x, у) — достаточно ровная функция собственных доводов, то k1(h), …, kq(h) и ?(h) —
ровные функции параметра h. Предположим, что
?(0)= ?’(0))= … = ?(s)(h) =0
при произвольных достаточно ровных функциях f(x,у), а
?(s+1)(0) ? 0 для некоей ровной функции f(x,у). В соответствии с формуле
Тейлора, справедливо равенство
(7-9)
где 0 ?
Порядок исполнения работы.
I. Делим отрезок [а, b] на n частей:
a = x0 x1 x2 … xi xi+1 … xn = 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 выяснен без его перечня параметров на решающем устройстве ОДУ, и функции yn(x) выяснены в D без их параметра. Второй параметр D должен быть задан сценарием как вектор, даже в том случае, если имеется лишь одиночный элемент.
Rkadapt — применяет неравномерные размеры внутреннего шага, в то время, когда решается дифференциально-разностное уравнение, добавляется большее количество шагов в области большего трансформации ответа, но возвращает ответ в равномерно распределенных точках.
Исполнение работы.
Ниже приведен пример ответа задачи охлаждения нагретого тела до температуры воздуха способом Рунге-Кута. Ответ ищется, в виде численного уравнения вида:
yn+1= yn+x/6(k1(yn)+1k2(yn)+2k3(yn)+k4(yn))
Температура тела изменяется от начальной температуры заданной в переменной y0 изменяющейся по экспоненциальному закону вида
r(at+ts).
Выполняется блок — схема метода на отдельном странице и строится график, как пример, рассмотрено ответ дифференциального уравнения df(t,y) = -r(y0 – Ts), прилагается распечатка ответа данной задачи с таблицей данных и графиком. В пакете Mathcad соответствующие блоки будут выглядеть следующим образом:
1. Вводим начальные условия:
Ts = 220С — температура воздуха;
y0 = 900С — температура нагретого тела;
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
Назад
ЛАБОРАТ0РНАЯ РАБОТА № 8