Линейное уравнение Фредгольма второго рода

Линейное интегральное неоднородное уравнение Фредгольма второго рода имеет вид:

где ядро определено в квадрате V = [a,b]*[a,b]. Кроме того, полагается, что ядро непрерывно в V. При = 1, используя квадратурную формулу трапеций с постоянным шагом h, получим:

где

• n = (b-a)/h+1 - целое
• A_j  = 1 при j, не равном 1 или n
• A_j  = 0.5 при j, равном 1 или n.

Для сравнения с нолем в алгоритм передается малое число epsilon, и любое число, по модулю меньшее epsilon, считается нолем.

Если нашли ошибку в алгоритме - сообщите!

Реализации алгоритма на различных языках:

Блоксхемы:

Реализация алгоритма на AlgoPascal:

unit Fredholm2Unit;

declaration
    function F(X : Real):Real;
    function K(X : Real; S : Real):Real;

interface SolveFredholm2;

implementation

(******************************************************************
Процедура  решает  интегральное  уравнение Фредгольма второго рода,
заданное   ядром   интегрирования   K(X,S) и правой частью F(X), на
отрезке [A, B].

Результат помещается в массив Y с номерами элементов от 1 до N, где
1 соответствует A, N соответсвует B.

Epsilon - малое число, передаваемое для сравнения с нолем   в  ходе
решения получаемой системы уравнений.
******************************************************************)
function SolveFredholm2(
    const A : Real;
    const B : Real;
    const N : Integer;
    out   Y : array of Real;
    const   Epsilon : Real):Boolean;
var
    h   :   Real;
    t   :   Real;
    m1  :   Real;
    x   :   Real;
    
    SMat:   array of array of Real;
    
    i   :   Integer;
    j   :   Integer;
    u   :   Integer;
    k1  :   Integer;
    m   :   Integer;
begin
    SetBounds(SMat, [1, N], [1, N+1]);
    SetBounds(Y,    [1, N]);
    
    h:=(b-a)/(n-1);
    i:=1;
    repeat
        x:=a+(i-1)*h;
        SMat[i,n+1]:=F(x);
        j:=1;
        repeat
            SMat[i,j]:=-h*K(x,a+(j-1)*h);
            if (j=1)or(j=n) then
            begin
                SMat[i,j]:=SMat[i,j]/2;
            end;
            if j=i then
            begin
                SMat[i,j]:=1+SMat[i,j];
            end;
            j:=j+1;
        until not(j<=n);
        i:=i+1
    until not(i<=n);


    SetBounds(Y, [1,N]);
    Result:=True;
    for i:=1 to n do
    begin
        k1:=i;
        m1:=AbsReal(SMat[i,i]);
        for j:=i+1 to n do
        begin
            if m1<AbsReal(SMat[j,i]) then
            begin
                m1:=AbsReal(SMat[j,i]);
                k1:=j;
            end;
        end;
        if AbsReal(m1)>=Epsilon then
        begin
            for j:=i to n+1 do
            begin
                t:=SMat[i,j];
                SMat[i,j]:=SMat[k1,j];
                SMat[k1,j]:=t;
            end;
            for k1:=i+1 to n do
            begin
                t:=SMat[k1,i]/SMat[i,i];
                SMat[k1,i]:=0;
                for j:=i+1 to n+1 do
                begin
                    SMat[k1,j]:=SMat[k1,j]-t*SMat[i,j];
                end;
            end;
        end
        else
        begin
            Result:=False;
            Break;
        end;
    end;
    if Result then
    begin
        i:=n;
        repeat
            y[i]:=SMat[i,n+1];
            j:=i+1;
            while j<=n do
            begin
                y[i]:=y[i]-SMat[i,j]*y[j];
                j:=j+1;
            end;
            y[i]:=y[i]/SMat[i,i];
            i:=i-1
        until not(i>=1);
    end;
end;


end.

Бочканов Сергей, Быстрицкий Владимир
Copyright © 1999-2004
При поддержке проекта MANUAL.RU