Некоторые виды уравнений, интегрируемых в квадратурах
По этой ссылке вы найдёте полный курс лекций по математике:
общем случае, даже зная, что решение уравнения существует, отыскать его довольно трудно. Однако существуют некоторые виды дифференциальных уравнений, методы получения решений которых особенно просты (при помощи интегралов от элементарных функций). Рассмотрим некоторые из них. 5.1. Уравнения с разделяющимися переменными Уравнение вида называется дифференциальным уравнением с разделенными переменными. Здесь f\(y), /2 (я) — известные непрерывные функции своих аргументов.
Покажем, как найти решение этого уравнения. Пусть F\(y) и F2(x) — первообразные функции f\(y) и /2(2) соответственно. Равенство (1) равносильно тому, что дифференциалы этих функций должны совпадать Отсюда следует, что где С — произвольная постоянная. Разрешая последнее уравнение (2) относительно у, получим функцию (может быть, и не одну) которая обращает уравнение (1) в тождество и значит, является его решением Например, — уравнение с разделенными переменными.
Записав его в виде и интегрируя обе части, найдем общий интеграл данного уравнения: Уравнение вида в котором коэффициенты при дифференциалах распадаются на множители, зависящие только от х и только от уу называется дифференциальным уравнением с разделяющимися переменными, так как путем деления на 4>\(y)fi(x) £ 0 оно приводится к уравнению с разделенными переменными Пример 1.
Проинтегрировать уравнение Некоторые виды уравнений, интегрируемых в квадратурах Уравнения с разделяющимися переменными Уравнения, однородные относительно х и у Линейные дифференциальные уравнения Уравнение Бернулли Уравнения в полных дифференциалах Деля обе части уравнения на , приведем его к виду Интегрируя обе части полученного равенства, найдем откуда Заметим, что деление на может привести к потере решений, обращающих в нуль произведение .
Например, разделяя переменные в уравнении получаем а после интегрирования — откуда (здесь С может принимать как положительные, так и отрицательные значения, но С Ф 0). При делении на у потеряно решение которое может быть включено в общее решение у = Сх, если постоянной С разрешить принимать значение С = 0. Если считать переменные х и у равноправными, то уравнение теряющее смысл при х = 0, надо дополнить уравнением которое имеет очевидное решение х = 0. где /j (ж, у) = у^у, используя уравнение (4′) там, где уравнение (4) не имеет смысла, а уравнение (4′) имеет смысл.
Некоторые дифференциальные уравнения путем замены переменных могут быть приведены к уравнениям с разделяющимися переменными. Например, уравнение вида где f(z) — непрерывная функция своего аргумента, а, Ь, с — постоянные числа, подстановкой z = ах + by + с преобразуется в дифференциальное уравнение с разделяющимися переменными: откуда После интегрирования получаем Заменяя в последнем соотношении найдем общий интеграл уравнения (5). Пример 2.
Проинтегрировать уравнение В общем случае наряду с дифференциальным уравнением следует рассматривать уравнение 4 Положим z = x + y, тогда откуда Интегрируя, находим или Подставляя вместо z величину х + у, получаем общее решение данного уравнения Пример 3. Известно, что скорость радиоактивного распада пропорциональна количеству х еще не распавшегося вещества. Найти зависимость х от времени t, если в начальный момент t = to имелось х = х0 вещества.
Дифференциальное уравнение процесса dx Здесь к > 0 — постоянная распада — предполагается известной, знак указывает на уменьшение х при возрастании t. Разделяя переменные в уравнении и интегрируя, получаем откуда х = Се
и. Учитывая начальное условие = xq, находим, что С — х0ек1с, поэтому Любой процесс (не только радиоактивный распад), при котором скорость распада пропорциональна количеству еще не прореагировавшего вещества, описывается уравнением (*).
Уравнение отличающееся лишь знаком правой части от уравнения (*), описывает лавинообразный процесс размножения, например «размножение» нейтронов в цепных ядерных реакциях или размножение бактерий в предположении, что скорость их размножения пропорциональна наличному числу бактерий. Решение уравнения , удовлетворяющее условию имеет вид и в отличие от решения уравнения ) возрастает с возрастанием t.
Уравнения можно объединить в одно которое дает простейшую математическую модель динамики популяций (совокупности особей того или иного вида растительных или животных Организмов). Пусть y(t) — число членов популяции в момент времени t. Если предположить, что скорость изменения популяции пропорциональна величине популяции, то мы приходим к уравнению . Положим к —т — я, где m — коэффициент относительной скорости рождаемости, an — коэффициент относительной скорости умирания.
Тогда к > 0 при п и к при тп. Если в момент t = О величина популяции равна уо, то уравнение приводит к экспоненциальному закону изменения популяции Предположение, что величины шип являются постоянными, не выполняется для больших популяций. Действительно, большое число членов популяции приводит к уменьшению соответствующих ресурсов, что снижает скорость рождаемости и увеличивает скорость умирания. Это можно задать простейшими законами положительные постоянные.
Возможно вам будут полезны данные страницы:
Тогда Уравнение динамики популяции в этой модели имеет вид Это так называемое логистическое уравнение — фундаментальное уравнение в демографии и в математической теории экологии. Оно применяется в математической теории распространения слухов, болезней и других проблемах физиологии и социологии. Разделяя переменные в последнем уравнении, получаем откуда и далее .
Потенцируя Некоторые виды уравнений, интегрируемых в квадратурах Уравнения с разделяющимися переменными Уравнения, однородные относительно х и у Линейные дифференциальные уравнения Уравнение Бернулли Уравнения в полных дифференциалах и выражая у через t, окончательно получаем Считая, , найдем уравнение логистической кривой При получаем, что Логистическая кривая содержит два параметра Л и а. Для их определения надо иметь два дополнительных значения y(t) при каких-то 5.2.
Уравнения, однородные относительно х и у Функция f(x, у) называется однородной функцией п-го измерения относительно переменных х иу, если при любом допустимом t справедливо тождество .
Например, для функции так что однородная функция относительно переменных х и у второго измерения. Для функции имеем есть однородная функция нулевого измерения. Дифференциальное уравнение первого порядка называется однородным относительно х и у, если функция f(x, у) есть однородная функция нулевого измерения относительно переменных х и у.
Пусть имеем дифференциальное уравнение однородное относительно переменных хну. Положив t = ^ в тождестве получим т. е. однородная функция нулевого измерения зависит только от отношения аргументов. Обозначая / (I, |) через (*), видим, что однородное относительно переменных ж и у дифференциальное уравнение всегда можно представить в виде При произвольной непрерывной функции переменные не разделяются. Введем новую искомую функцию u(z) формулой и = J, откуда у = хи.
Подставляя выражение ^ = и 4- в уравнение (6), получаем Деля обе части последнего равенства на и интегрируя, находим Заменяя здесь и на его значение \, получаем общий интеграл уравнения (6). Пример 4. Проинтегрировать уравнение Положим уравнение преобразуется к виду Интегрируя, найдем или Пример 5. Найти форму зеркала, собирающего пучок параллельно падающих на него лучей в одну точку.
Прежде всего, зеркало должно иметь форму поверхности вращения, так как только для поверхности вращения все нормали к поверхности проходят через ось вращения. Выберем прямоугольную декартову систему координат так, чтобы лучи были параллельны оси Ох и чтобы точкой, в которой собирались бы отраженные лучи, явилось бы начало координат. Найдем форму сечения зеркала плоскостью хОу. Пусть уравнение сечения есть ) (рис.6).
В точке Af(x,y) падения луча L на зеркало проведем касательную BN к сечению и обозначим ее угол с осью Ох через а. Пусть N — точка пересечения этой касательной с осью Ох. По закону отражения углы NMO и BML равны. Нетрудно видеть, что угол МОР равен 2а. Так как tga = у , tg2a = J то во всякой точке кривой у — ip(x) выполняется соотношение — дифференциальное уравнение, определяющее требуемый ход луча. Разрешая это уравнение относительно производной, получаем два однородных уравнения:
Первое из них путем замены | = и преобразуется к виду — произвольная постоянная. Решение неоднородного уравнения (10) ищем в виде — новая неизвестная функция. Вычисляя производную и подставляя значения и у в исходное уравнение (10), получаем dx откуда — новая произвольная постоянная интегрирования. Следовательно, Это есть общее решение линейного неоднородного дифференциального уравнения (10).
В формуле (14) общего решения неопределенные интегралы можно заменить определенными интегралами с переменным верхним пределом: Здесь С = у(жо) = уо, поэтому общее решение уравнения (10) можно записать в виде (15) где роль произвольной постоянной играет начальное значение уо искомой функции у (ж). Формула (15) является общим решением уравнения (10) в форме Коши.
Отсюда следует, что если р(х) и q(x) определены и непрерывны в интервале то и решение у(х) уравнения (10) с любыми начальными данными у(ж0) = уо будет непрерывным и даже непрерывно дифференцируемым при всех конечных значениях ж, так что интегральная кривая, проходящая через любую точку (жо, уо), будет гладкой кривой в интервале Поэтому сначала интегрируем соответствующее однородное уравнение общее решение которого имеет вид Пример б.
Проинтегрировать уравнение Ч Однородное уравнение соответствующее данному, проинтегрируем, разделяя переменные: Решение исходного уравнения будем искать в виде где С(х) — неизвестная функция. Находя ^ и подставляя , последовательно получаем: где С — постоянная интегрирования. Из формулы находим общее решение уравнения Частное решение неоднородного уравнения легко усматривается. Вообще, если удается «угадать» частное решение линейного неоднородного уравнения, то разыскание его общего решения значительно упрощается. Пример 7.
Рассмотрим дифференциальное уравнение, описывающее изменение силы тока при замыкании цепи постоянного электрического тока. 4 Если R — сопротивление цепи, Е — внешняя ЭДС, то сила тока / = /(f) постепенно возрастает от значения, равного нулю, до конечного стационарного значения ^. Пусть L — коэффициент самоиндукции цепи, роль которой такова, что при всяком изменении си- лы тока в цепи появляется электродвижущая сила, равная L-^ и направленная противоположно внешней ЭДС.
На основании закона Ома, по которому в каждый момент t произведение силы тока на сопротивление равно фактически действующей ЭДС, получаем ) Уравнение (*) есть линейное неоднородное уравнение относительно I(t). Нетрудно видеть, что его частным решением является функция Общее решение соответствующего однородного уравнения откуда общее решение неоднородного уравнения («):
При t = 0 имеем , поэтому так что окончательно Отсюда видно, что сила тока при включении асимптотически приближается при к своему стационарному значению ожет быть проинтефировано также следующим приемом. Будем искать решение у(х) уравнения (10) в виде Линейное неоднородное дифференциальное уравнение где и<х) и v(x) — неизвестные функции, одна из которых, например v(x), может быть выбрана произвольно.
| Подставляя у(х) в форме (16) в уравнение |
(10), после элементарных преобразований получим Выберем в качестве v(ar) любое частное решение v(«) ^ 0 уравнения Тогда в силу (17) для и(ж) получим уравнение которое без труда интегрируется в квадратурах. Зная v(x) и и(х), найдем решение у(х) уравнения (10). Пример. Найти общее решение уравнения Будем искать решение у(х) данного линейного неоднородного уравнения в виде Подставляя исходное уравнение, получим.
Определим функцию v(x) как решение уравнения Разделяя переменные, найдем Выберем любое частое решение, например, отвечающее ) получим Для общего решения исходного уравнения получаем выражение Преимущество метода вариации постоянной заключается в том, что он переносится на линейные неоднородные дифференциальные уравнения высшего порядка. Уравнение Бернулли Некоторые дифференциальные уравнения путем замены переменных могут быть сведены к линейным.
К числу таких уравнений относится уравнение Бернуми Уравнение это предложено Я. Бернулли в 1695 г., метод решения опубликовал И. Бернулли в 1697 г. При а = 1 получаем однородное линейное уравнение При а = 0 — неоднородное линейное уравнение Поэтому будем предполагать, что (для а нецелого считаем, что у > 0). Подстановкой z = y_a+1 уравнение Бернулли приводится к линейному уравнению относительно функции z(x). Однако уравнение Бернулли можно проинтегрировать сразу методом вариации постоянной.
Это делается так. Сначала интегрируем уравнение Его общее решение Решение уравнения Бернулли будем искать в виде где С(х) — новая неизвестная функция. Подставляя это выражение для у(х) в уравнение Бернулли, получаем — уравнение с разделяющимися переменными относительно С(х). Интегрируя это уравнение, находим где С — постоянная интегрирования. Тогда из формулы получаем обший интеграл уравнения Бернулли Замечание.
При а > 0 уравнение Бернулли имеет очевидное решение у = 0. Для интегрирования уравнения Бернулли можно также воспользоваться подстановкой Некоторые виды уравнений, интегрируемых в квадратурах Уравнения с разделяющимися переменными Уравнения, однородные относительно х и у Линейные дифференциальные уравнения Уравнение Бернулли Уравнения в полных дифференциалах где в качестве берется любое нетривиальное решение уравнения а функция и(ж) определяется как решение уравнения Пример.
Найти решение уравнения Бернулли Ищем решение у(х) уравнения в виде Подставляя исходное уравнение, получим . Выберем в качестве v(x) какое-нибудь ненулевое решение уравнения и проинтегрируем его, Поскольку нас интересует какое угодно частное решение, положим . возьмем v = Тогда для и(х) получим уравнение интегрируя которое, найдем Общее решение у(х) исходного уравнения определится формулой 5.5.
Уравнения в полных дифференциалах Уравнение называется уравнением в полных дифференциалах, если левая часть уравнения представляет собой полный дифференциал некоторой функции гх(ж, у) двух независимых переменных ж и у, т. е. В этом случае «(ж, у) — С будет общим интегралом дифференциального уравнения (18). Будем предполагать, что функции М(ж, у) и N(x, у) имеют непрерывные частные производные соответственно по у и по ж в некоторой односвязной области D на плоскости хОу. Теорема 4.
Для того чтобы левая часть М(ж, у) dx + N(xt у) dy уравнения (18) была полным дифференциалом некоторой функции и(х, у) двух независимых переменных х и у, необходимо и достаточно, чтобы выполнялось тождество 4 Необходимость. Предположим, что левая часть уравнения (18) есть полный дифференциал некоторой функции Дифференцируем первое соотношение по у, а второе по х: ду вудх дх дхду Отсюда, в силу равенства смешанных производных, вытекает тождество Необходимость (19) доказана.
Достаточность. Покажем, что условие (19) является и достаточным, а именно, предполагая его выполненным, найдем функцию tx(x, у) такую, что или, что то же, Найдем сначала функцию и(ху у), удовлетворяющую первому условию (20). Интегрируя это равенство по х (считаем у постоянной), получаем где произвольная функция от у. Подберем так, чтобы частная производная по у от функции и, определяемой формулой (21), была равна N(x> у). Такой выбор функции (р(у) при условии (19) всегда возможен.
В самом деле, из (21) имеем Приравняв правую часть полученного равенства к N(x, у), найдем Левая часть последнего равенства не зависит от х. Убедимся в том, что при условии (20) в его правую часть также не входит х. Для этого покажем, что частная производная по х от правой части (22) тождественно равна нулю. Имеем Теперь, интегрируя равенство (22) по у, получим, что где С — постоянная интегрирования.
Подставляя найденное значение для в формулу (21), получим искомую функцию полный дифференциал которой, как нетрудно проверить, равен Приведенный прием построения функции м(х, у) составляет метод интегрирования уравнения (18), левая часть которого есть полный дифференциал. Пример 8. Проверить, что уравнение является уравнением в полных дифференциалах, и проинтегрировать его. А В данном случае откуда Следовательно, уравнение есть уравнение в полных дифференциалах.
Теперь находим и или Находя от функции и из и приравнивая ^ функции , получаем следовательно, Подставив найденное выражение для , найдем Таким образом, — общий интеграл исходного уравнения. Иногда можно найти такую функцию , что будет полным дифференциалом, хотя М dx + N dy может им и не быть. Такую функцию у) называют интегрирующим множителем.
Можно показать, что для уравнения первого порядка при определенных условиях на функции М(х, у) и N(x, у) интегрирующий множитель всегда существует, но отыскание его из условия в общем случае сводится к интегрированию уравнения в частных производных, что составляет, как правило, задачу еще более трудную. Задача. Найти интегрирующий множитель для линейного дифференциального уравнения dy Указание. Искать множитель в виде ц = ц(х).
Присылайте задания в любое время дня и ночи в ➔ 
Официальный сайт Брильёновой Натальи Валерьевны преподавателя кафедры информатики и электроники Екатеринбургского государственного института.
Все авторские права на размещённые материалы сохранены за правообладателями этих материалов. Любое коммерческое и/или иное использование кроме предварительного ознакомления материалов сайта natalibrilenova.ru запрещено. Публикация и распространение размещённых материалов не преследует за собой коммерческой и/или любой другой выгоды.
Сайт предназначен для облегчения образовательного путешествия студентам очникам и заочникам по вопросам обучения . Наталья Брильёнова не предлагает и не оказывает товары и услуги.
20. Численное интегрирование методом квадратур
Численное интегрирование методом квадратур
Приведенные ниже функции осуществляют интегрирование и двойное интегрирование, используя квадратурную формулу Симпсона или метод Гаусса-Лобатто. Квадратура — численный метод нахождения площади под графиком функции/(т), т. е. вычисление определенного интеграла вида
В приведенных ниже формулах подынтегральное выражение fun обычно задается в форме дескриптора функции, поэтому с дидактическими целями используем нотацию @fun.
Функции quad и quadl используют два различных алгоритма квадратуры для вычисления определенного интеграла. Функция quad выполняет интегрирование по методу низкого порядка, используя рекурсивное правило Симпсона [4, 40]. Но она может быть более эффективной при негладких подынтегральных функциях или при низкой требуемой точности вычислений. Алгоритм quad в MATLAB 6 изменен по сравнению с предшествовавшими версиями, точность по умолчанию по сравнению с версиями 5.3х повышена в 1000 раз (с 10- 3 до 10- 6 ). Новая функция quadl (квадратура Лобатто) использует адаптивное правило квадратуры Гаусса— Лобатто очень высокого порядка. Устаревшая функция quads выполняла интегрирование, используя квадратурные формулы Ньютона—Котеса 8-го порядка [40]. Достижимая точность интегрирования гладких функций в MATLAB 6 поэтому также значительно выше, чем в предшествующих версиях.
quad(@fun,a.b) — возвращает численное значение определенного интеграла от заданной функции @fun на отрезке [а Ь]. Используется значительно усовершенствованный в MATLAB 6 адаптивный метод Симпсона;
quad(@fun,a.b.tol) — возвращает численное значение определенного интеграла с заданной относительной погрешностью tol. По умолчанию to1=l.e-6. Можно также использовать вектор, состоящий из двух элементов tol =[rel_tol abs_tol], чтобы точно определить комбинацию относительной и абсолютной погрешности;
quad(@fun,a.b,tol .trace) — возвращает численное значение определенного интеграла и при значении trace, не равном нулю, строит график, показывающий ход вычисления интеграла;
quad(@fun,a,b.tol,trace,PI,P2. ) — возвращает численное значение определенного интеграла по хот подынтегральной функции fun, использует дополнительные аргументы P1, P2. которые напрямую передаются в подынтегральную функцию: G=fun(X.Pl,P2. ). Примеры:
dblquad(@fun,inmin,inmax.outmin,outmax) — вычисляет и возвращает значение двойного интеграла для подынтегральной функции fun (Inner, outer), по умолчанию используя квадратурную функцию quad. Inner — внутренняя переменная, изменяющаяся на закрытом интервале от inmin до inmax, a outer — внешняя переменная, изменяющаяся на закрытом интервале от outmin до outmax. Первый аргумент @fun — строка, описывающая подынтегральную функцию. Этс может быть либо дескриптор функции, либо объект inline (в последнем случае символ «@» в ее записи отсутствует). Обычная запись в апострофах тепер недопустима. Эта функция должна быть функцией двух переменных вид. fout=fun(inner.outer). Функция должна брать вектор inner и скаляр outer возвращать вектор fout, который является функцией, вычисленной в outer и каждом значении inner [ Функция inime(‘expr’, ‘argl’. ‘argn’) так же создает объект, но без дескриптора, ‘ехрг’ — выражеши Строки ‘argx’ —входные аргументы. При их отсутствии по умолчанию подставляется х. Если вмест ‘arg’ — скаляр, то он означает количество дополнительных переменных Р. Примеры запис; g = inline(exp); g — inline(‘t A 2’); g = inline(‘sin(2*pi*f + theta)’); g = inline(‘sin(2*pi*f + theta)’, : ‘theta’); g — inline(‘x A Pl+P2’, 2). — Примеч. ред. ];
dblquad(@fun,inmin.inmax.outmin,oiitmax,tol .trace) — передает в функцию dblquad параметры tol и trace. Смотрите справку по функции quad для получения информации о параметрах to! и trace;
dblquad(@fun,inmin,inmax,outmin,outmax.tol .trace,order) — передает параметры tol и trace для функции quad или quadl в зависимости от значения строки order. Допустимые значения для параметра order — @quad , @quadl или имя любого определенного пользователем квадратурного метода с таким же вызовом и такими же возвращаемыми параметрами, как у функций quad и quad!. (Например, при проверке старых программ можно использовать @quad8 для большей совместимости с прежними версиями MATLAB). По умолчанию (без параметра order) вызывается @quad. поскольку подинтегральные функции могут быть негладкими.
Пример: пусть m-файл integl.m описывает функцию 2*y*sin(x)+x/2*cos(y), тогда вычислить двойной интеграл от той функции можно следующим образом:
И. О. Арушанян. Применение метода квадратур для численного решения интегральных уравнений второго рода. Учебное пособие
- Анна Калантаева 3 лет назад Просмотров:
1 МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ им. М.В. ЛОМОНОСОВА Механико математический факультет Кафедра вычислительной математики И. О. Арушанян Применение метода квадратур для численного решения интегральных уравнений второго рода Учебное пособие Москва, 18
2 И О. Арушанян ПРИМЕНЕНИЕ МЕТОДА КВАДРАТУР ДЛЯ ЧИСЛЕННОГО РЕШЕНИЯ ИНТЕГРАЛЬНЫХ УРАВНЕНИЙ ВТОРОГО РОДА Учебное пособие ОГЛАВЛЕНИЕ Введение Численное интегрирование Квадратурные формулы интерполяционного типа Элементарные формулы трапеций, средних прямоугольников и Симпсона Составные формулы трапеций, средних прямоугольников и Симпсона Построение формул Ньютона Котеса методом неопределенных коэффициентов Квадратурные формулы Гаусса Правило Рунге практической оценки погрешности квадратурных формул Процесс Эйткена оценки фактической точности квадратурной формулы Способы построения практическихалгоритмов численного интегрирования. Адаптивные алгоритмы Простейший неадаптивный) алгоритм Модификация простейшего алгоритма Адаптивный алгоритм последовательного передвижения слева-направо Адаптивный алгоритм с контролем точности по глобальной ошибке Интегральные уравнения Фредгольма второго рода Метод квадратур решения интегральных уравнений Фредгольма второго рода Сходимость метода квадратур Практический алгоритм численной реализации метода квадратур. 37
3 .. Метод последовательных приближений Первый алгоритм численной реализации метода последовательных приближений Второй алгоритм численной реализации метода последовательных приближений Интегральные уравнения Вольтерра второго рода Метод квадратур решения интегральных уравнений Вольтерра второго рода Сходимость метода квадратур Построение приближенного решения в виде непрерывной функции Способы построения квадратурных формул для решения уравнений Вольтерра второго рода Практический алгоритм построения приближенного решения с заданной точностью Варианты выполнения задания. Содержание отчета Тестовые задачи для отладки программ Литература
4 Введение Данное учебное пособие посвящено рассмотрению метода квадратур для решения линейных интегральных уравнений второго рода, предлагаемого для практической реализации на ЭВМ при выполнении заданий студенческого практикума на механико-математическом факультете Московского университета. В пособии рассматриваются два типа интегральных уравнений. Это интегральные уравнения Фредгольма второго рода ux) b Kx,s)us)ds = fx), x [, b], ) и интегральные уравнения Вольтерра второго рода ux) x Kx,s)us)ds = fx), x [, b]. ) Здесь ux) искомая функция с областью определения [, b]. Заданные функции Kx, s) и fx) называются соответственно ядром и правой частью интегрального уравнения. Формально уравнение **) является частным случаем уравнения *), если ядро уравнения отвечает условию Kx, s) при s > x. Однако это свойство функции K имеет настолько глубокие следствия, что делает теории интегральных уравнений **) и *) принципиально различными. В рамках пособия имеет смысл остановиться на том, какое влияние оказывают эти отличительные особенности на методы численного решения интегральных уравнений. Из всего многообразия методов численного решения интегральных уравнений вида *) и **) в пособии рассматривается только метод квадратур и некоторые его модификации. Входящие в уравнения интегралы заменяются квадратурными суммами, а полученные конечные соотношения могут быть как вспомогательными, так и имеющими самостоятельный характер в качестве окончательных расчетных выражений. В уравнениях Вольтерра **) значение искомой функции u в точке x определяется ее значениями в предшествующих точках s x не влияют на ux). Для приближенного решения уравнения на отрезке [, b] берут сетку из точек s 5 чтобы по известным значениямu n,u n 1. u n 1 могло быть построено уравнение для нахождения u n. В случае уравнения Фредгольма *) ситуация принципиально иная. Интеграл в *) зависит от значений неизвестной функции ux) на всем отрезке. Поэтому в каждое уравнение аппроксимирующей системы должны входить значения u n во всех точках сетки и система должна решаться не последовательно, а в целом. В учебной литературе описание численных методов решения интегральных уравнений состоит, как правило, в построении аппроксимирующей системы линейных уравнений, алгоритме ее решения и далеко не всегда) оценке погрешности в точках сетки. В данном пособии в качестве задания практикума для решения интегральных уравнений предложено несколько численных методов, связанных общим принципом построением последовательности непрерывных функций, сходящихся в интегральной норме к точному решению. Для реализации подобных методов существенное значение имеет изучение алгоритмов численного интегрирования с автоматическим выбором распределения узлов, описанию которых посвящен первый раздел пособия. Во втором и третьем разделах рассматриваются методы численного решения интегральных уравнений второго рода. В четвертом разделе излагаются этапы выполнения задания и дан перечень вариантов. Заключает пособие набор тестовых задач для отладки программ студенческого практикума. 5
6 1. ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ Пусть требуется вычислить значение определенного интеграла I = b fx)dx, где функция fx) непрерывна на отрезке [, b]. Числа и b называют пределами интегрирования; они могут быть как конечными, так и бесконечными тогда интеграл называется несобственным). Как известно, через элементарные функции интегралы выражаются далеко не всегда. По этой и другим причинам для их вычисления применяют численные методы, основанные на замене fx) такой аппроксимирующей функцией, чтобы интеграл от нее вычислялся в элементарных функциях достаточно просто Квадратурные формулы интерполяционного типа В качестве аппроксимирующих функций возьмем интерполяционные многочлены Лагранжа L n x). Тогда искомый интеграл заменяется линейной комбинацией значений подынтегральной функции fx): I = b fx)dx = n c i fx i )+R = n c i f i +R. 1) Выписанное соотношение называется квадратурной формулой или просто квадратурой), в которой величины x i называют узлами, c i весами, а R погрешностью остаточным членом). Таким образом, интеграл приближенно заменяется суммой, похожей на интегральную сумму, причем узлы и веса не зависят от fx). Рассмотрим квадратуру 1) более подробно. Пусть узлы x i образуют равномерную сетку с шагом h: x i+1 = x i + h, i = 1. n. Если пределы интегрирования и b входят в состав сетки, то h = b )/n 1) и x 1 =, x i = x 1 +i 1)h, x n = b. В этом случае формула 1) называется формулой замкнутого типа и имеет по крайней мере два узла. Исходный отрезок интегрирования разбивается на n 1 подотрезков длины h. Если и b не входят в состав сетки, то h = b )/n+1) и x i = +ih. Тогда формула 1) называется формулой открытого типа и может иметь один узел. Количество подотрезков равно n+1. Если либо, либо b включены в состав сетки, то соответствующие квадратуры носят название полуоткрытых или полузамкнутых). 6
7 Теперь заменим fx) на многочлен ЛагранжаL n x), интерполирующий fx) в точках x i,f i ) на отрезке [,b]: fx) = L n x) = n f i l i x)+ ω nx) n! f n) ξ), ξ [,b], где ω n x) = x x 1 ). x x n ) и l i x) = j=1 j i x x j x i x j = = x x 1). x x i 1 )x x i+1 ). x x n ) x i x 1 ). x i x i 1 )x i x i+1 ). x i x n ). Функцииl i x) называют фундаментальными многочленами Лагранжа. Подставив выписанное приближенное представление fx) в 1), получим где I = n b = b ) l i x)dx f i + n c i f i +R, b ω n x) n! f n) ξ)dx = c i = 1 b l i x)dx, R = 1 b ω n x)f n) ξ)dx. b n! Если сетка
8 Действительно, при нечетном n один из узлов сетки совпадает с серединой отрезка интегрирования x = + b)/, а остальные узлы лежат симметрично относительно x. Рассмотрим многочлен qx) = x x) n. Этот многочлен является нечетным относительно x; следовательно b qx) = n c i q i =. Отсюда R =, т.е. формула ) точна для qx). Покажем, что эта формула точна и для любого многочлена p n x) степени n: p n x) = 1 x n + x n n x+ n+1. Для этого представим p n x) в виде p n x) = 1 x x) n +p n 1 x) = 1 qx)+p n 1 x) и подставим его в ). Поскольку формула ) по построению точна для p n 1 x) и точна для qx), то она точна и для p n x). Главный член погрешности формулы Ньютона Котеса с n узлами при условии достаточной гладкости fx) имеет порядок O h [n 1)/]+3), где [ ] целая часть числа. Порядок по h главного члена погрешности называется порядком точности сходимости) квадратурной формулы. Для формул замкнутого типа коэффициенты Котеса c i положительны при1 n 8, а приn = 9 иn 11 среди них есть отрицательные, что приводит к увеличению ошибок, содержащихся в fx). Действительно, пусть ошибка задания функции fx) в каждом узле сетки оценивается сверху по модулю некоторой величиной ε. Погрешность, которая может получиться в сумме c i f i, оценивается величиной ε c i. Поскольку c i = 1, то наличие отрицательных c i приводит к увеличению c i. Коэффициент увеличения ошибок иллюстрируется следующими цифрами: ci 1.45, n = 9; ci 3.1, n = 11; ci 8.3, n = 15; ci 56, n =. Значения c i при больших n быстро растут по абсолютной величине, а c i будет большой и быстро возрастающей величиной. Коэффициент увеличения ошибок становится неприемлемым, а соответствующие формулы непригодными для расчетов. Для формул открытого типа коэффициенты Котеса положительны при n = 1,,4, а при остальных n среди них есть отрицательные. 8
9 1.. Элементарные формулы трапеций, средних прямоугольников и Симпсона Рассмотрим первые три формулы Ньютона Котеса. Сначала приблизим fx) на [,b] многочленом Лагранжа L x) с узлами x 1 = и x = b. Это означает, что вместо кривой fx) мы взяли многочлен первой степени, проходящий через точки,f)) и b,fb)). Теперь искомый интеграл, равный площади криволинейной фигуры, заменим на площадь трапеции с основаниями f) и fb) и высотой h = b : I T = h f)+fb) = b f)+fb) ). 3) Мы получили квадратурную формулу трапеций, которую можно также получить из ) при n = с сеткой для формул замкнутого типа. Заметим, что алгебраический порядок точности этой формулы равен единице, т.к. середина отрезка [,b] не входит в состав сетки. Если на [, b] взять единственный узел квадратурной формулы, то fx) аппроксимируется многочленом Лагранжа L 1 x), т.е. многочленом нулевой степени. Возьмем в качестве этого узла середину отрезка x = +b)/. Тогда получим формулу средних прямоугольников I P = hf x). 4) Геометрический смысл этой формулы состоит в том, что площадь криволинейной фигуры заменяется на площадь прямоугольника с основанием h = b и высотой f x). Ту же самую формулу можно получить из ) при n = 1 с сеткой для формул открытого типа. Алгебраический порядок формулы средних прямоугольников равен единице, поскольку середина отрезка интегрирования входит в состав узлов сетки. Теперь получим выражения для остаточных членов R T = I T и R P = I P. Для этого разложим fx) в ряд Тейлора относительно точки x = + b)/ в предположении о достаточной гладкости функции fx): fx) = f x)+x x)f x x) x)+ f x)+ x x)3 + f x x)4 x)+ f IV x) ) Воспользовавшись этим разложением, получим представление остаточного члена формулы средних прямоугольников: I = hf x)+ h3 4 f x)+ h5 19 fiv x)+ = P +R P, 6) 9
10 поскольку b x x) j dx = h, j = ;, j = 1; h3 h5, j = ;, j = 3; 1 8, j = 4. Подставим в 5) x = и x = b и учтем, что x = h/ и b x = h/: f) = f x) h f x)+ h 8 f x) h3 48 f x)+ h4 384 fiv x)+, fb) = f x)+ h f x)+ h 8 f x)+ h3 48 f x)+ h4 384 fiv x)+. Сложим эти два равенства: Отсюда получим f)+fb) = f x)+ h 8 f x)+ h4 384 fiv x)+. hf x) = h f)+fb) Подставим это выражение вместо hf x) в 6): I = h f)+fb) h3 8 f x) h5 384 fiv x)+. h3 1 f x) h5 48 fiv x)+ = T +R T. 7) Итак, главные члены погрешности в формулах средних прямоугольников и трапеций равны h3 4 f x) и h3 1 f x) и имеют противоположные знаки. Это означает, что точное значение интеграла лежит в вилке между ними. Объединяя 6) и 7), можно записать напомним, что здесь h = b ): I = + h3 4 f x)+ h5 19 fiv x)+, I = T h3 1 f x) h5 48 fiv x)+. Умножим первое равенство на /3, а второе на 1/3 и сложим: I = T h5 88 fiv x)+ = S h5 88 fiv x)+. Мы видим, что новая формула S, называемая формулой Симпсона или формулой парабол), имеет вид S = T = h ) ) +b f)+4f +fb), h = b. 8) 6 1
11 Главный член погрешности этой формулы равен h5 88 fiv x). Отметим важное обстоятельство: из самого вывода этой формулы видно, что ее остаточный член R S на равномерной сетке имеет разложение по нечетным степеням h = b, начиная с h 5. Следовательно, формула Симпсона точна для многочленов третьей степени. Коэффициенты формулы Симпсона и оценку остаточного члена можно вывести из ) при n = 3 для сетки замкнутого типа. Из этого построения следует, что полученная формула должна быть точна для многочленов второй степени, однако она обладает повышенной точностью в силу своей симметрии. Изменим запись формулы Симпсона 8), рассматривая ее как частный случай формулы Ньютона Котеса для равномерной сетки из трех узлов с шагом h = b )/: S = h ) +b f)+4f 3 ) +fb) Составные формулы трапеций, средних прямоугольников и Симпсона В общем случае длина отрезка [, b] не мала, а потому остаточный член в рассматриваемых формулах может быть велик. Для повышения точности на отрезке интегрирования вводят достаточно густую сетку = x 1 12 составная формула средних прямоугольников + 1 n 1 4 I = n 1 i + n 1 R h 3 if xi +x i+1 составная формула Симпсона = i n 1 ) xi +x i+1 h i f + ) n 1 ) h 5 if IV xi +x i+1 + ; n 1 n 1 I = S i + Ri S = 1 n 1 ) xi +x i+1 h i f i +4f )+f i n 1 ) h 5 88 if IV xi +x i+1 +. На равномерной сетке остаточные члены этих квадратурных формул могут быть представлены следующим образом отбрасываем члены, содержащие более высокие степени h): R T h n 1 ) hf xi +x b i+1 h f x)dx; 1 1 R h n 1 ) hf xi +x b i+1 h f x)dx; 4 4 R S h4 n 1 ) hf IV xi +x b i+1 h4 f IV x)dx Приведенные оценки являются асимптотическими, т.е. выполняются при h с точностью до членов более высокого порядка малости. Но для справедливости этих оценок необходимо существование непрерывных производных подынтегральной функции соответствующих порядков. Если эти производные кусочно-непрерывны, то можно сделать только мажорантные оценки: R T b ) 1 b ) h M, R 4 M = mx [,b] h M, R S b ) 88 h4 M 4, f x), M 4 = mx f IV x). [,b] 1
13 1.4. Построение формул Ньютона Котеса методом неопределенных коэффициентов Для повышения точности интегрирования можно увеличивать число точек n, по которым подынтегральная функция аппроксимируется интерполяционным многочленом Лагранжа. В этом случае коэффициенты c i и остаточный член R определяются из ). Более наглядным способом мы построили формулы замкнутого типа для n = и n = 3 и формулу открытого типа для n = 1. Рассмотрим другой способ построения квадратурных формул, называемый методом неопределенных коэффициентов. Будем рассматривать только формулы замкнутого типа. Пусть n = 4. Формула не будет симметричной, т.к. середина отрезка не входит в состав сетки. Построим формулу вида где I = b fx)dx c 1 fx 1 )+c fx )+c 3 fx 3 )+c 4 fx 4 ), x 1 =, x = + b 3, x 3 = + b 3, x 4 = b. Коэффициенты c i будем подбирать так, чтобы эта формула была точна для многочленов как можно более высокой степени. Для упрощения выкладок проведем стандартную замену переменных x = b+ + b в результате которой получим b fx)dx = b t, 1 t 1, x b, 1 f xt) ) dt = b 1 gt) dt. Теперь, с учетом проведенной замены, исходная задача свелась к поиску таких коэффициентов d i, чтобы квадратурная формула 1 gt)dt d 1 g 1)+d g 1/3)+d 3 g1/3)+d 4 g1) была точна для многочленов наиболее высокой степени. Погрешность квадратуры имеет вид Rg) = gt)dt d 1 g 1)+d g 1/3)+d 3 g1/3)+d 4 g1) ). 1 13
14 m Подставим в это выражение многочлен gt) = j t j степени m: j= Rg) = n j R t j). j= Подбором d i мы попытаемся добиться выполнения равенств R1) =, Rt) =. Rt m ) = при возможно большем значении m. Подставим в эти равенства последовательно gt) 1, gt) t, gt) t, gt) t 3. Получим систему из четырех линейных уравнений, решение которой дает следующие значения коэффициентов квадратурной формулы: d 1 = d 4 = 1 4, d = d 3 = 3 4. Проверкой убеждаемся, что построенная квадратурная формула не будет точна для gt) t 4. Учитывая, что c i = d i b )/ и b = 3h, получим искомую квадратурную формулу с введенными выше узлами x i : b fx)dx 3 8 hf 1 +3f +3f 3 +f 4 ). 9) Эту квадратуру называют правилом или формулой) трех восьмых. Поскольку она точна для многочленов третьей степени, то ее остаточный член имеет порядок Oh 5 ) т.е. тот же, что и остаточный член формулы Симпсона), однако раскладывается в ряд Тейлора по последовательным степеням h. Ту же формулу можно получить из ) при n = 4. Погрешность формулы 9) представляется в виде 3h 5 f IV ξ)/8, ξ [,b]. Составная формула трех восьмых может быть построена разбиением [, b] на элементарные подотрезки с последующим применением 9) на каждом из них. Мы же построим ее в несколько другом виде. Возьмем число элементарных подотрезков кратным трем, т.е. равным 3m, m = 1. Тогда число узлов полученной сетки n = 3m + 1, а ее шаг h = b )/n 1) = b )/3m. Возьмем строенный отрезок [+h, + +3)h], =,3,6. и применим на нем правило 9): ++3)h fx)dx = 3 8 hf 1 +3f +3f 3 +f 4 )+R +1, где +h R +1 = 3h5 8 fiv ξ +1 ), ξ +1 [+h, + +3)h]. 14
15 Если эти равенства расписать для всех остальных строенных подотрезков и сложить их почленно, то получим составную формулу трех восьмых: b fx)dx = 3h 8 [ f1 +f n )+f 4 +f 7 + +f n 3 )+ +3f +f 3 +f 5 +f 6 + +f n +f n 1 ) ] +R, где погрешность R ведет себя следующим образом: R = h4 [ 3h f IV ξ 1 )+ +f IV ξ m/3 ) ]) b h4 f IV x)dx. 8 8 Мажорантная оценка имеет вид R b ) 8 h 4 M 4, M 4 = mx f IV x). [,b] Хотя формула Симпсона по точности и количеству вычислений значений подынтегральной функции предпочтительнее, правило трех восьмых имеет самостоятельное значение, т.к. оно может быть использовано в случае таблично заданной функции и нечетного числа элементарных отрезков, когда формула Симпсона неприменима. Возьмем теперь n = 5, h = b )/4. В этом случае квадратурная формула будет симметрична, т.к. середина отрезка входит в состав сетки; следовательно, она будет точна не только для многочленов четвертой степени, но и пятой. Применяя описанный выше метод неопределенных коэффициентов, получим следующую формулу ее иногда называют формулой Боде): b fx)dx b ) 7 45 f f f f ) 45 f 5 = 14 =h 45 f f f f ) 45 f 5. Главный член ее погрешности равен h7 f VI ξ), ξ [,b], а остаточный член разлагается по нечетным степеням h, начиная с h 7. Может быть построена составная формула Боде с главным членом погрешности Oh 6 ). 15
16 1.5. Квадратурные формулы Гаусса При построении квадратурных формул Гаусса важную роль играют ортогональные многочлены Лежандра, имеющие вид: P n x) = 1 n n! d n dx n [ x 1) n], n =,1. Эти многочлены обладают следующими свойствами: 1) P n 1) = 1, P n 1) = 1) n, n =,1. ; ) 17 Для того чтобы квадратурная формула была точна для многочленов степени до n 1 включительно, необходимо и достаточно, чтобы она была точна для одночленов ft) 1, t, t. t n 1. Будем подбирать параметры t i и c i методом неопределенных коэффициентов, который мы использовали ранее для построения формул Ньютона Котеса. Получим систему n уравнений с n неизвестными: n c i =, n c i t i =. n c i t n i = n c i t n 1 i =. n 1, 1) Эта система нелинейная, поэтому ее решение затруднено уже при небольших значениях n. Кроме того, нет никакой гарантии, что она вообще разрешима, или что ее решенияt i вещественны и принадлежат отрезку[ 1,1]. Докажем разрешимость системы 1). Рассмотрим семейство многочленов вида f t) = t P n t), =,1. n 1, где P n t) многочлен Лежандра. Так как степени многочленов f t) не превосходят n 1, то, если t i и c i удовлетворяют системе 1), должны быть выполнены равенства 1 t P n t)dt = n c i t i P nt i ), =,1. n 1. В силу свойства ортогональности многочленов Лежандра справедливы равенства t P n t)dt = 18 Эти равенства заведомо будут выполняться при любых c i, если положить P n t i ) =, i = 1. n, т.е. если в качестве узлов t i взять корни многочлена Лежандра P n t). Напомним, что эти корни вещественны, различны и симметрично относительно нуля расположены на 1, 1). Если мы подставим t i в 1), то из первых n линейных уравнений найдем коэффициенты c i. Они определяются при этом однозначно, поскольку определителем соответствующей линейной системы является определитель Вандермонда. Осталось показать, что построенная квадратурная формула точна для всех многочленов степени до n 1 включительно. Действительно, представим произвольный многочлен Q n 1 t) степени n 1 в виде Q n 1 t) = G n 1 t)p n t)+f n 1 t), где G n 1 и F n 1 многочлены степени n 1. Предложенная квадратурная формула по своему построению точна для каждого из слагаемых в этом разложении, следовательно, точна и для произвольного многочлена степени не выше n 1. В случае применения квадратуры Гаусса для вычисления интеграла на произвольном отрезке [, b] делаем замену переменных в результате которой получим b x = b+ fx)dx = b 1 + b b+ f t, + b ) t dt. Последний интеграл заменим квадратурной формулой Гаусса: b fx)dx b n c i fx i ), где x i = b+ + b t i, i = 1. n, а t i корни многочлена Лежандра P n t) на [ 1, 1]. Остаточный член формулы Гаусса с n узлами имеет вид R n = b )n+1 n!) 4 f n) ξ) [n)!] 3 n+1) 18
19 с учетом формулы Стирлинга) В частности, b.5 n ) n b f n) ξ), ξ [,b]. 3n R R ) 5 b f 4) ξ), ) 7 b f 6) ξ). Кроме того, имеем t i = t n i+1 в силу свойств корней многочлена Лежандра) и c i = c n i+1 в силу равенства коэффициентов при симметричных узлах). Могут быть построены составные формулы Гаусса, если отрезок интегрирования разбить на подотрезки и на каждом из них применить формулу Гаусса. Однако в отличие от составных формул Ньютона Котеса замкнутого типа, в которых можно экономить одно обращение к подынтегральной функции на каждой граничной точке подотрезка кроме концов всего отрезка интегрирования), в составных формулах Гаусса такой возможности нет Правило Рунге практической оценки погрешности квадратурных формул Из п следует, что чем выше степень многочлена Лагранжа L n x), аппроксимирующего заданную подынтегральную функцию, тем выше точность соответствующей квадратуры. Это верно для достаточно малых длин отрезков интегрирования. Кроме того, не следует забывать, что для аппроксимации, как правило, используют L n x) для n, не б ольших4или5. Поэтому применяется аппроксимация не сразу на всем отрезке, а кусочнополиномиальная аппроксимация многочленами невысокой степени, что приводит к построению составных квадратур. При этом возможны два способа: либо берут достаточно густую равномерную сетку и на ней строят составную формулу; либо разбивают отрезок на подотрезки и затем на каждом из них применяют элементарную квадратуру с одновременным суммированием полученных приближений интеграла на этих подотрезках. Шаг сетки и длины подотрезков следует подбирать из того требования, чтобы значение интеграла вычислялось с некоторой заранее предписанной точностью. Для оценки достигнутой при вычислении интеграла точности мы имеем пока только асимптотические или мажорантные оценки остаточных членов квадратурных формул, которые либо трудно использовать на 19
20 практике, либо дают излишне малые значения для шага интегрирования и длины подотрезков. Поэтому применяется другой практически эффективный и легко реализуемый способ контроля точности интегрирования, основанный на оценке главного члена погрешности квадратуры при помощи правила Рунге. Этот метод называют еще методом двойного пересчета или экстраполяцией по Ричардсону. Пусть на [, b] взята равномерная сетка
21 Выписанное представление интеграла носит название второй формулы Рунге. Таким образом, мы описали еще один способ построения квадратурных формул, который может быть использован наряду с методом неопределенных коэффициентов и аппарата интерполирования. Поскольку можно ожидать, что приближение I h/ более точно, чем I h, то разумно в качестве нового приближения брать I = I h/ + I h/ I h +c p h p+1 +O h p+), 14) 1 > <<>I h,h/ которое получается, если в 1) подставить оценку ) p h c I h/ I h p 1 Повторив процесс с вдвое меньшим шагом, получим формулу с главным членом погрешности порядка h p+ и т.д. Таким образом, мы получили не только способ контроля точности при вычислении интеграла, дающий возможность выбора надлежащей величины шага h, но и простой метод построения все более точных квадратурных формул без явного выписывания их в виде квадратурной суммы с коэффициентами Котеса). В качестве исходной формулы можно взять простую формулу невысокой точности например, трапеций или средних прямоугольников). Покажем, что операция построения формулы I h,h/ носящей имя Ричардсона) есть операция экстраполирования. Для этого надо показать, что, если I h I h/, то I h,h/ всегда лежит вне отрезка с границами I h и I h/. Действительно, если I h/ > I h, то т.к. p 1. Если I h/ I h/, I h,h/ = I h/ I h I h/ p 1 22 увеличивается на два, а не на единицу, как это имеет место в случае несимметричных формул. Наиболее употребительным является случай, когда в качестве первоначальной формулы выбирается формула трапеций T, для которой, как было показано выше, имеет место следующее разложение остаточного члена если подынтегральная функция имеет s + непрерывные производные на [,b]): I = T h +c 1 h +c h 4 +c 3 h 6 + +c s h s +O h s+). Для целей, которые будут видны из последующего изложения, переобозначим T h T h, и T h/ T h/,. Тогда по второму правилу Рунге исключим первый член c 1 h/) разложения в ряд Тейлора остаточного члена у T h/, и получим новую формулу T h/,1 = T h/, + T h/, T h,, 1 у которой остаточный член раскладывается также по четным степеням, но начиная с h 4 : I = T h/,1 +b h 4 +b 3 h 6 + +b s h s +O h s+). Новая формулаt h/,1 имеет порядок точности, превышающий на два порядок точности формулы трапеций, т.е. она точна для многочленов третьей степени. Покажем, что новая формула T h/,1 есть формула Симпсона, записанная в неявном виде, т.е. без непосредственного выписывания коэффициентов Котеса. Для простоты вывода будем полагать, что h = b и формула трапеций строится по двум узлам. Тогда T h, = h f 1 +f ), T h/, = h 4 f 1 +f 1/ +f ), T h/,1 = h 4 f 1 +f 1/ +f )+ + 1 h 3 4 f 1 +f 1/ +f ) h ) f 1 +f ) = = h 6 f 1 +4f 1/ +f ) = S. Получили формулу Симпсона, выписанную по сетке с тремя узлами. Теперь мы можем применить правило 14) уже к формуле Симпсона T h/,1, использовав уже вычисленное T h/4,1 : T h/4, = T h/4,1 + T h/4,1 T h/,1. 4 1
23 Здесь в знаменателе двойка возводится в четвертую степень, т.к. при вычислении T h/4, мы исключили член b h 4 в разложении по степеням h остаточного члена формулы T h/,1. В результате для новой формулы T h/4, остаточный член разлагается по четным степеням h, начиная с h 6 : I = T h/4, + 3 h s h s +O h s+). Следовательно, алгебраический порядок точности формулы T h/4, превышает порядок точности формулы Симпсона на два и равен пяти т.е. она точна для многочленов пятой степени). Можно показать, что она совпадает с ранее выписанной формулой Боде, которая может быть получена методом неопределенных коэффициентов на сетке с пятью узлами. Если мы применим правило Рунге к формуле Боде T h/4,, то получим формулу по девяти узлам, имеющую седьмой порядок точности; в разложении остаточного члена этой формулы присутствуют четные степени h, начиная с восьмой. С каждым следующим применением правила Рунге к вновь полученным формулам мы будем получать формулы с порядком точности 9, 11 и т.д. Отметим, что, начиная с формулы, построенной по 14) на основе формулы Боде, порядок точности будет меньше, чем если бы эти формулы строились в качестве формул Ньютона Котеса, т.е. как формулы интерполяционного типа. Например, для n = 9 формула Ньютона Котеса имеет девятый порядок точности, а не седьмой, как формула, полученная из формулы Боде по 14). Поэтому начиная с этого момента по второму правилу Рунге нельзя построить формулы Ньютона Котеса это будут уже другие формулы, не являющиеся точными для многочленов наиболее высокой степени Процесс Эйткена оценки фактической точности квадратурной формулы У всех рассмотренных квадратурных формул остаточный член можно разложить в ряд по степеням шага сетки. Следовательно, к ним применим метод Рунге при этом требуется, чтобы был заранее известен порядок точности исходной формулы). Кроме того, разложение остаточного члена проводится при условии достаточной гладкости подынтегральной функции. Если же функция не обладает требуемой гладкостью, то это разложение уже не имеет смысла, т.е. нам уже не известен фактический реальный) порядок точности формулы. Рассмотрим процесс Эйткена, которой на основе повторных просчетов на нескольких сетках дает возможность: 3
24 определить фактический порядок точности p квадратурной формулы для заданной подынтегральной функции; уточнить результат, полученный на исходной сетке. В упрощенном виде процесс Эйткена можно описать следующим образом. Выберем три сетки с шагами h 1 = h, h = h/, h 3 = h/4. Вычислим приближения I h, I h/ и I h/4 к интегралу I по выбранной квадратурной формуле. Тогда, если учитывать только главный член погрешности, имеем три уравнения для определения I, и p, где p фактический, заранее неизвестный порядок точности формулы для данной подынтегральной функции: I = I h +ch p, I = I h/ + 1 p chp, I = I h/ p chp, в которой значения I, c и p неизвестны. Из первого и второго уравнений имеем ch p 1 1 ) = I p h/ I h. Из второго и третьего уравнений получим ) p chp = I p h/4 I h/. Из последних двух равенств получаем уравнение для определения p: p = I h/ I h I h/4 I h/. Оценка для главного члена погрешности имеет вид ch p = Ih/ I h ) I h/ I h I h/4. Полученную оценку можно использовать для уточнения I h. Пример 1. Пусть вычисляется следующий интеграл методом трапеций: xdx. Тогда для этого интеграла фактический порядок точности p будет равен 1/. 4
25 Пример. Пусть вычисляется следующий интеграл методом трапеций: n xdx. Тогда для этого интеграла фактический порядок точности p будет равен n 1)/n Способы построения практических алгоритмов численного интегрирования. Адаптивные алгоритмы Ранее мы рассмотрели различные квадратурные формулы и способ Рунге контроля достигнутой при их применении точности, однако оставили в стороне вопрос построения практически пригодных алгоритмов. Всегда можно подобрать такую подынтегральную функцию fx), для которой любая заранее выбранная программа численного интегрирования дает совершенно неверные результаты. Сказанное иллюстрирует следующий простейший пример. Возьмем fx) 1 и вычислим интеграл по данной программе. Запомним выбранные этой программой узлы интегрирования x i и доопределим функцию fx) следующим образом: fx i ) = 1, fx) = A 1, x x i. Очевидно, что при повторном запуске для видоизмененной функции программа даст прежний результат, хотя верный результат может значительно отличаться от вычисленного при надлежащем задании A. Таким образом, при построении практических алгоритмов класс подобных примеров должен быть максимально сужен, насколько это возможно без неоправданного усложнения логики или потери эффективности на распространенных типах задач. Рассмотрим несколько алгоритмов такого рода, в которых для контроля достижения предписанной точности применяется автоматический выбор шага интегрирования. При этом будем учитывать, что б ольшая часть машинного времени в реальных расчетах тратится на вычисление подынтегральной функции. Поэтому программа считается тем эффективней, чем меньше таких вычислений она производит для достижения заданной точности Простейший неадаптивный) алгоритм Зададим на отрезке интегрирования n узлов и вычислим по выбранной квадратурной формуле приближенное значение I n интеграла I по этим узлам. Затем увеличим число узлов вдвое и вычислим I n. По первому правилу Рунге вычислим оценку погрешности ошибки, которая, конечно, не учитывает погрешность суммирования R Σ, а дает лишь 5
26 оценку главного члена погрешности самой формулы, т.е. метода интегрирования: δ = I n I n p 1, где p порядок точности используемой формулы. Если δ ε, то мы опять удваиваем текущее значение n и повторяем описанный процесс. Данный алгоритм не учитывает локальный характер поведения подынтегральной функции fx). Общее количество узлов равномерной сетки подбирается с учетом наихудшего поведения fx) на отдельных участках которые могут иметь небольшую суммарную длину), при этом не учитываются участки, где fx) меняется медленно. Такие алгоритмы называют неадаптивными. Рассмотренный алгоритм удобен при применении формул Ньютона Котеса, для которых нетрудно реализовать составные формулы. Для формул Гаусса этот алгоритм неудобен, т.к. надо хранить в программе большое число узлов и коэффициентов из-за того, что заранее неясно, какое n окажется пригодным Модификация простейшего алгоритма Отрезок[, b] разбивается на n подотрезков частей), а затем на каждом подотрезке применяется выбранная квадратурная формула. В процессе вычислений приближенные значения интегралов по подотрезкам суммируются, что дает в результате приближение I n по всему отрезку. Затем n удваивается, вычисляется I n, после чего поступают как в алгоритме Данная модификация обладает теми же недостатками, но более удобна для применения формул Гаусса, т.к. требует хранения узлов и коэффициентов только одной выбранной формулы. Контроль точности разумно вести по критерию: δ 27 использованы при вычислениях с n узлами или на n подотрезках). При применении формул Гаусса такой возможности нет, поскольку в обоих случаях при повторном счете узлы будут совершенно другими. Поэтому применение правила Рунге для получения оценки погрешности при использовании формул Гаусса является неэффективным. Вместо правила Рунге для формул Гаусса могут быть использованы различные подходы. Простейший из них следующий. Выполним вычисления по двум формулам Гаусса с n и n + 1 узлами. На это требуется n + 1 вычисление fx). Разность между ними примем за оценку погрешности, допущенную формулой с n узлами. Однако формула с n + 1 узлами ненамного точнее формулы с n узлами, поэтому вычисление оценки погрешности примерно столь же трудоемко, как и вычисление самого интеграла. Поэтому более эффективен другой подход. Возьмем формулу Гаусса с n узлами и построим другую формулу она называется формулой Гаусса Кронрода) с теми же n узлами и с n+1 дополнительными узлами: b fx)dx = n c i fx i )+ n d j fξ j )+R n+1 f). Узлы ξ j и коэффициенты c i, d j подбираются так, чтобы данная формула имела полиномиальный порядок точности 3n + 1. Выполним оба просчета по формуле Гаусса с n узлами и по новой формуле с n+1 узлами, затратив на это n + 1 вычислений значений подынтегральной функции. Разность между этими просчетами принимается за оценку погрешности. Таким образом при просчетах по формулам Гаусса с n и n+1 узлами и по формуле Гаусса Кронрода с n и n + 1 узлами требуется одинаковое количество вычислений fx) в n+1 точках, однако результат, полученный по формуле Гаусса Кронрода, гораздо более точен. Пользуясь свойствами формул Гаусса Кронрода, алгоритм можно представить следующим образом. Возьмем формулу Гаусса с тремя узлами, добавим четыре узла по Кронроду, а затем выполним вычисления по обеим формулам. Если их разность, являющаяся оценкой погрешности, больше предписанной точности, то добавляем еще 8 узлов всего получается 15), вычисляем приближенное значение интеграла по новой формуле Гаусса Кронрода и сравниваем полученный результат с результатом по предыдущей формуле. Если точность не достигнута, то добавляем еще 16 узлов всего получается 31 узел) и т.д. Данное видоизменение алгоритма обеспечивает механизм получения все более точных результатов без потери ранее вычисленных значений подынтегральных функций. Если заданная точность так и не оказалась достигнутой, то исходный отрезок [, b] можно разбить на два подотрезка и на каждом из них применить описанный алгоритм. 7 j=1
28 Адаптивный алгоритм последовательного передвижения слева направо Рассмотрим адаптивный алгоритм т.е. алгоритм, учитывающий при контроле точности локальное поведение подынтегральной функции), который широко применяется при интегрировании обыкновенных дифференциальных уравнений. Зададимся некоторым шагом h b и вычислим приближенное значение I h интеграла I по выбранной элементарной квадратурной формуле на отрезке [, +h]. Затем вычислим по той же формуле два приближения на отрезках [, +h/] и [+h/, +h] и сложим их. Полученную сумму обозначим I h/. Вычислим локальную погрешность на шаге h по первой формуле Рунге: δ h = I h/ I h p 1. Если δ h > mxε A, ε O I h/ ), то положим h = h/ и повторим вычисления. В противном случае считаем, что локальная погрешность обеспечивает заданную точность и текущий шаг, который обозначим через h 1, признается удачным. Теперь возникает вопрос: каким должен быть следующий шаг? Здесь возможны два случая: либо следующий шаг выбирается равным предыдущему, либо, если текущая локальная погрешность очень мала например, мы находимся в области гладкости подынтегральной функции), то следующий шаг следует выбрать б ольшим, чем предыдущий. Часто выбирают следующую стратегию. Если δ h > 1 p mx ε A, ε O I h/ ), то следующий шаг h полагается равным h 1, в противном случае берем h = h 1. Такое удвоение обосновывается возможным прогнозом о гладкости подынтегральной функции на следующем шаге. Показатель p выбирается из таких соображений. Для следующего шага h имеет место соотношение h ) p, I I h / c причем на предыдущем шаге h 1 было выполнено h1 ) p 1 c 29 мало изменится из-за гладкости fx), возьмем h = h 1. Тогда возможно соотношение h ) p c = c h1 ) p mxε A, ε O I h1/ ), т.е. предположительно с некоторой вероятностью) шаг h = h 1 окажется удачным. Выбрав следующий шаг, проверяем, не выйдет ли он за правый предел отрезка [, b]. Если нет, то повторяем процесс, передвигаясь слева направо, если да, то предварительно меняем следующий шаг, полагая его равным разности между b и достигнутой точкой отрезка. Итак, строится последовательность шагов h i, значений I hi /, которые суммируются по достижении конца отрезка их сумма дает приближенное значение всего интеграла), и локальных оценокδ i достигнутой точности на шаге h i. Теперь мы видим, что там, где подынтегральная функция резко меняется, значения h i уменьшаются, а на гладких участках увеличиваются, и лишние вычисления не выполняются. Таким образом, алгоритм как бы адаптируется к характеру поведения fx). Критерием достигнутой точности мы выбираем i δ i 30 Поэтому при использовании критерия 15) значение E может оказаться излишне меньшим ε, т.е. результат будет получен с завышенной точностью. Вследствие этого часто использую более слабый критерий, задавая для локальной ошибки значение ε примерно равным ε/b ), где ε требуемая точность для всего отрезка интегрирования. При этом достигается экономия вычислительной работы. Для формул Ньютона Котеса замкнутого типа на равномерной сетке дополнительная экономия достигается также, если запоминать значения подынтегральной функции в общих, уже использованных узлах. Например, если в качестве расчетной формулы была выбрана формула Симпсона, то на шаге h i она она использует три узла, а на двух шагах h i / пять узлов, причем три из них уже использованы. Следовательно, во втором применении формулы Симпсона требуется лишь два новых значения fx). Ясно, что при передвижении вправо с шагом h i+1 должны использоваться значения fx), уже вычисленные в правом конце i-го подотрезка. Кроме того, при каждом делении h i пополам следует запоминать узлы и значения fx) из правой половины для последующего использования. Поскольку надо заранее учитывать максимальный объем памяти ЭВМ для хранения этих величин, то задается предел HALFMAX на уровень деления отрезка пополам. Когда этот максимум достигнут, то текущий подотрезок считается в любом случае приемлемым, даже если локальная ошибка на нем слишком велика. Длина такого подотрезка может быть сделана совсем небольшой, если выбрать HALFMAX достаточно большим: например, если HALFMAX=3, то длина подотрезка равна b )/ 3. Отметим возникающее при этом важное обстоятельство. Появление такого рода подотрезка заданной минимальной длины означает, что fx) имеет внутри него некоторую особенность например, неинтегрируемую особенность). Поэтому необходимо сообщить пользователю границы подотрезка для последующего осмысления результата. Поскольку длина такого подотрезка мала, то можно надеяться, что локальная ошибка на нем все же не очень большая и при дальнейшем интегрировании она может быть компенсирована другими локальными ошибками, имеющими противоположные знаки. Поэтому процесс разумно продолжить, а не прерывать расчеты. При таком подходе мы имеем возможность не только получить в конечном счете значение интеграла с глобальной ошибкой, меньшей предписанного ε, но и выявить все точки особенностей fx), препятствующие получению локальных погрешностей с приемлемой точностью. Такой подход к составлению программы численного интегрирования позволяет всегда получить приближенное значение интеграла; при этом, если предписанная точность не достигнута, то выдается наилучший возможный для нее результат, оценка E реально достигнутой точности и, если это имеет место, точки особенностей неинтегрируемости) fx). Если в случае формул Ньютона Котеса замкнутого типа удается до- 3
31 стигнуть максимальной экономии вычислений fx), то для формул Гаусса этого сделать нельзя, поскольку узлы при их применении на всем подотрезке и на двух его половинах будут совершенно различными. Поэтому применение правила Рунге для оценки локальной погрешности для формул Гаусса будет неэффективно с точки зрения производимых затрат. Для уменьшения трудоемкости применяют формулы Гаусса Кронрода либо выбирают формулу Гаусса с небольшим числом узлов например, n = 5). В последнем случае поступают так: на текущем подотрезке вычисляют значения интеграла по пяти и шести узлам, а в качестве оценки локальной ошибки берут разность между ними. При этом вследствие высокого порядка точности формул Гаусса можно надеяться, что продвижение по [, b] слева направо будет идти с относительно большими шагами, что позволит снизить общие трудозатраты, обусловленные двумя независимыми расчетами при n = 5 и n = 6 на каждом подотрезке для оценки локальной ошибки. Дополнительное улучшение может быть сделано из практически подтвержденного наблюдения, что для многих подотрезков достигается б ольшая точность, чем это требуется, тогда как для некоторых подотрезков заданная граница погрешности превышается ненамного. Поэтому может быть использован прием, получивший название банкирование : мы вкладываем в банк разность между предписанной точностью и оценкой локальной погрешности для приемлемого подотрезка с тем, чтобы использовать ее впоследствии для уменьшения точности на других интервалах, на которых предписанная точность не достигается. Это может быть сделано, например, так: если i δ i > mx ε A, ε O то проверяется менее жесткое условие j=1 j=1 ) I hj /, i ) δ i 32 начало α = ; β = b; I = ; E = ; = ; mx = ; пока α mx = ; печать точка особенности α, β)); I = I +I hi / +δ i ; E = E +δ i ; α = β; β = b; печатьi, E); конец Эта модификация выписана для формул Ньютона Котеса с контролем локальной точности по правилу Рунге. В случае формул Гаусса можно положить I hi = G 5, I hi / = G 6, где G 5 и G 6 результаты вычисления интеграла на шаге h i по формулам Гаусса с 5 и 6 узлами соответственно и δ i = G 6 G 5, либо применить формулы Гаусса Кронрода Адаптивный алгоритм с контролем точности по глобальной ошибке Рассмотренный алгоритм при движении слева направо реализует локальный способ контроля точности, поскольку вывод о том, удовлетворяет ли полученная локальная оценка интеграла на каждом подотрезке предписанной точности или нет, делается независимо от локальных оценок на других подотрезках. Это означает, что порядок, в котором выбираются подотрезки, не играет никакой роли для получения конечного результата. Таким образом, выбор очередного подотрезка для разбиения следует делать целиком из соображений легкости программной реализации, т.к. он не отражается на эффективности алгоритма. Алгоритм использует принцип слева направо, когда разбивается самый левый подотрезок, для которого локальная ошибка велика. Такой принцип нетруден для программирования и при относительно небольших затратах памяти дает возможность запоминать, а затем использовать ранее вычисленные значения подынтегральной функции. Однако можно предложить алгоритм, в котором проверка локальных ошибок не проводится вовсе. Вместо этого контроль точности осуществляется по глобальной погрешности. Выпишем схему алгоритма: 3
http://lib.qrz.ru/node/8921
http://docplayer.com/114131014-I-o-arushanyan-primenenie-metoda-kvadratur-dlya-chislennogo-resheniya-integralnyh-uravneniy-vtorogo-roda-uchebnoe-posobie.html