Решение интегральных уравнений Вольтерра неявным методом Рунге-Кутты тема автореферата и диссертации по математике, 01.01.07 ВАК РФ

Российская Академия наук Сибирское отделение Институт вычислительной математики и математической геофизики

РГ6 ОД На правах рукописи

УДК 518.5 51-2.83

1 )£В К-

Савченко Александр Оливерович

Решение интегральных уравнений Вольтерра неявным методом Рунге-Кутты

01.01.07 - вычислительная математика

Автореферат диссертации на соискание ученой степени кандидата, физико-математических наук

Новосибирск, 1997

Работа выполнена в Институте вычислительной математики и математической геофизики СО РАН.

Научные руководители: д.ф.-м.н. Кузнецов Ю.И.,

д.ф.-м.н. Цецохо В.А.

Официальные оппоненты: д.ф.-м.н. Григорьев Ю.Н.

к.ф.-м.н. Левыкин А.И.

Ведущая организация - Институт Математики СО РАН.

Защита состоится " // года в " " час.

на заседании специализированного совета К002.10.01 по присуждению ученой степени кандидата физико-математических наук в Институте вычислительной математики и математической геофизики по адресу: 630090, Новосибирск 90, просп. Академика Лаврентьева, 6.

С диссертацией можно ознакомиться в читальном зале библиотеки ИВМ и МГ (пр. Академика Лаврентьева, 6).

Автореферат разослан " " Зк^а^Лг 199<| г.

Ученый секретарь специализированного совета д.ф.-м.н. Ю.И. Кузнецов

Общая характеристика работы

Актуальность темы. Интегральные уравнения Вольтерра находят широкое прменение в задачах механики сплошной среды, электродинамики, астрономии, экологии, сейсмики и т.д. Это объясняет интерес, который имеется к их теории и методам решения. Поскольку точное аналитическое решение данного класса уравнений возможно лишь в некоторых частных случаях, особую актуальность приобретают численные методы их решения.

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

В работах Де Хуга и Вейсса (1973-75 гг.), впервые предложен аналог неявного РК - метода для решения интегральных уравнений Вольтерра I и II рода, доказана устойчивость и сходимость предложенного метода. При этом элементы квадратурной матрицы находились как интегралы от коэффициентов многочлена Лагранжа для наперед заданного набора узлов разбиения шага интегрирования. Однако в реальных расчетах условия сходимости численного решения к искомому являются неконструктивными, т.к. шаг интегрирования в большинстве случаев либо фиксирован, либо может быть уменьшен до некоторой конечной величины. В данной диссертации предложен другой подход к построению матрицы квадратурных коэффи-

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

Цель работы. - Построение, анализ и апробация высокоэффективных численных методов решения интегральных уравнений Вольтерра на основе неявных методов типа Рунге-Кутты.

Научная новизна и практическая ценность. Все

результаты диссертации, начиная с параграфа 3 главы 1, являются новыми. Доказано, что при надлежащем выборе квадратурных коэффициентов и узлов, значение константы погрешности решения будет убывать экспоненциально с ростом числа стадий метода, что позволяет проводить интегрирование с большим временным шагом. Таким образом, высокая точность получаемых результатов достигается не столько за счет малости временного шага, сколько из-за малости константы в выражении для погрешности решения. Этим предложенный в диссертации подход значительно отличается от предшествующих.

Другая важная идея, реализованная в диссертации, связана с включением ядра, или его части, в квадратурные коэффициенты. Эта идея, применимая именно к решению интегрального уравнения, предлагается впервые. Разработанная на этой основе технология применена также к сингулярному интегральному уравнению - уравнению Абеля. Кроме того, дана модификация алгоритма для случая, обычно вызывающего затруднения при численных расчетах, когда ядро интегрального оператора имеет свойство <) = 0.

Апробация работы. Основные результаты работы докладывались на конференциях молодых ученых ВЦ СО РАН (1995, 1996, 1997 гг.), семинарах секции Вычислительной ма-

