Методическое пособие 521

Решение.

На промежутке [0, 2] рассмотрим три узловые точки: на концах x1 =0, x3 =2 и в его середине x2 =1. Для вычисления интеграла воспользуемся простой формулой Симпсона:

2

(x t)y(t)dt 13 xy1 4(x 1)y2 (x 2)y3 ,

0

и исходное уравнение примет вид

квадратурной формулы x1, x2 и x3, придем к системе трёх уравнений относительно узловых переменных y1, y2 и y3:

y1 (4y2 2y3)3, y2 1 (y1 y3)3, y3 8 (2y1 4y2)3.

или

3y1 4y2 2y3 0, y1 3y2 y3 3,

2y1 4y2 3y3 24.

Найдем решение этой линейной системы: y1 =4, y2 =1, y3 =4.

Используя соотношение (8), получим приближенное решение задачи

y(x) x3 4x 4.

Поскольку данное интегральное уравнение может быть решено аналитически точно, то есть возможность дать оценку точности полученного решения. Точное решение легко полу-

160

чается путем двукратного дифференцирования исходного уравнения и последующего интегрирования линейного дифференциального уравнения; оно равно

y(x) = x3 –(132/35)x+128/35=x3 –3,7714x+3,6571.

В узловых точках имеем: y(0)=3,657, y(1)=0,886, y(2)=4,114. Максимальная относительная погрешность превышает 12 %. Чтобы повысить точность, промежуток [0,2] разобьем на два равных отрезка и на каждом при интегрировании будем использовать простую формулу Симпсона, в результате чего получим так называемую составную формулу Симпсона:

Обозначая y1 =y(0), y2 =y(0,5), y3 =y(1), y4 =y(1,5), y5 =y(2)

и применяя данную квадратурную формулу, получим

y(x) x3 16 xy1 4(x 0,5)y2 2(x 1)y3 4(x 0,5)y4 (x 2)y5

и тем самым систему линейных уравнений для нахождения узловых переменных

y1 (13)y2 (13)y3 y4 (13)y5 0,

(112)y1 y2 (16)y3 (23)y4 (14)y5 18, (16)y1 (13)y2 y3 (13)y4 (16)y5 1,

(14)y1 (23)y2 (16)y3 y4 (112)y5 278, (13)y1 y2 (13)y3 (13)y4 y5 8.

Решение этой системы y1 =3,679, y2 =1,911, y3 =0,893, y4 =1,375, y5 =4,107 при подстановке в (7.8) дает функцию

161

y(x) = x3 –3,786x+3,679,

являющуюся значительно более точным приближением задачи, чем полученное ранее с тремя узлами. Максимальная погрешность уже не превосходит 1%.

Если взять в 2 раза больше частичных отрезков, т.е. n=4, и провести аналогичную процедуру, получим еще более точное решение задачи

y(x) = x3 –3,7723x+3,6585.

Очевидно, метод сходится. Составим таблицу, в которой представлены значения численного решения при различных n, а также точного решения в ряде точек отрезка [0,2].

Порядок решения в системе Maple

Составим программу решения интегрального уравнения общего вида

b

y(x) f (x) K(x,s)y(s)ds

a

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

> restart;

Введем переменную _meth, указывающую метод численного интегрирования: её значение 1 соответствует формуле трапеций, 2 – формуле Симпсона. Задает число делений отрезка в составной квадратурной формуле.

> _meth:=2; n:=4;

162

Смысл следующих переменных ясен из обозначений: a и b – границы основного промежутка – пределы интегрирования, f – свободный член интегрального уравнения, K – ядро уравнения

>a:=0.;b:=2.; f:=x->x^3;:=(x,t)->x-t; lambda:=1;

Шаг – расстояние между соседними узлами

>h:=(b-a)/n; m:=n+1;

>if(_meth=2) then m:=2*n+1; h:=(b-a)/(2*n); fi;

>z:=vector(n+1); X:=[(a+(k-1)*h)$k=1..n];

Далее вычисление конечной суммы – по методу трапеций

> if(_meth=1) then finsum:=x->add((K(x,X[j])*z[j] +K(x,X[j+1])*z[j+1])*h/2,j=1..m-1);fi;

и по методу Симпсона

if(_meth=2) then finsum:=x-> add((K(x,X[j])*z[j]+ 4*K(x,X[j+1])*z[j+1]+K(x,X[j+2])*z[j+2])*h/3, j=[seq(l,l=1..m-2,2)]);fi;

Составляем дискретные уравнения

> for i from 1 to m do g[i]:=z[i]-f(X[i])- lambda*finsum(X[i]);od;

Решаем полученную систему относительно узловых переменных

r:=solve({g[k]$k=1..m},{z[k]$k=1..m});

Y:=subs(r,[z[k]$k=1..m]);

Вывод окончательного решения интегрального уравнения phi:=f(x)+lambda*subs(r,finsum(x));

Графическое представление решения

>gr2:=plot([[X[j],Y[j]]$j=1..m],x=X[1]..X[m], style=POINT,symbol=CIRCLE,color=black):

>gr6:=plot(phi,x=a..b):

>plots[display](gr2,gr6);

Чтобы провести анализ точности, зададим как готовую функцию аналитического решения

