使用高斯方法在 Turbo Prolog 2.0 上查找三次样条的系数

我会回答我自己的问题。

3点例子的决策逻辑：

1. 让我们确定段数。
   相关性很明显：段数 = 点数 - 1。
   三个点的段数 - 2（从 x1 到 x2 以及从 x2 到 x3）。
   
2. 让我们确定样条线的数量。
   既然有 2 段，那么：
   (1) a1 x^3 + b1 x^2 + c1 x + d1 for x in [x1,x2);
   (2) a2 x^3 + b2 x^2 + c2 x + d2 对于 [x2,x3] 中的 x；
   我们需要找到系数 a1, b1, c1, d1, a2, b2, c2, d2。
   
3. 我们形成具有以下行的方程组：
   （1）三次多项式的行（按横坐标），等于相应纵坐标的值。
   a1 x1^3 + b1 x1^2 + c1 x1 + d1 = y1
   a1 x2^3 + b1 x2^2 + c1 x2 + d1 = y2
   a2 x2^3 + b2 x2^2 + c2 x2 + d2 = y2
   a2 x3 ^3 + b2 x3^2 + c2 x3 + d2 = y3
   (2) 两个多项式在公共点的一阶导数相等的行。
   3 a1 x2^2 + 2 b1 x2 + c1 = 3 a2 x2^2 + 2 b2 x2 + c2
   (3) 两个多项式在公共点的二阶导数相等的行。
   6 a1 x2 + 2 b1 = 6 a2 x2 + 2 b2
   (4) 二阶导数在边界点处设置为零的行。（在我们的例子中，边界点 (x1, y1), (x3, y3)）。
   6 a1 x1 + 2 b1 = 0
   6 a2 x3 + 2 b2 = 0
   
4. 我们根据上述方程形成一个 AX=B 形式的矩阵方程：
   
   
5. 我们用高斯方法求解这个方程。

Turbo Prolog 代码：

NOWARNINGS

% Раздел описания доменов
DOMAINS
    file = datafile % Файл datafile
    point = point(real, real) % Структура, описывающая точку
    points = point* % Список точек
    list_of_real = real* % Список вещественных чисел
    matrix = list_of_real* % Список списков вещественных чисел

% Раздел описания предикатов
PREDICATES
    length(points, integer)
    input(points)
    input_action(points, integer)
    read_points_from_console(points)
    read_points_from_file(points)
    quicksort(points, points)
    partition(points, real, points, points)
    merge(points, points, points)
    merge(list_of_real, list_of_real, list_of_real)
    merge(matrix, matrix, matrix)
    distinct(points, points)
    form_system(points, matrix)
    create_list(real, integer, list_of_real)
    create_equation(list_of_real, real, integer, integer, list_of_real)
    create_equation_1(point, integer, integer, list_of_real)
    create_equation_2(point, integer, integer, list_of_real)
    create_equation_3(point, integer, integer, list_of_real)
    create_equation_4(point, integer, integer, list_of_real)
    get_last(points, point)
    stage_1(integer, points, matrix)
    stage_1_1(integer, point, list_of_real)
    stage_1_2(integer, points, matrix)
    stage_1_2(integer, points, matrix, integer)
    stage_1_3(integer, point, list_of_real)
    stage_2(integer, points, matrix)
    stage_2(integer, points, matrix, integer)
    stage_3(integer, points, matrix)
    stage_3(integer, points, matrix, integer)
    stage_4(integer, points, matrix)
    gaussian_elimination_with_backsubstitution(matrix, list_of_real)
    gaussian_elimination(matrix, matrix, matrix)
    backsubstitution(matrix, list_of_real, list_of_real)
    backsubstitution_auxillary(list_of_real, list_of_real, real, real)
    pivot_row(matrix, list_of_real, matrix)
    normalized_pivot_row(list_of_real, list_of_real)
    normalized_pivot_row_auxillary(list_of_real, real, list_of_real)
    normalized_other_rows(matrix, list_of_real, matrix)
    normalized_other_row(list_of_real, real, list_of_real, list_of_real)
    output(points, list_of_real)
    output_action(points, list_of_real, integer)
    write_coefficients(points, list_of_real)

% Раздел описания внутренней цели
GOAL
    input(Points), % Ввод точек
    quicksort(Points, SortedPoints), % Быстрая сортировка списка точек [O(N log N)]
    distinct(SortedPoints, SortedPointsWithoutDuplicates), % Удаление точек с повторяющимися X [O(N)]
    form_system(SortedPointsWithoutDuplicates, EquationSystem), % Формирование системы уравнений [O(N^2)]
    gaussian_elimination_with_backsubstitution(EquationSystem, Result), % Решение системы уравнений [O(N^3)]
    output(SortedPointsWithoutDuplicates, Result), % Вывод результата
    readchar(_). % Завершение программы после нажатия клавиши клавиатуры

