Материал: Спецглавы высшей математики. численные методы. Пантелеев И.Н

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

Рис. X.2. Численное решение интегрального

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

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

Это замечание иллюстрирует рис. X.3, на котором представлено решение уравнения из примера X.1, полученное при использовании квадратурной формулы Симпсона. Далее приведен листинг файла Solve2_g11.m (листинг X.4), содержащего описание соответствующей функции.

Листинг X.4. Файл Solve2_g11.m function [X,Y]=solve2_g11(x1,x2,N,Lambda) % задание временной сетки h=(x2-x1)/(N-1);

i=1:N; t(i)=x1+h*(i-1); s=t;

% задание коэффициентов квадратурной формулы метода Симпсона

A(1)=1/3;

A(N)=1/3;

k=0;

for i=2:N-1 if k==0

A(i)=4/3;

k=1; else

A(i)=2/3;

k=0;

end;

end;

%вычисление значений функции Q(t,s) в узлах сетки for i=1:N

for j=1:N q(i,j)=Q11(t(i),s(j));

end;

end;

%вычисление значений функции f(t) в узлах временной сетки

F=f11(t);

for i=1:N for j=1:N

if i==j

M(i,j)=1-Lambda*A(i)*q(i,j)*h; else

M(i,j)=-Lambda*A(i)*q(i,j)*h; end;

end;

end;

X=t;

Y=M^-1*F'; % нахождение решения интегрального уравнения

Для получения численного решения интегрального уравнения следует создать файлы q11.m, f11.m, Solve1_g11 и затем выполнить следующую последовательность команд:

>>x1=1; % левая граница отрезка поиска решения

>>x2=2; % правая граница отрезка поиска решения

>>Lambda=1;

>>N=300; % число узлов разбиения отрезка

>>[X,Y]=solve2_g11(x1,x2,N,Lambda);

>>plot(X,Y) % визуализация решения интегрального уравнения

%(рис. X.3)

Рис. X.3. Численное решение интегрального уравнения

Симпсона

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

Так как линейные интегральные уравнения Вольтерра в отличие от уравнения Фредгольма имеют единственное непрерывное решение при любых значениях параметра λ, при нахождении численного решения уравнения

a

где t [а, b], можно положить λ = 1.

Учитывая, что уравнение Вольтерра формально можно считать уравнением Фредгольма вида

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

Введем в рассмотрение временную сетку s j [a, b] , состоящую

из п узлов, и выберем конкретную квадратурную формулу с весами Аj, тогда приближенное решение интегрального уравнения принимает вид (X.17). Составим систему линейных алгебраических уравнений, аналогичную системе (X.18), которая в силу свойств ядра интегрального уравнения (X.21) вырождается в треугольную:

Из (X.22) видно, что искомые значения x1, x2,…,хп находятся последовательными вычислениями по следующим формулам:

x1 f1 , (X.23) 1 A Q

1 11

i 1

fi Aj Qij x j

Пример X.2. Решение в пакете МАТLАВ интегрального уравнения

1.Создайте файл q11_2.m (листингX.5), содержащий описание функции, возвращающей значения подынтегральной функции. Листинг X.5. Файл q11_2.m

function z=Q11_2(t,s) z=t*cos(t*s.^3).^2;

2.Создайте файл F11_2.m (листинг X.6), содержащий описание функции, возвращающей значения функции f(t).

Листинг X.6. Файл F11_2.m function z=f11_2(t) z=t.^2-1/3*tan(t.^4);

3.Создайте файл Solve3_g11.m (листинг X.7), содержащий описание функции, возвращающей решение интегрального уравнения.

Листинг X.7. Файл Solve3_g11.m function [T,Y]=Solve3_g11(t1,t2,N) % задание временной сетки h=(t2-t1)/(N-1);

i=1:N; t(i)=t1+h*(i-1);

241

s=t;

% задание коэффициентов квадратурной формулы метода трапеций

A(1)=0.5; m=2:N-1; A(m)=1; A(N)=0.5;

%вычисление значений функции Q(t,s) в узлах сетки for i=1:N

for j=1:N Q(i,j)=Q11_2(t(i),s(j));

end;

end;

%вычисление значений функции f(t) в узлах временной сетки

F=f11_2(t);

%вычисление решения интегрального уравнения

x(1)=F(1)./(1-A(1)*Q(1,1)); for m=2:N

S=0;

for j=1:m-1 S=S+h.*A(j).*Q(m,j).*x(j);

end; x(m)=(F(m)+S)./(1-h.*A(m).*Q(m,m));

end;

T=t;

Y=x;

4. Выполните следующую последовательность команд:

>>t1=0; % левая граница отрезка поиска решения

>>t2=0.5; % правая граница отрезка поиска решения

>>N=300; % число узлов разбиения отрезка

>>[X,Y]=solve3_g11(x1,x2,N,Lambda);

>>plot(X,Y) % визуализация решения интегрального уравнения

%(рис. X.4)

242

Рис. X.4. Численное решение интегрального уравнения

%оценка параметров решения интегрального уравнения

%методом наименьших квадратов

>>Z=[ones(size(X'))X'X'.^2];

>>a=Z\Y'

a=

-0.0000 -0.0000 1.0003

>>format long

>>a

a=

1.0е+002*

-0.00000187595514 -0.00000470664690 1.00030838604155

243

Получим расчетные формулы для решения уравнения Вольтерра первого рода (X.5) при использовании метода средних прямоугольников. Решения уравнения будем находить в узлах временной сетки

a

Рис. X.5. Пространственно-временная сетка, используемая для решения уравнения Вольтерра

Из (X.26) видно, что в данном случае условие совпадения узлов квадратур si с узлами временной сетки ti не является обязательным, поэтому их можно выбрать посередине элементарных промежутков интегрирования [ti-1, ti] (рис. X.5). Выбор данной сетки означает, что

Учитывая выбор квадратурной формулы и условия (X.27), запишем (X.26) в следующем виде:

hQ(t1 , s1 )x1 f (t1 ),

hQ(t2 , s1 )x1 hQ(t2 , s2 )x2 f (t2 ),

244

Листинг X.10. Файл Volterra11.m function [T,Z]=Volterra11(t1,t2,N)

%задание временной сетки h=(t2-t1)/(N-1);

i=1:N;

t(i)=t1+h*i;

s=t+h/2;

%вычисление значений ядра интегрального уравнения

%в узлах временной сетки

fоr i=1:N fоr j=1:N

q(i+1,j+1)=Q11_4(t(i),s(j)); end;

end;

%вычисление значений функции F(t) в узлах временной сетки

F(1)=F11_4(t(1)); F(i+1)=F11_4(t(1));

%Вычисление значений решения интегрального уравнения в соответствие с (X.29),(X.30)

х(1)=F(1)/(h*q(1,1)); for m=2:N

s=0;

fоr j=2:m-1 s=s+q(m,j)*m(j);

еnd;

х(m) = (F(m) -h*s)/( h*q(m,m)); end;

Т=t; Z=х;

4. Выполните следующую последовательность команд:

>>t1=1; % левая граница интервала поиска решения

>>t2=1.3;%правая граница интервала поиска решения

>>N=300; % число интервалов разбиения отрезка

%[t1,t2]

>>[Т,Z]=Voltегга11(t1,t2,N); % численного решения

% интегрального уравнения Вольтерра