Практикум j практическое примщенше численных методов

3 = 1
где yi — приближенное решение в узле x_i} i = 1,2,..., п.
В проекционных методах приближенное решение интегрального уравнения (11.2) ищется в виде
п
у(^х) = (Н.8)
г=1
где (fi(x), i = 1,2,..., п — заданные линейно независимые функции, которые называются координатными. Часто удобнее ориентироваться па несколько отличное от (11.8) представление приближенного решения:
у(х) = f(x) + ^2dtpi{x), (11.9)
г=1
Метод проекционного типа характеризуется выбором координатных функций
г = 1,2, ...,п и способом определения вектора неизвестных коэффи­циентов с = {с_ьс₂,... ,с_п}. Отметим некоторые возможности по нахождению коэффициентов в представлении (11.8), (11.9).
При использовании представления (11.8) определим невязку
п
г(ж,с) =
рЬ ^п
я)-А / K(x,s)'£_lc_jtp_j(s)da- f( х).
i=i ^Ja j=i в методе наименьших квадратов постоянные с*, г = 1,2,..., п находятся из минимума квадрата нормы невязки в L₂(a, 6), т.е.
f^b
J(c) = / r²(x,c)dx —> min, с € Rⁿ. J a
Для определения c*, i = 1,2,...,n получим систему линейных алгебраиче­ских уравнений
(11.10)
3 = 1
где
рь рь
ац = / K(x,s)ipi(s)ds)x
J a J а
x(— А / K(x,s)(pj(s)ds)dx,
Ja
Jpb pb
/(^)((/?i(x) — Л / K(x,s)
1,2, ...,n.
a «/a
Тем самым матрица системы (11.10) симметрична.
В методе Галеркина коэффициенты с*, г = 1,2,..., п определяются из усло­вия ортогональности в L₂(a, 6) невязки г (а;, с) функциям
1,2,..., n:
/ г(а?,с)<дДа;)^я, г = 1,2, ...,п.
•/ a
В этом случае имеем систему линейных уравнений (НЛО), в которой
^aij— [ (ч>з(^х) - А [ K(x,s)(f_j(s)ds)ipi(x)dx,
«/a Ja
bi = / f{x)ipi(x)dx, г = 1,2,... ,n.
J a
Отметим среди проекционных методов и метод коллокации. В этом случае на отрезке [а, 6] выбирается п точек коллокации ж*, г = 1,2,... , п и коэффи­циенты , г = 1,2,... ,п в представлении (11.8) (или (11.9)) выбираются так, Нто невязка обращалась в нуль в точках коллокации, т.е.
r(xi, с) = 0, г = 1, 2,... , п.
Для коэффициентов матрицы и правой части системы (11.10) при использо- вании представления (11.8) получим
Oij = Vj(xi) -A / К(х_г, s)ip_j(s)ds,
J а
bi = f(xi), i = 1,2,,.. ,п.
Для решения системы линейных алгебраических уравнений (11.10) примени- ются прямые или итерационные методы.
Интегральные уравнения
с переменными пределами интегрирования
При приближенном решении интегрального уравнения Вольтерра второго рода (11.5) используется как метод квадратур, так и проекционные мето­ды. Для определенности, будем считать, что х_г = а, х_п = Ь. Для точек i = 1,2,..., п из (11.5) получим
u(xi)-X K(x_ils)u(s)ds = f(xi), г = 1,2, ...,n. (11.11)
J a
Принимая во внимание то, что интегрировать необходимо по отрезку пере­менной длины, запишем используемую квадратурную формулу в виде
/Xi i
6(x)dx « ^2 ^c^Q(^xj)i i = 2,3,..., n.
3 = 1
Применение к (11.11) дает систему линейных уравнений
Vi-^J2^cf^K(^xi>^s})y3 = г = \,2,...,п. (11.12)
3 = 1
Отличительная особенность системы уравнений (11.12) состоит в том, что матрица ее коэффициентов треугольная. Это позволяет найти приближен­ное решение интегрального уравнения у_и у₂,..., у_п последовательно друг за другом по рекуррентным формулам в предположении, что все диагональные элементы матрицы ненулевые. Наиболее простые расчетные формулы при решении интегрального уравнения Вольтерра второго рода мы получим при использовании квадратурной формулы трапеций.
При численном решении интегрального уравнения первого рода (11.6) мож­но ориентироваться на использование метода квадратур. Подобно (11.11) из (11.6) будем иметь
f
i — 1,2,..., n,

K(xi, s)u(s)ds — ^((Tj),
что дает систему линейных алгебраических уравнений
Y cfK(xi, Sj)yj = f(xi), i = 1,2,..., n.
3 = 1
Для того чтобы решение этой системы существовало необходимо потребовать выполнение условия К(х,х) Ф 0.
П
и(х) + Г
J а


1 dK(x,s)
'I
К(х,х) дх


\(s)ds =


1 df
К(х, х) dx


(х),

ри численном решении интегральных уравнений часто полезно провести предварительное преобразование исходной задачи. Типичным примером яв­ляется приведение интегрального уравнения Вольтерра первого рода к инте­гральному уравнению второго рода. Будем считать, что ядро и правая часть дифференцируемы К(х,х) ф 0. Тогда от уравнения (11.6) можно перейти к уравнению
которое представляет собой интегральное уравнения Вольтерра второго рода.
Интегральное уравнение Фредгольма первого рода
Интегральное уравнение (11.4) есть наиболее характерный пример некор­ректно поставленной задачи. Некорректность обусловлена тем, что при ма­лых возмущениях правой части f(x) не гарантируется малого возмущения решения.
Помимо (11.4) рассмотрим уравнение с возмущенной правой частью
f К(х, s)u(s)ds = /(.т), #б[а, Ь]. (11.13)
J а
Ядро K(x,s) есть вещественная непрерывная функция двух аргументов, а f(x),f(x) е L₂{a,b), причем
||/(х) -/(х)|| < 5,
при использовании обозначений
\\и(х)\\ = ((и,и)у/², (г*,г?)= f u(x)v(x)dx.
J а
При 8 —> 0 норма погрешности решения ||£t(or) — гл(ж)|| не стремиться к нулю. Определим линейный интегральный оператор
Ay = j K(x,s)y(s)ds, xe[a,b\. (11.14)
J a
Задачу с неточной правой частью (11.13) запишем в виде операторного урав­нения первого рода
В методе регуляризации Тихонова приближенное решение задачи (11.15) на­ходится из минимума сглаживающего функционала:
J_a(y) -> min, у е Ь₂(а, Ь), (11.16)
где
Ш = \\Ау-1\\² + Ф\\²,
аа>0 - параметр регуляризации.
Обозначим решение задачи (11.16) через у_а. Оно може^гг быть найдено как решение уравнения Эйлера для вариационной задачи (11.16)
<*у_а + А* Ау_а = A*f,
где
А*у = / K(s,x)y(s)ds, х Е [а,Ь].
J а
Тем самым приходим к интегральному уравнению Фредгольма
ау_а + ! G(x,s)y_a(s)ds = ф(х), же [а,6] J а
G
с симметричным ядром

(x,s)= f K(t,x)K(t,s)dt J а
ф
и правой частью

(х)= f K(s,x)f(s)ds.
J а
Принципиальный момент в методе регуляризации связан с выбором пара­метра регуляризации а, его согласованием с погрешностью входных данных. При использовании принципа невязки параметр регуляризации выбирается из условия
I\Ау_а - /|| = <5.
При таком выборе а = а(5) норма погрешности ||у_а — ?х|| —> 0 при 5 —> 0, т.с. приближенное решение стремится к точному решению задачи.

3. 
   Упражнения
   

Упражнение 11.1 Напишите программу для числеппого решения инте­грального уравнения Фредгольма второго рода методом квадратур с исполь­зованием квадратурной формулы трапеций при равномерном разбиении ин-
т
ервала интегрирования. С помощью этой программы найдите приблиэюен- пое решение интегрального уравнения
при различном числе частичных отрезков.
В модуле fredholm функция fredholmO обеспечивает приближенное решение интегрального уравнения
и (х) - K(x,s)u(s) = f(x), xe[a,b]
при использовании квадратурной формулы трапеций в методе квадратур. Для решения системы линейных уравнений используется функция solveLUO из модуля lu.
import пшпру as пр
from lu import solveLU
def fredholm(k, f, a, b, n):
и и и
Solution Fredholm integral equation of the second kind. k(x,s) is the kernel of the integral equation, f(x) is the right part, 0 < x,s < b.
Method with trapezoidal quadrature formula.
и и и
h = (b - a) / n A = np.identity(n+1, ’float’) r = np.zeros((n+1), ’float’) for i in range(n+l): x = a + i*h
A[i,0] = A[i,0] - k(x,a)*h/2 for j in range(l,n): _s = _a + j*h
A[i,j] = A[i,j] - k(x,s)*h A[i,n] = A[i,n] - k(x,b)*h/2 r [i] = f (x) у = solveLU(A, r) return у
Для приближенного решения модельного интегрального уравнения при ис­пользовании различного числа узлов г = 1,2,п используется следующая программа.
д If : _;;_:м ^^:^ v: ■:^::; ;; ■ ■^::::_:: ;_:-
i



Рис. 11.1 Приближенное решение при числе узлов гг = 6,11,51

mport numpy as np import matplotlib.pyplot as pit from fredholm import fredholm def k(x,s):
return l./(np.pi*(l. + (x-s)**2)) def f (x):
return 1. a = -1. b = 1.
nList = [5, 10, 50] sglist = [’-’, ’]
for kk in range(len(nList)): n = nList[kk]
x = np.linspace(a, b, n+1) у = fredholm(k, f, a, b, n) si = ^,n=^> + str(n+l) sg = sglist[kk] pit.plot(x, y, sg, label=sl) pit ^label^x’)
pit.grid(True) pit.legend(loc=0) pit.show()
Сходимость приближенного решения при увеличении числа разбиений п де­монстрируется рис. 11.1.
У
пражнение 11.2 Напишите программу для численного решения инте­грального уравнения Фредгольма первого рода методом квадратур с исполь­зованием квадратурной формулы прямоугольников при равномерном разбие­нии интервала интегрирования с симметризацией матрицы системы урав­нений и ее итерационном решении методом сопряэ1сениых градиентов. С по­мощью этой программы найдите приблиэюенное решение интегрального уравнения
при различной точности приблиэюенного решения системы линейных урав­нений и сравните его с точным решением интегрального уравнения.
После дискретизации приходим к системе уравнений
Ay = f.
Итерационный метод сопряженных градиентов применяется для приближен­ного решения системы уравнений
A* Ay = A* f.
В
>b


модуле fredholml функция fredholml О обеспечивает приближенное реше­ние интегрального уравнения
при использовании квадратурной формулы прямоугольников в методе квад­ратур. Для решения системы линейных уравнений используется функция eg () из модуля eg.
import numpy as np from eg import eg
def fredholml(k, f, a, b, n, tol = 1.0e-9):
II II II
Solution Fredholm integral equation of the first kind. k(x,s) is the kernel of the integral equation, f(x) is the right part, 0 < x, s < b.
CG iterative method with rectangle quadrature formula.