% Раздел описания предложений
CLAUSES
    % Вспомогательные общие предикаты
    % Длина списка
    length([], 0).
    length([_|Xs], N) :-
        length(Xs, P),
        N = P + 1.

    % Ввод точек
    input(Points) :-
        write("-== MENU ==-"), nl,
        write("1. Read from console;"), nl,
        write("2. Read from file."), nl,
        write("Another button to exit"), nl,
        readint(C), % Считывание целого числа
        input_action(Points, C), % Ввод списка точек
        length(Points, Length), Length >= 3;
        write("There must be at least three points!"), fail. % Если меньше 4 точек

    % Ввод из диалогового окна
    input_action(Points, 1) :-
        read_points_from_console(Points).

    % Ввод из файла
    input_action(Points, 2) :-
        write("File name: "),
        readln(FileName), % Ввод названия файла
        existfile(FileName), % Существует ли файл
        openread(datafile, FileName), % Открытие файла для чтения
        readdevice(datafile), % Перенаправление ввода на файл
        read_points_from_file(Points), % Чтение точек из файла
        closefile(datafile), !; % Закрытие файла
        write("Error reading file!"), nl, fail.

    % Чтение точек с консоли
    read_points_from_console([point(X, Y)|Tail]) :-
        write("X: "), readreal(X), % Чтение X
        write("Y: "), readreal(Y), % Чтение Y
        read_points_from_console(Tail).
    read_points_from_console([]).

    % Чтение точек с файла
    read_points_from_file([Point|Tail]) :-
        readterm(point, Point), % Чтение терма point(X, Y)
        read_points_from_file(Tail).
    read_points_from_file([]).
        % Пример содержимого файла:
        % point(2, 2)
        % point(1, 1)
        % point(3, 3)

    % Быстрая сортировка
    quicksort([point(X, Y)|Xs], Zs) :-
        partition(Xs, X, Left, Right), % Разбиение
        quicksort(Left, Ls), % Обработка левой части
        quicksort(Right, Rs), % Обработка правой части
        merge(Ls, [point(X, Y)|Rs], Zs). % Объединение списков
    quicksort([], []).
        % Пример (для integer):
        %    quicksort([5, 3, 1, 2, 4], X)
        % -> X = [1, 2, 3, 4, 5]

    % Разбиение списка
    partition([point(X, Y)|Xs], Z, Ls, [point(X, Y)|Rs]) :-
        X > Z,
        partition(Xs, Z, Ls, Rs).
    partition([point(X, Y)|Xs], Z, [point(X, Y)|Ls], Rs) :-
        X <= Z,
        partition(Xs, Z, Ls, Rs).
    partition([], _, [], []).
        % Примеры (для integer):
        %    partition([1, 2, 3, 4, 5], 2, L, X)
        % -> L = [1, 2], X = [3, 4, 5]
        %    partition([1, 2, 3, 4, 5], 4, L, X)
        % -> L = [1, 2, 3, 4], X = [5]

    % Объединение двух списков
    merge([H|Xs], Zs, [H|Ts]) :- merge(Xs, Zs, Ts).
    merge([], Zs, Zs).
        % Пример (для integer):
        %    merge([1, 2, 3], [4, 5, 6], X)
        % -> X = [1, 2, 3, 4, 5, 6]

    % Удаление дубликатов точек по X [O(N) - нужна предварительная сортировка]
    distinct([], []). % В пустом списке нет дубликатов
    distinct([Point], [Point]). % В списке из одной точки нет дубликатов
    distinct([point(X, Y1), point(X, Y2)|Tail], Result) :- % Если точки одинаковы по X
        distinct([point(X, Y1)|Tail], Result). 
    distinct([Point1, Point2|Tail], [Point1|Result]) :- % Если точки разные по X
        distinct([Point2|Tail], Result).
        % Пример (для integer):
        %    distinct([1, 1, 2, 2, 2, 3, 3, 4], X)
        % -> X = [1, 2, 3, 4]
    % Комментарий: без предварительной сортировки потребуется алгоритм со сложностью O(N^2)

    % Формирование системы уравнений
    form_system(Points, Equations) :-
        length(Points, Length),
        % Строки кубических полиномов (по абсциссам), приравненных значениям соответствующих ординат
        stage_1(Length, Points, Equations1),
        % Строки с приравненными значениями первых производных двух полиномов в общей точке
        stage_2(Length, Points, Equations2),
        % Строки с приравненными значениями вторых производных двух полиномов в общей точке
        stage_3(Length, Points, Equations3),
        % Строки с приравненными нулю вторыми производными в граничных точках
        stage_4(Length, Points, Equations4),
        % Объединение четырех систем уравнений
        merge(Equations1, Equations2, Equations12),
        merge(Equations3, Equations4, Equations34),
        merge(Equations12, Equations34, Equations).

    % Создание списка длины N, заполненного значением X
    create_list(X, N, [X|L]) :-
        N > 0, % Выполняем, пока N > 0, записывая в голову значение Х
        Next = N - 1, % Получаем следующее значение N
        create_list(X, Next, L).
    create_list(_, N, []) :- N <= 0. % Рекурсия останавливается, когда/если N меньше 1

    % Создание уравнения
    create_equation(List, Y, ZerosBeforeCount, ZerosAfterCount, Equation) :-
        create_list(0, ZerosBeforeCount, ZerosBefore), % Создание левой части уравнения
        create_list(0, ZerosAfterCount, ZerosAfter), % Создание правой части уравнения
        merge(ZerosBefore, List, EqPart1), % Слияние левой части с заполненной
        merge(EqPart1, ZerosAfter, EqPart2), % Последующее слияние с правой частью
        merge(EqPart2, [Y], Equation). % Добавление Y в конец

    % Создание уравнения первого этапа
    create_equation_1(point(X, Y), ZerosBeforeCount, ZerosAfterCount, Equation) :-
        AX = X * X * X, BX = X * X,
        List = [AX, BX, X, 1], % List = [X^3, X^2, X, 1]
        create_equation(List, Y, ZerosBeforeCount, ZerosAfterCount, Equation).

    % Создание уравнения второго этапа
    create_equation_2(point(X, _), ZerosBeforeCount, ZerosAfterCount, Equation) :-
        AX = 3 * X * X, BX = 2 * X, NAX = -AX, NBX = -BX,
        List = [AX, BX, 1, 0, NAX, NBX, -1, 0], % List = [3*X^2, 2*X, 1, 0, -3*X^2, -2*X, -1, 0]
        create_equation(List, 0, ZerosBeforeCount, ZerosAfterCount, Equation).

    % Создание уравнения третьего этапа
    create_equation_3(point(X, _), ZerosBeforeCount, ZerosAfterCount, Equation) :-
        AX = 6 * X, BX = 2, NAX = -AX, NBX = -BX,
        List = [AX, BX, 0, 0, NAX, NBX, 0, 0], % List = [6*X, 2, 0, 0, -6*X, -2, 0, 0]
        create_equation(List, 0, ZerosBeforeCount, ZerosAfterCount, Equation).

    % Создание уравнения четвертого этапа
    create_equation_4(point(X, _), ZerosBeforeCount, ZerosAfterCount, Equation) :-
        AX = 6 * X, BX = 2,
        List = [AX, BX, 0, 0], % List = [6*X, 2, 0, 0]
        create_equation(List, 0, ZerosBeforeCount, ZerosAfterCount, Equation).

    % Получение последней точки списка
    get_last([Head], LastPoint) :- Head = LastPoint.
    get_last([_|Tail], LastPoint) :- get_last(Tail, LastPoint).

    % Строки кубических полиномов (по абсциссам), приравненных значениям соответствующих ординат
    stage_1(Length, [Point|Tail], Equations) :-
        get_last(Tail, LastPoint),
        stage_1_1(Length, Point, Equation1), % Уравнение для первой точки
        stage_1_2(Length, Tail, Equations2), % Уравнения для точек между первой и последней
        stage_1_3(Length, LastPoint, Equation3), % Уравнение для последней точки
        merge([Equation1], Equations2, Eq12), % Слияние уравнений
        merge(Eq12, [Equation3], Equations).

    % Уравнение для первой точки
    stage_1_1(Length, Point, Equation) :-
        ZerosAfterCount = 4 * (Length - 2),
        create_equation_1(Point, 0, ZerosAfterCount, Equation).

    % Уравнения для точек между первой и последней
    stage_1_2(Length, Tail, Equations) :- stage_1_2(Length, Tail, Equations, 1).
    stage_1_2(_, [_], [], _).
    stage_1_2(Length, [Point|Tail], [Equation1, Equation2|Equations], CurrentIndex) :-
        ZerosBeforeCount1 = 4 * (CurrentIndex - 1),
        ZerosAfterCount1 = 4 * (Length - CurrentIndex - 1),
        ZerosBeforeCount2 = 4 * CurrentIndex,
        ZerosAfterCount2 = 4 * (Length - CurrentIndex - 2),
        create_equation_1(Point, ZerosBeforeCount1, ZerosAfterCount1, Equation1),
        create_equation_1(Point, ZerosBeforeCount2, ZerosAfterCount2, Equation2),
        NewCurrentIndex = CurrentIndex + 1,
        stage_1_2(Length, Tail, Equations, NewCurrentIndex).

    % Уравнение для последней точки
    stage_1_3(Length, Point, Equation) :-
        ZerosBeforeCount = 4 * (Length - 2),
        create_equation_1(Point, ZerosBeforeCount, 0, Equation).

    % Строки с приравненными значениями первых производных двух полиномов в общей точке
    stage_2(Length, Points, Equations) :- stage_2(Length, Points, Equations, 0).
    stage_2(Length, [_|Tail], Equations, 0) :- stage_2(Length, Tail, Equations, 1).
    stage_2(_, [_], [], N) :- not(N = 0).
    stage_2(Length, [Point|Tail], [Equation|Equations], CurrentIndex) :-
        ZerosBeforeCount = 4 * (CurrentIndex - 1),
        ZerosAfterCount = 4 * (Length - CurrentIndex - 2),
        create_equation_2(Point, ZerosBeforeCount, ZerosAfterCount, Equation),
        NewCurrentIndex = CurrentIndex + 1,
        stage_2(Length, Tail, Equations, NewCurrentIndex).

    % Строки с приравненными значениями вторых производных двух полиномов в общей точке
    stage_3(Length, Points, Equations) :- stage_3(Length, Points, Equations, 0).
    stage_3(Length, [_|Tail], Equations, 0) :- stage_3(Length, Tail, Equations, 1).
    stage_3(_, [_], [], N) :- not(N = 0).
    stage_3(Length, [Point|Tail], [Equation|Equations], CurrentIndex) :-
        ZerosBeforeCount = 4 * (CurrentIndex - 1),
        ZerosAfterCount = 4 * (Length - CurrentIndex - 2),
        create_equation_3(Point, ZerosBeforeCount, ZerosAfterCount, Equation),
        NewCurrentIndex = CurrentIndex + 1,
        stage_3(Length, Tail, Equations, NewCurrentIndex).

    % Строки с приравненными нулю вторыми производными в граничных точках
    stage_4(Length, [Point|Tail], [Equation1, Equation2]) :-
        First_ZerosAfterCount = 4 * (Length - 2),
        create_equation_4(Point, 0, First_ZerosAfterCount, Equation1),
        get_last(Tail, Last),
        Last_ZerosBeforeCount = 4 * (Length - 2),
        create_equation_4(Last, Last_ZerosBeforeCount, 0, Equation2).

    % Решение системы уравнений методом Гаусса
    gaussian_elimination_with_backsubstitution(Ass, Xs) :-
        gaussian_elimination(Ass, [], ReversedReducedAss), % Прямой ход метода Гаусса
        backsubstitution(ReversedReducedAss, [], Xs). % Обратный ход метода Гаусса
        % Пример:
        %    gaussian_elimination_with_backsubstitution([[1, 2, 3], [2, 2, 3]], X)
        % -> X = [0, 1.5]

    % Прямой ход метода Гаусса
    gaussian_elimination([], Xss, Xss).
    gaussian_elimination(Lower, Upper, Xss) :-
        pivot_row(Lower, PivotRow, OtherRows), % Получаем строку с опорным элементом
        normalized_pivot_row(PivotRow, NormalizedPivotRow), % Приводим строку с опорным элементом к нормальному виду
        normalized_other_rows(OtherRows, NormalizedPivotRow, NewLower), % Приводим другие строки к нормальному виду
        gaussian_elimination(NewLower, [NormalizedPivotRow|Upper], Xss). % Повтор
        % Пример:
        %    gaussian_elimination([[1, 2, 3], [2, 2, 3]], [], X)
        % -> X = [[(3-1*(3/2))/(2-1*(2/2))], [2/2, 3/2]] = [[1.5], [1, 1.5]]

    % Получение опорного элемента
    pivot_row([PivotRow], PivotRow, []).
    pivot_row([[X1|Row1],[X2|Row2]|Rows], PivotRow, [[X2|Row2]|Rest]) :-
        abs(X1) > abs(X2), !, % Ищем максимальный по модулю элемент
        pivot_row([[X1|Row1]|Rows], PivotRow, Rest). % Ищем опорный элемент в первом столбце матрицы
    pivot_row([Row1,Row2|Rows], PivotRow, [Row1|Rest]) :- % Перебираем все строки
        pivot_row([Row2|Rows], PivotRow, Rest). % Для нахождения опорного элемента

    % Приведение строки к нормальному виду
    normalized_pivot_row([Pivot|As], Bs) :-
        not(Pivot = 0.0), % Недопускаем деления на нуль
        normalized_pivot_row_auxillary(As, Pivot, Bs).

    % Приведение остальных строк к нормальному виду
    normalized_other_rows([], _, []).
    normalized_other_rows([[A|As]|Ass], PivotRow, [Bs|Bss]) :-
        normalized_other_row(As, A, PivotRow, Bs),
        normalized_other_rows(Ass, PivotRow, Bss).

    % Вспомогательный предикат для приведения строки к нормальному виду
    normalized_pivot_row_auxillary([], _, []).
    normalized_pivot_row_auxillary([A|As], Pivot, [B|Bs]) :-
        B = A / Pivot, % Делим значение головы списка на значение опорного элемента
        normalized_pivot_row_auxillary(As, Pivot, Bs).

    % Приведение строки к нормальному виду
    normalized_other_row([], _, [], []).
    normalized_other_row([A|As], X, [P|Ps], [B|Bs]) :-
        B = A - X * P, % Вычитаем строки из строки с опорным элементом
        normalized_other_row(As, X, Ps, Bs). % Берем следующие числа в списках

    % Обратный ход метода Гаусса
    backsubstitution([], Xs, Xs).
    backsubstitution([As|Ass], Ys, Xs) :- % Обратная подстановка
        backsubstitution_auxillary(Ys, As, 0, Y),
        backsubstitution(Ass, [Y|Ys], Xs). % Обратная подстановка для оставшейся части
        % Пример:
        %    backsubstitution([[1.5], [1, 1.5]], [], X)
        % -> X = [1.5-(0+(1.5-0)*1), 1.5-0] = [0, 1.5]

    % Вспомогательный предикат для обратного хода метода Гаусса
    backsubstitution_auxillary([], [B], Acc0, Acc) :- % Ищем один неизвестный свободный член в строке
        Acc = B - Acc0. % Вычитаем правую часть
    backsubstitution_auxillary([X|Xs], [A|As], Acc0, Acc) :-
        Acc1 = Acc0 + X * A, % Подставляем предыдущие значения
        backsubstitution_auxillary(Xs, As, Acc1, Acc).

    % Вывод результата
    output(Points, Result) :- 
        write("1. Write to console;"), nl,
        write("2. Write to file."), nl,
        write("Another button to exit"), nl,
        readint(C), % Считывание числа
        output_action(Points, Result, C); % Вывод
        write("The program is complete.").

    % Вывод в диалоговое окно
    output_action(Points, Result, 1) :-
        write_coefficients(Points, Result).

    % Вывод в файл
    output_action(Points, Result, 2) :-
        write("File name: "),
        readln(FileName), % Ввод названия файла
        openwrite(datafile, FileName), % Открытие файла для записи
        writedevice(datafile), % Перенаправление вывода на файл
        write_coefficients(Points, Result), % Запись данных
        closefile(datafile). % Закрытие файла

    % Вывод коэффициентов для каждого отрезка
    write_coefficients([_], []).
    write_coefficients([point(X1, Y1), point(X2, Y2)|PointTail], [A, B, C, D|CoefTail]) :-
        write("From X=", X1, " To X=", X2), nl,
        write("A: ", A), nl,
        write("B: ", B), nl,
        write("C: ", C), nl,
        write("D: ", D), nl, nl,
        write_coefficients([point(X2, Y2)|PointTail], CoefTail).

笔记：

1. 可以在此处找到在线高斯方法计算器（用于验证） 。
2. 可以在此处找到在线三次样条计算器（用于验证目的） 。（上述程序是使用此在线服务测试的。）
3. 三次样条可以很好地处理龙格现象。
4. 上述求三次样条系数的算法很容易理解，但时间和内存都非常复杂。这可能会以某种方式进行优化......
5. Turbo Prolog 2.0 早就过时了，尽管我的大学将它用作“逻辑和函数式编程”学科的主要版本。如果可能（您可以说服老师或其他人），请切换到 Visual Prolog 或其他人。
6. 代码是一天写的，所以各位读者不要严格判断。