тематики (1996 г.). Полиостью диссертация докладывалась на секциях Физики атмосферы и океана (1997 г.) и Вычислительной математики (1997 г.) ВЦ СО РАН, на семинаре института Вычислительных технологий СО РАН (1997 г.).

Публикации. По результататам диссертации опубликовано 5 печатных работ и 3 работы находятся в печати.

Содержание работы

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

В первой главе приводится общая схема численного решения уравнений Вольтерра на основе неявного метода Рунге-Кутты. Значения квадратурных коэффициентов находятся из условия точного равенства интеграла и квадратурной формулы на классе полиномов, порядок котороых равен числу точек подсеточного разбиения. Найдена погрешность аппроксимации всех интегральных составляющих. Для уравнения Вольтерра II рода получен повышенный порядок аппроксимации в узлах сетки при дополнительном выборе квадратурных коэффициентов.

В первом параграфе этой главы рассматрена общая постановка задачи решения интегральных уравнений Вольтерра:

<4/(0 + /(0 = /' в, »(«))&, ¿о < г < Т. (1)

Jto

где а = 1 для уравнений Вольтерра II рода и а — 0 для уравнений Вольтерра I рода. Приведены достаточные условия на существование и единственность решения уравнения (1).

Во втором параграфе проведены дискретизация интервала интегрирования, выбор квадратур общего вида с неопределенными пока коэффициентами для аппроксимации элементарных интегралов и даны основные результаты, полученные другими авторами, применившими неявный метод Рунге-Кутты для решения интегральных уравнений. На отрезке [¿о, Т] вводится сетка с узлами ,

= г/г; г = 0, • • •, /; к =— подсеточное разбиение

= Ь + 0 < Аг < Л2 < • • • < А,„ < 1

и в предположении, что для для уравнений Вольтерра (1) выполняется:

Ут € [«о, Г] : |у(т)| < £, К[т, •, •) € х [~Ь,1]),

в точках ¿¿^ дан вид уравнения (1):

^ Г*к+1 /-г<>

(2)

Система уравнений, аппроксимирующая (2), выбирается в виде:

» — 1 Ш

2 /г^/фи, 1к1,Ук1) + £ ЬЬЯК (г,,, Ци У и) = + аУч

к-О 1=1 1=1

г = 0, •••,/- 1; j = 1,- • - ,пг ,

(3)

где Д,- = а Урд обозначает конечномерную аппроксима-

цию для у{1рч).

В третьем параграфе первой главы определяются элементы квадратурной матрицы В = {Ь^/), путем приравнивания

к нулю первых т слагаемых в погрешности аппроксимации интеграла на отрезке Погрешность аппроксимации р^

находится путем разложения подынтегрального выражения по формуле Тейлора. В этом случае

итп Ш»

где

г,- € [£,-, ¿,-+1], а сеточная функция <5^, } = 1, • • т имеет вид:

1 m

(4)

Система уравнений для определения матрицы В будет следующего вида:

BW = AWL, (5)

где Л = diag(\i, • • •, Am), W = (AjZ}) - матрица Вандермонда,

Для численного решения уравнения Вольтерра II рода узловые коэффициенты Aj выбираются равными узлам квадратурной формулы Гаусса, а уравнения Вольтерра I рода-узлам квадратурной формулы Радо с Am = 1.

Погрешность 8m = (¿im, ■ ■•, Smm) может быть выражена в явном виде как функция от узлов подсеточного разбиения Ау.

где M = 2m и M = 2m — 1 уравнений Вольтерра II и I рода соответственно.

В четвертом параграфе найдена погрешность приближения интеграла на отрезке [tb^fc+i]-

