Тема: «Вычисление определенных интегралов»
Цель работы: изучение методов вычисления определенных интегралов, оценка точности этих методов. Получение практических навыков программной реализации методов вычисления определенных интегралов.
Технические средства: IBM PS AT.
Программное средство: Visual Studio или Delphi
Теоретические сведения
При решении многих физических и геометрических задач приходится искать неизвестную функцию по данному соотношению между неизвестной функцией, ее производными и независимыми переменными. Такое соотношение называется дифференциальным уравнением, а отыскание функции, удовлетворяющей уравнению, называется решением или интегрированием данного уравнения.
Простейшее дифференциальное уравнение первого порядка у' = f(у, t) имеет семейство решений у = у(t), являющихся интегральными кривыми второго порядка, чаще всего вида у = Сеt, с произвольной постоянной С. Тогда выбор начального значения у(0)=у0 определяет одну из кривых семейства решений, которая и будет считаться решением поставленной задачи.
Основной считается задача Коши в разделе высшей математики "решение дифференциальных уравнений приближенными методами", формулирующаяся традиционно так: требуется найти решение уравнения у' = f(у, t) в виде у = у(t), удовлетворяющее начальному условию у(0) = у0. Иными словами, требуется найти среди семейства интегральных кривых, являющихся частными решениями дифференциального уравнения у' = f(у, t), интегральную кривую у(t), которая проходит через заданную точку М(t0, у0).
|
|
Для двумерного случая задача Коши формулируется так: если f(х,у) непрерывна в некоторой области R, определенной как |х - х0| < а; |у - у0| < b, то существует, по крайней мере, одно решение у = у (х, t), определенное в окрестности |х - х0| < h, где h > 0, и это решение будет единственным, если выполняется условие Липшица
|f(х,y) - f(х,у)| < N . |y - у|,
где N - некоторая константа Липшица, зависящая в общем случае от а и b.
Заметим, что решение системы дифференциальных уравнений первого порядка
при условии, что производные df/dу, df/dх, dg/dу и dg/dх существуют на каждом интервале интегрирования, содержит уже две постоянные интегрирования и, следовательно, нужны два начальных условия, чтобы однозначно определить эти константы.
Как правило, аналитическими методами можно решить дифференциальные уравнения только определенного вида. Если уравнение не сводится никакими тождественными преобразованиями ни к одному из известных типов дифференциальных уравнений, то решение такой задачи может быть найдено только специальными численными методами. Рассмотрим наиболее часто применяемые на практике.
|
|
Приближенное решение обыкновенного дифференциального уравнения методом Эйлера.
При обсуждении методов решения задачи Коши ограничимся рассмотрением дифференциального уравнения первого порядка у' = f(х,у), удовлетворяющего начальному условию у(х0) = у0 .
Численное решение задачи состоит в построении таблицы приближенных значений у1, у2, ..., уn, являющихся решениями уравнения у = у(х) в точках х1, х2, ..., хn. Чаще всего точки хi расположены равномерно: хi = х0 + h.i, где i = =1, 2, ..., n. Точки хi называют узлами сетки, h - шагом сетки. Ясно, что h > 0.
Для получения решения задачи Коши воспользуемся разложением функции в ряд Тейлора. Для этого продифференцируем исходное уравнение у' = f(х,у) по х n раз. Тогда с учетом правила дифференцирования сложной функции получим следующие соотношения:
у" = fx (х, у) + fx (х, у) × у';
у"'= fxx (х,у)+2fxy(х,у) × у' + fyy(х,у) × (у') + fy(х,у) × у";. . .
|
|
Полагая х = х0 и у = у0, получаем ряд у'(х0); у"(х0); у"'(х0); ... ; у(n)(х0), который сходится к решению поставленной задачи, т.е.
. (4.1)
Надо заметить, что если |х - х0| больше радиуса сходимости ряда у'(х0); у"(х0); у"'(х0); ...; у(n)(х0), то погрешность |y - у| не стремится к нулю при n , т.е. ряд расходится. Тогда поступают следующим образом. Отрезок [х0,х0+Х], на котором ряд расходится, разбивают на меньшие отрезки [хi, хi+1], где i = 1, 2, ..., n,и получают новые последовательные приближения уj к решению у(хj )при j = 1, 2, ..., m, используя следующую схему:
1) уi считаем найденным (для первого шага им может быть у0);
2) вычисляем в точке хi производные уi(k);
3) полагаем ;
4) считаем уi+1= zi(хi+1);
5) повторяем с первого пункта.
Если же в формуле (4.1) положить n=1, то получим основную расчетную формулу метода Эйлера: уi+1=уi+h´´f(х,у).
Этот метод относится к группе одношаговых, в которых для расчета точки (хi+1, уi+1) требуется информация только о последней вычисленной точке (хi, уi). Он допускает простую геометрическую интерпретацию. Получается, что на каждом новом приближении решение переходит на другую кривую из семейства решений. Этот факт для некоторых дифференциальных уравнений может привести к большим ошибкам и решения задачи Коши начнет расходиться. Заметим, что хотя метод Эйлера относится к итерационным, но ошибки, сделанные на ранних этапах, не уменьшаются из-за неустойчивости вычислительной схемы и ее чувствительности к шагу разбиения отрезка, на котором ищется решение. Чтобы обойти эту трудность, применяют другие вычислительные схемы или методы.
|
|
Для оценки погрешности метода Эйлера на одном из шагов сетки разложим точное решение в ряд Тейлора в окрестности узла хi:
у(хi+1) = у(хi+h) = у(хi) + у'(хi).h + Q(h2) =
= у(хi) + f(хi,уi).h + Q(h2).
Сравнение разложения с основной расчетной формулой метода показывает, что они совпадают до членов первого порядка по h, а погрешность формулы равна Q(h2), т.е. метод Эйлера относится к методам первого порядка.
Формальные параметры процедуры. Входные: х (тип real) - очередное значение хi; h (тип real) - величина шага разбиения отрезка, на котором ищется решение; у (тип real) - предыдущее значение уi; FUNC(х) (тип real) - процедура-функция, вычисляющая значение правой части уравнения по заданным х и у. Результатом работы данной процедуры-функции явится очередное значение уi+1.
Приближенное решение обыкновенного дифференциального уравнения методом Эйлера с уточнением
Метод Эйлера относится к методам первого порядка, потому что обладает не только неустойчивостью решения, но и низкой точностью. Однако в смысле погрешности он дает удовлетворительные результаты даже при не очень больших h. К сожалению, ошибки при вычислении накапливаются от шага к шагу, и к концу отрезка решение содержит значительные погрешности. Заметим также, что источником ошибок очень часто являются и исходные данные. Существуют, однако, методы, позволяющие использовать все достоинства метода Эйлера и одновременно повышающие точность вычислений.
Рассмотрим основную расчетную формулу метода Эйлера: уi+1 = уi+ h.f(х,у). Одна из модификаций метода заключается в том, что сначала вычисляют промежуточные значения
хi+1/2 = хi + h/2и уi+1/2 = уi + fi . h/2
и находят направление поля интегральных кривых в средней точке (хi+1/2, уi+1/2) отрезка [хi, хi+1], т.е.
fi+1/2 = f(х i+1/2, у i+1/2),
а затем полагают уi+1 = уi+1/2+ fi . h/2.
Другой модификацией этого же метода является нелинейная вычислительная схема: зная предыдущее значение уi, вычисляют у ~i+1 как у ~i+1 = уi + h уi , а затем определяют поправку к решению D у ~ i+1по формуле D у ~i+1 = f(хi+1, уi+1) . Решение находится из уравнения
уi+1 = уi + h/2. [ f(хi,уi) + D у ~i+1 ] .
Эту вычислительную схему еще называют методом Эйлера - Коши.
Другая вычислительная схема, которую называют методом Эйлера с упорядочением решения, заключается в том, что каждое значение уi+1= у(хi+1), где у(х) - искомое решение, а х i+1 = х0 + h. (I+1), i = 0, 1, ..., N, определяется по следующей схеме:
1) за начальное приближение берется
,
2) найденное значение уточняется по формуле
, где i = 1, 2, ..., N;
3) уточнение продолжают до тех пор, пока в пределах требуемой точности (заранее заданной) два последовательных приближения не совпадут, т.е.
,
где e - погрешность или любое малое число. После этого полагают уk+1 = уk+1-(i), где уk+1-(i) общая часть полученного ряда решений уk+1(i).
Иногда описанная выше схема называется методом Эйлера с итерационным уточнением решения.
Если после трех итераций решения в пределах требуемой погрешности не совпали, то следует либо уменьшить шаг, либо увеличить погрешность.
Приближенное решение обыкновенного дифференциального уравнения методом Адамса
Если дифференциальное уравнение у'= f(х, у)имеет в правой части сложное аналитическое выражение, то тогда применяют экстраполяционный метод Адамса, который не требует многократного подсчета правой части.
Пусть для дифференциального уравнения заданы начальные условия х = х0 , у = у0. Требуется найти решение уравнения у' = f(х, у) на отрезке [а, b].
Разобьем отрезок [а, b] равномерно на n частей точками хi = х0 + h.i, i = 0, 1, ..., n, h = (b - а)/n. Выберем произвольно элементарный отрезок, на котором проинтегрируем дифференциальное уравнение
или .
Для нахождения производной воспользуемся второй интерполяционной формулой Ньютона, ограничившись разностями третьего порядка с учетом t = (х - хi)/h:
Подставим полученное выражение для у' в интегральное уравнение и, учитывая, что dх = hdt, имеем
Обозначим через qi = уi'; h = f(хi, уi).h, , тогда для любой разности Dm qi = Dm (уi'h) имеем выражение
Dyi = qi + 1/2 . Dqi-1 + 5/12.D2qi-2 + 3/8 . D3qi-3,
используемое для получения решения уравнения
уi+1 = уi + Dуi.
Две последние формулы являются основными в экстраполяционном методе Адамса.
Для начала процесса вычисления нужны четыре начальных значения у0, у1, у2 и у3, которые можно определить любым известным методом. Далее, зная у0, у1, у2 и у3 , находят q0 = =hy0’ = h f(x0, y0); q2 = hy2’ = h f(x2, y2); q3 = hy3’ =h f(x3, y3); q4 = hy4’ = h f(x4, y4)и составляют таблицу конечных разностей величин q.
Таблица 4.1
№ п/п | xi | yi | Dyi | yi’=f(x0,y0) | qi = hyi’ | Конечные разности | ||
0 | x0 | y0 | f(x0,y0) | — | — | — | — | |
1 | x1 | y1 | f(x1,y1) | q0 | Dq0 | D2q0 | D3q0 | |
2 | x2 | y2 | f(x2,y2) | q1 | Dq1 | D2q1 | ||
3 | x3 | y3 | Dy3 | f(x3,y3) | q2 | Dq2 | ||
4 | x4 | y4 | f(x4,y4) | q3 | ||||
5 | x4 | y5 | ... | ... | ||||
... | ... | ... |
Метод Адамса заключается в продолжении данной таблицы разностей с помощью формулы дляDуi. Используя уже вычисленные q3, Dq2, D2q1 иD3q0, расположенные в таблице диагонально, по формуле для Dуi получают, полагая n = 3,
Dу3 = q3 + 0.5Dq2 + (5/12) . D2q1 + (3/8) . D3q0 ,
Dу3 вносят в таблицу и находят у4 = у3+Dу3. Затем, используя х4 и у4 находят f(х4,у4), q4, Dq3, D 2q2 иD3q1, т.е. новую диагональ. По этим данным определяют значение Dу4, которое тут же вносят в таблицу, и находят у5 = у4 + Dу4.
Таблицу продолжают по описанному алгоритму до ее заполнения, вычисляя правую часть формулы при этом только один раз. Чтобы оценить погрешность полученного результата, можно применить правило Рунге или просто следить за третьими разностями D3qi, которые считаются постоянными. Этого можно добиться, выбирая h каждый раз такой, чтобы выражение для оценки погрешности было |D3qi-1 - D3qi| < e. На практике h выбирают из неравенства h4< e, где e- заданная точность решения.
Метод Рунге состоит в том, что сначала находится решение дифференциального уравнения при шаге h, а затем значение h удваивается и находится решение при новом шаге. Погрешность оценивается по формуле
e = (2m - 1) . |yn~ - y~2n| ,
где yn~ - значение приближенного вычисления при двойном шаге; m - порядок метода.
Формальные параметры процедуры. Входные: а, b (тип real) - начало и конец отрезка, на котором ищут решение; h (тип real) - шаг разбиения отрезка [а, b]; ЕРS (тип real) - заранее заданное малое число (используется для определения первых четырех значений уi методом Эйлера с уточнением); уN - начальное у0. Перед началом работы следует определить процедуру-функцию, по которой вычисляют правую часть дифференциального уравнения f (х, у). Выходные: массив уi (тип real) значений уi, найденных методом Адамса.
Приближенное решение обыкновенного дифференциального уравнения методом Рунге - Кутта
Методы Рунге - Кутта образуют семейство методов решения задачи Коши. Oни относятся к многошаговым методам повышенной точности (от второго порядка и выше). Здесь будет рассмотрен одношаговый метод Рунге - Кутта четвертого порядка точности. Заметим попутно, что методы Рунге - Кутта отличаются от разностных тем, что все они допускают вычисление функции не только в заранее заданных узлах сетки, но и в некоторых промежуточных точках, тем самым позволяя уменьшать шаг вычислений во время самого процесса построения решетки.
Итак, пусть требуется найти решение дифференциального уравнения у' = f(х, у), удовлетворяющего начальному условию у(х0) = у0 .
Численное решение задачи состоит в построении таблицы приближенных значений у0, у1, ..., уn решения уравнения у' = f(х,у) в точках х0, х1, ..., хn, которые называют узлами сетки. Для краткости обозначим хi = х0+ ih; уi = у(хi), где h - шаг сетки (ясно, что h > 0).
Обычно методом Рунге-Кутта в литературе называют метод, согласно которому у искомой функции вычисляют по формулам:
уi+1 = уi + D/6,
где обозначено:
D = k1+2k2 +2k3 +k4 ; k1 = hf(xi, yi);
k2 = h f(xi + h/2 ,yi + k1/2); k3 = h f(xi + h/2 ,yi + k2/2);
k4 = h f(xi + h ,yi + hk3).
Погрешность метода на одном шаге сетки равна М.h5 , где М - некоторое число. Но на практике даже приближенно очень трудно оценить М, поэтому при оценке погрешности используют правило Рунге. Для чего сначала проводят вычисления с шагом h, а затем повторяют их с числом h/2. Если уi(h) - приближение, вычисленное с шагом h, а уi(h/2) - с шагом h/2, то справедлива оценка
.
Формальные параметры процедуры. Входные: N (тип integer) - количество разбиений отрезка [а, b]; Х, Y (тип real) - начальные значения, соответствующие x0 и y0; h (тип real) - шаг интегрирования системы; FUNС - имя процедуры-функции, по которой вычисляют значения правой части уравнения у' = f(х, у). Выходные: YR (тип real) - массив, содержащий решения системы в точках хi = х0 + hi.
Задание
Написать программу для интегрирования обыкновенного дифференциального уравнения (методы задаются преподавателем).
Варианты
1. ,
2. ,
3. ,
4. ,
5. ,
6. ,
7. ,
8. ,
9. ,
10. ,
Требования к программе
1. Программа должна быть написана в среде Visual Studio или Delphi.
2. В программе должно быть предусмотрено следующее:
- ввод исходных данных и вывод результата в удобной для пользователя форме;
- возможность пошагового просмотра решения задачи;
- справка с комментариями по работе с программой;
- корректировка исходных данных;
- сохранение исходных данных и результата в файле.
Порядок выполнения работы
1. Получить задание у преподавателя.
2. Разработать алгоритм решения задачи.
3. Реализовать полученный алгоритм.
4. Проанализировать результаты работы алгоритма.
5. Оформить отчет по лабораторной работе.
Содержание отчета
1. Номер и тема лабораторной работы.
2. Цель выполнения работы.
3. Описание алгоритма.
4. Руководство пользователю.
5. Пример работы программы.
6. Анализ полученных результатов и вывод по работе.
7. Текст программы с комментариями привести в приложении.
Библиографический список
1. Белашов В.Ю., Чернова Н.М. Эффективные алгоритмы и программы вычислительной математики. Магадан: СВКНИИ ДВО РАН, 1997. 160 с.
2. Амосов А.А., Дубинский Ю.А., Копченова Н.В. Вычислительные методы для инженеров. М.: Высш. шк., 1994. - 544 с.
3. Турчак Л.И. Основы численных методов. М.: Наука, 1987. – 320 с.
4. Демидович Б.Н., Марон И.А. Основы вычислительной математики. М.: Высш. шк., 1994. - 172 с.
МЕТОДИЧЕСКИЕ УКАЗАНИЯ
к выполнению лабораторных работ
по дисциплине «Вычислительная математика»
для студентов специальностей
230104 – «Системы автоматизированного проектирования»,
230202 – «Информационные технологии в образовании»
всех форм обучения
Составители
Бобров Александр Иванович
Федорков Евгений Дмитриевич
Чекменев Анатолий Николаевич
В авторской редакции
Компьютерный набор А.И. Боброва
Подписано в печать _____________.
Формат 60х84/16. Бумага для множительных аппаратов.
Усл. печ. л. __. Уч.-изд. л.__. Тираж ___ экз. “C”.
Зак. №
ГОУВПО «Воронежский государственный технический университет»
394026 Воронеж, Московский просп., 14
Дата добавления: 2018-04-15; просмотров: 212; Мы поможем в написании вашей работы! |
Мы поможем в написании ваших работ!