> u:=unapply(x^3-(132/35.)*x+128/35.,x);

Протабулируем функцию phi вместе с точным решением; – относительная погрешность.

> for s from 1 to n do printf(`x=%5.3g z=%8.6g u=%8.6g delta=%g\n`,X[s],Y[s],u(X[s]),

163

> a:=0.; b:=3.; m:=8; h:=(b-a)/m; f:=x->1; K:=(x,t)->(x-t); lambda:=1;

if(_meth=1)then n:=m+1;else n:=2*m+1;h:=h/2;fi;

Чтобы свести задачу к уравнению Фредгольма, ядро берем равным нулю при x < t:

> Kt:=(x,t)->(`if`(x-t>=0,K(x,t),0.));

z– массив узловых значений искомой функции.

>z:=vector(n+1);

>X:=[(a+(k-1)*h)$k=1..n]; # узловые точки xk

Формируем систему линейных уравнений и решаем её:

> for i from 1 to n do

if(_meth=1) then g[i]:=z[i]-f(X[i])-lambda* add((Kt(X[i],X[j])*z[j]+

Kt(X[i],X[j+1])*z[j+1])*h/2,j=1..i-1); else g[i]:=z[i]-f(X[i])-lambda*

( add((Kt(X[i],X[j])*z[j]+ 4*Kt(X[i],X[j+1])*z[j+1]+

Kt(X[i],X[j+2])*z[j+2])*h/3, j=[seq(l,l=1..i-2,2)])+

`if`(i mod 2=0,(Kt(X[i],X[i-1])*z[i-1]+ Kt(X[i],X[i])*z[i])*h/2.,0) );

fi; od;

>r:=solve({g[k]$k=1..n},{z[k]$k=1..n});

>Y:=subs(r,[z[k]$k=1..n]);

Получили дискретные значения функции y(x) в равноотстоящих точках интервала [0, 3]. Сравним этот результат с точным решением, которое получим операционным методом.

>r1:=inttrans[laplace](y(x)=f(x)+lambda* int((x-t)*y(t),t=0..x),x,p);

>inttrans[invlaplace](solve(r1,

laplace(y(x),x,p)),p,x);

165

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

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

7.4.2.Метод взвешенных невязок

Вметоде взвешенных невязок приближенное решение ~y(x) интегрального уравнения ищется в виде суммы f(x) и

линейной комбинации линейно независимых на отрезке [a, b] базисных функций N1(x), N2(x), …, Nm(x), т.е.

i 1

где 1, 2, …, m – параметры аппроксимации – некоторые числа.

Взяв для определенности уравнение (7.3) и подставив в него (7.9), получим невязку

где

x

i (x) K(x,t)Ni (t)dt , i = 1, …, m.

a

Задача состоит в том, чтобы подобрать значения параметров 1, 2, …, m, чтобы невязка была как можно ближе к нулю.

Для этого согласно методу взвешенных невязок ставится условие ортогональности невязки к (весовым) функциям W1(x),

167

симых функций. Эти условия дают следующую систему линейных уравнений

b

относительно параметров 1, 2, …, m. Найденные в результате решения этой системы числовые значения параметров после подстановки в (9) дают приближенное решение исходного интегрального уравнения.

В качестве весовых функций Wi(x), i = 1, …, m должны выступать линейно независимые функции, а, следовательно, за такие функции можно взять и базисные функции Ni(x). Именно такой вариант метода взвешенных невязок чаще всего и используется на практике. При этом в литературе он может иметь название метода Галёркина либо метода моментов.

Замечания. 1) Метод взвешенных невязок сходится, если с увеличением с m функция ym(x) становится ближе к точному решению интегрального уравнения y(x), т.е.

ym (x) y(x) 0 . Но поскольку точное решение, как прави-

m

ло, неизвестно, на практике для проверки сходимости прове-

ряют условие ym 1(x) ym (x) 0.

m

2) Метод будет сходиться и в случае, когда в аппроксимации (9) нет функции f(x).

3) Если решается уравнение Вольтерра (с переменным верхним пределом x), то верхний предел интегрирования b в условии (10) в известной степени определяется произвольно как правая граница промежутка a x b, на котором ищется решение уравнения.

Пример. Решить уравнение

x

y(x) 1 x (x t)y(t)dt .

0

168

Решение.

Полагая m=3, в соответствии с (7.9) используем аппроксимацию

~y y3(x) 1 x 1 2x 3x2 .

Тогда невязка имеет вид

x

R 1 2x 3x2 (x t)(1 x 1 2x 3x2)dt

0

1 2x 3x2 12 1x2 16 2x3 121 3x4 12 x2 16 x3.

Будем искать решение уравнения на отрезке [0, 1]. Из условия ортогональности невязки R к функциям 1, x, x2 на этом отрезке имеем

Подставляя сюда выражение для невязки R, а после вычисляя интегралы и совершая некоторые преобразования, получим систему линейных алгебраических уравнений для определения 1, 2 и 3:

Тем самым приближенное решение исходного интегрального уравнения принимает вид

y3(x) 1,0131 0,8510x 0,8392x2 .

Если в аппроксимации (9) взять четыре (m=4) параметра, то получим

169