Квадратурные коэффициенты Ь/ выбираются равные квадратурным коэффициентам формулы Гаусса для решения уравнения Вольтерра II рода, и квадратурным коэффициентам формулы Радо для решения уравнения I рода. Отсюда погрешность аппроксимации тщь интеграла на отрезке ¿А.+1] квадратурной формулой будет иметь вид:

. при аппроксимации соответствующих интегралов в уравнении Вольтерра II рода, и

»,2т

при аппроксимации соответствующих интегралов в уравнении Вольтерра I рода; где гк £ а

*<»>-(«Г Р)

п М — 1т для квадратуры Гаусса, М = 2т—1 для квадратуры Радо.

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

Значения численного решения уравнения Вольтерра II рода в узлах сетки находятся путем введения дополнительной, т + 1-й точки подсеточного разбиения, совпадающей с правым концом соответствующего отрезка:

= ¿«+1» 1 = 0, •••,/— 1.

При дополнительных ограничениях на гладкость ядра исходного уравнения в узловых точках сетки £;, г = 0,•••,/:

из условия равенства нулю первых т слагаемых разложения погрешности аппроксимации интеграла на отрезке [t;,t|,m+i]) получено равенство нулю и последующих m — 1 слагаемых разложения. Отсюда погрешность аппроксимации на этом отрезке будет:

где сеточная функция Sjk определена в формуле (4).

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

В первом параграфе второй главы находится оценка погрешности решения где е,- = m). £ij = Yij -

y(tij), которая имеет следующий вид:

II ft ||<|| (Д - аЕ)~* II (1+ II (Д - аЕГ1 II • || А II Г1-•(P{m)iQ^+ || 6т || (QT+ || А || • || (Во - <хЕ)~1 || QS1))

где Е - единичная матрица, а элементы матрицы Д определяются как покомпонентное произведение элементов матрицы f?,-и производной ог ядра по последней переменной в соответствующих точках подсеточного разбиения. 6т = (5jm,-",5mm), P(m) определена в (7), а сеточные функции QfVQf, и элементы матрицы Di определяются по значениям производных от ядра, а также зависят и от шага Л сетки.

Во втором параграфе второй главы проводится численный анализ оценки погрешности решения £;.

Поскольку значения || (Я,— аЕ)'1 ||, || Я,-1|, и <2,^ в (8) зависят от ядра интегрального оператора, для них требуется только их ограниченность.

Значения двух оставшихся функций в формуле (8), а именно Р(т) и || &т || можно получить численно, предварительно найдя значения узлов Aj при ] = 1,---,т и вычислив и Р(т) по формулам (6) и (7) соответственно.

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

|| 6т ||и ехр(-1.5т - 0.7)

0.001 -10.00 ] -20.00 -30.00 4 -40.001 -50.00 -60.001

-70.00

0.00

5.00

10.00 15.00 20.00 25.00 т

Для функции Р(т) показатель экспоненты получается применением формулы Стирлинга:

Р(тп) = ЖТП02

—4т

где /? < еагр(т/3).

и

Таким образом, уменьшение третьего сомножителя в правой части неравенства (8) с ростом числа стадий тп носит экспоненциальный характер, и уменьшение всей нормы погрешности || £{ || должно быть аналогичным (за исключением, как показывают практические вычисления, достаточно редких случаев, когда норма матрицы Bf* будет иметь аналогичный характер экспоненциального роста с увеличением числа стадий).

В третьем параграфе находится оценка погрешности решения задачи с возмущениями ядра и правой части уравнения.

Рассмотрен случай, когда вместо K(t,s, y{s)) в формуле (1) задано /?(i, s, y(s)), а вместо J(t) - /(f), причем

||K(M,v(S))-K(M,yW)|| C([liiT] , [1ь,1Г1 < »,

II'<«>-/<«> II CfcTJ^

В предположении, что на функцию K(t,s,y(s)) наложены те же ограничения на гладкость, что и на исходную функцию Ii(t,s,y(s)), оценка погрешности || е; ||, г = 0, •••,/- 1, где ё« = • • ■,£im)T, iij = Yij - Uij, будет иметь вид:

II ii ||< qi(9m(i + 1 + gi) + + ffi)/A+ || Ri || +gi || po || ),

где элементы матрицы R{ состоят из сумм погрешностей аппроксимации в точках сетки:

1-1

Vij = Pij +

к=О

Ро = (роъ • • ■) Рот) - вектор погрешностей квадратуры на первом временном шаге, а коэффициенты д,- и <?,• определяются по элементам квадратурной матрицы В и значениям производных от ядра интегрального уравнения.

При малых возмущениях правой части и интегрального ядра при некотором фиксированном значении шага h, погрешность решения задачи с возмущениями будет близка, к оценке (8) решения задачи без возмущений.

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

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

В третьей главе предложен новый метод численного решения линейных уравнений Вольтерра на основе представления искомой функции по формуле Тейлора и использования значений ядра в точках подсеточного разбиения в качестве весовых коэффициентов квадратур. Особо рассмотрен случай с ядром интегрального оператора, для которого ¿) = 0. Предложенная технология применена также к сингулярному интегральному уравнению - уравнению Абеля.

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

Основная идея нового метода состоит в разложении функции у(£) (а не подынтегрального ядра, как ранее) по формуле Тейлора в точке Получив выражение для погрешности аппроксимации элементарного интеграла на отрезке [£¿,4^], значение матрицы квадратурных коэффициентов

находится, как и ранее, путем приравнивания к нулю первых т слагаемых этой погрешности. Отсюда получена следующая

система линейных уравнений:

Д(0 Ж = С(,) где С= (с^)> М7 - матрица Вандермонда и

к = 0, • ••, т - 1; ] = г = 0, ••*,/ — 1

Отсюда, с учетом (3), искомый вектор = (У^, • • •, У{т) будет находиться из решения системы уравнений:

(В^-аВД^-Х^У* > к=о

где

4 = (4^); 43 = м^а, «ы); я = (/¿1, • • •, )

- единичная матрица, а а = 1 для уравнения Вольтерра II рода и а = 0 для уравнения Вольтерра I рода.

Во втором параграфе рассмотрен метод решения уравнений Вольтерра первого рода с ядром вида к(Ь,Ь) = 0.

Для уравнений с ядром такого вида, диагональные элементы матрицы Врассмотренной в предыдущем параграфе, уже определены и равны нулю.

Для нахождения оставшихся т? — те элементов этой матрицы, введены дополнительные матрицы и вектрра где

Wj - матрица Вандермонда с пропуском ой строки и последнего столбца;

- прямоугольная матрица, полученная из С^ вычеркиванием последнего столбца;

Отсюда получена следующая система уравнений относительно

(О

неизвестных векторов иу:

И^ и^ = 4", з = \ (9)

где — У-я строка матрицы Б^У Найдя компоненты векторов и^р решением системы уравнений (9), искомые значения элементов матрицы В^ определяются из условий:

$ ~ при ^ = 1; з =

^¡1 = при ^ = ; +3 —

Ь^ = 0 при I = з\ з = 1, • • •, то

Порядок аппроксимации в этом случае будет на единицу меньше и равен тп — 1.

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

В четвертом параграфе рассмотрено применение метода, изложенного в первом параграфе этой главы, на примере решения уравнения Абеля:

{(П^-ЛО (0<«<1)

Ввиду того, что ядро этого уравнения имеет слабую сингулярность, очевидно, что применение аналога неявного метода Рунге-Кутты, которому посвящены две предыдущие главы, для ого решения будет неприемлимым.

Для нахождения элементов квадратурной матрицы, искомая функция у{0 разлагается по формуле Тейлора в точке

tjj, tij = t{ 4- Ajh. Найдя выражение для погрешности аппроксимации и приравняв первые т слагаемых ее к нулю, получим следующую систему уравнений для определения элементов квадратурной матрицы В = (Ьц): т

Eh(A,-A =

k = 1, ■ ■ ■, т; j = 1, m

При ее решении матрица Вандермонда будет меняться на каждом шаге j при нахождении j-он строки матрицы В.

Теоретические положения, развитые в предыдущих главах, подтверждены численным» экспериментами в четвертой главе. В модельных экспериментах были выбраны такие интегральные уравнения, решение которых было найдено другими численными методами. Применение изложенных в диссертации методов для решения этих уравнений уменьшило значение погрешности решения на несколько порядков по сравнению с другими методами. Проведен анализ полученных результатов. В частности, можно отметить, что несмотря на плохую обусловленность матриц Bj для уравнения Вольтерра I рода, особенно с ядром, для которого k(t,t) = 0, при увеличении числа стадий метода происходит уменьшение и погрешности решения до приемлемого уровня, что численно подтверждает характер изменения погрешности || £",• || в (8).

С другой стороны, для уравнения Вольтерра II рода матрица (В/ — Е) является хорошо обусловленной, так что норма ее обратной матрицы с увеличением числа стадий практически неизменна и небольшая по величине. По этой причине увеличение количества стадий метода в этом случае будет более существенным фактором уменьшения погрешности решения.

С увеличением числа шагов по времени также происходит уменьшение погрешности. Для уравнений Вольтерра I рода при этом происходит также и увеличение || BJ1 || и, тем

более, третьего сомножителя в правой части неравенства (8). Это обстоятельство несколько снижает эффект увеличения количества шагов по времени, однако полученные значение при этом тем не менее уменьшаются на несколько порядков.

Отметим преимущество использования неявного метода Рунге-Кутты высокого порядка точности особенно для решения уравнений Вольтерра I рода. Так, например, значение погрешности решения при применении 20-стадийного одноша-гового метода является меньшим, чем при использовании 10-стадийного четырехшагового метода.

Основные результаты диссертации

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

2. Предложен новый метод численного решения линейных уравнений Вольтерра на основе представления искомой функции по формуле Тейлора и использования значений ядра в точках подсеточного разбиения в качестве весовых коэффициентов квадратур. Разработанная на этой основе технология применена также к сингулярному интегральному уравнению - уравнению Абеля. Дана модификация алгоритма для случая, когда ядро интегрального оператора имеет свойство к(1,1) = 0.

3. Созданы программы численного решения интегральных

уравнений по всем указанным алгоритмам. Теоретические положения подтверждены численными экспериментами.

Публикации по теме диссертации

1. Савченко А.О. Решение интегральных уравнений Вольтер-ра I рода неявным методом Рулге-Кутты высокого порядка.

- Конф. молодых ученых ВЦ СО РАН, Новосибирск, 1996, стр.109-126.

2. Савченко А.О. Решение интегральных уравнений Воль-терра II рода неявным одношаговым методом Рунге-Кутты.

- Конф. молодых ученых ВЦ СО РАН, Новосибирск, 1995, стр.155-159.

3. Savchenko А.О. Numerical solution of Volterra integral equations of the second kind by implicit Runge-Kutta method.

- Bulletin of Novosibirsk Computer Center, Numerical Analysis. 7, 1996, p.87-92.

4. Savchenko A.O. Numerical solution of Volterra integral equations of the first kind by implicit Runge-Kutta method of high accuracy. - Тезисы доклада на международной конф."СМСР -96", Dubna, Russia, 1996. стр. 124.

5. Савченко А.О. Устойчивость решения интегральных уравнений Вольтерра неявным методом Рунге-Кутты. - Конференция молодых ученых ВЦ СО РАН, 1997, 9 стр. В печати.

6. Savchenko А.О. Numerical solution of Volterra integral equations of the first kind by implicit Runge-Kutta method of high accuracy - Bulletin of the Novosibirsk Computer Center. Numerical analysis, 20 p. To appear.