Как увеличить точность расчетов до long double - Matlab - Страница 6

@vittiv · 05.09.2020, 19:51

Студворк — интернет-сервис помощи студентам

Я для себя это понимаю так, что 25 мерная фигура очень вытянута вниз и в месте, где находится ее минимум она очень острая, как игла... и угол на кончике этой иглы очень не стабильный. Но при этом это все-равно около минимума... игла то острая... а значение угла касательной к вершине иглы не так иж и важно.

@Alexis · 05.09.2020, 20:01

Самое главное, предельная ошибка итеративного метода:

$https://www.cyberforum.ru/cgi-bin/latex.cgi?||\Delta x||=\frac{||x||}{1-||B||}eps$

т.е. эти методы намного точнее прямых методов, и QR в том числе.

На всякий случай привожу последний вариант алгоритма.
В нём используются максимальные, что позволяет определить максимальную ошибку решения:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
function [x, err]=Iter_Solve(A, b, x, n_iter)
% уравнение Ax=b;
% n_iter - число итераций алгоритма
% err - максимальная ошибка для любого элемента решения x
 
% проверка сходимости и определение оптимального шага
s=eig(A);
if (isreal(s) && min(s)>0) 
     step=2/(max(s)+min(s));   % шаг 
     N=size(A,1);
 
     g=b*step;
     B=eye(N)-A*step;
     nb=norm(B,inf);
 
    % итеративное уточнение
    for ii=1:n_iter;
            xold=x;
            x=B*x+g;
    end
     err=(nb*norm(x-xold,inf)+norm(x,inf)*eps)/(1-nb);  % максимальная ошибка x
else
     error('Алгоритм не сходится. Матрица системы вырождена.');
end
end

его можно повторять до тех пор, пока err не перестанет уменьшаться.

Добавлено через 6 минут

Сообщение от vittiv

она очень острая, как игла...

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

@vittiv · 05.09.2020, 20:07

Сообщение от Alexis333

последний вариант алгоритма.

Попробую, результат сообщу...

Но вот, что интересное я обнаружил...

Вот мой способ решения СЛАУ, который Вы раскритиковали...

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
[szyk,szxk]=size(K);
 
xo=double(zeros(szyk,szxk));
 
for i=1:szxk
 
    St=S;
    Kt=K(:,i);
    for j=1:szyk
        St(j,j)=St(j,j)+St(j,j)*dltm/100;
    end
 
    x=St\Kt;
 
    for its=1:10
 
        St=S;
        for j=1:szyk
            St(j,j)=St(j,j)+St(j,j)*dltm*(its^2/100)/abs(x(j));
        end
 
        x=St\Kt;
 
    end
 
    w=Kt-St*x;
    x=x+St\w;
 
    xo(:,i)=double(x);
 
end

дает лучший результат аппроксимации. Точное прямое решение по Гаусу с нормой невязки = 0.0002 дало точность RSME = 8.45191, а мой метод дал норму невязки = 250 но точность RMSE стала = 8.45177.

Как это можно объяснить? Возможно угол так колбасит возле минимума, что легко можно проскочить истинный минимум и выдать не минимум за минимум.

Кстати в этом методе, я как раз прибавляю к диагонале матрицы S только не константу... а вектор пропорциональный найденному промежуточному значению решению СЛАУ.

Добавлено через 58 секунд

Сообщение от Alexis333

многомерной параболой

тонкая по всем измерениям многомерная игла...

@Alexis · 05.09.2020, 20:18

Сообщение от vittiv

Как это можно объяснить?

Я объ этом уже писал раньше, возможно Вы не внимательно читали. Это объясняется тем, что между внутренней невязкой СЛАУ и точностью решения задачи МНК нет никакой связи. Невязка может увеличиваться, а точность только расти, что собственно Вы и наблюдаете.

Итеративный алгоритм, который я выложил, по идее, должен дать ещё меньшее значение rmse. Проверьте.

@vittiv · 05.09.2020, 20:23

Сообщение от Alexis333

Невязка может увеличиваться, а точность только расти,

Про внутреннюю невязку СЛАУ я не понял... это на уровне шаманства? )... мне более симпатичен вариант с иглой.

Как управлять этим процессом?

Ваша программа управляет этим процессом?

@Alexis · 05.09.2020, 20:42

Ну Вы попробуйте - и узнаете.

Добавлено через 2 минуты

Сообщение от vittiv

Про внутреннюю невязку СЛАУ я не понял... это на уровне шаманства? )

в каком смысле? если Вы о попытках её вычисления, то я склоняюсь к мысли, что да. Её вообще видимо не нужно было вычислять.

Добавлено через 13 минут

Сообщение от vittiv

RMSE стала = 8.45177

нужно было сразу смотреть по RMSE, т.к. именно она является целевой функцией, а не какая то сомнительная внутренняя невязка СЛАУ

@vittiv · 05.09.2020, 21:47

Сообщение от Alexis333

его можно повторять до тех пор, пока err не перестанет уменьшаться.

Чтобы не разбираться в Вашем коде и по причине того, что не сложно... Можно попросить Вас дописать программу, чтобы в программе показывалось err в теле цикла. И чтобы можно было задавать цель по err и максимальное число циклов.

Результат:

err = -5e-13
RMSE = 8.45191
норма невязки = 0.00292

Т.е. ничего не изменилось

Добавлено через 7 минут

Сообщение от Alexis333

его можно повторять до тех пор, пока err не перестанет уменьшаться.

Чтобы не разбираться в Вашем коде и по причине того, что не сложно... Можно попросить Вас дописать программу, чтобы в программе показывалось err в теле цикла. И чтобы можно было задавать цель по err и максимальное число циклов.

Результат:

err = -5e-13
RMSE = 8.45191
норма невязки = 0.00292... более точный RMSE результат был при норме невязки = ~250 - т.е. более 85 тысяч раз больше...

Т.е. ничего не изменилось

Все-таки непонятно, что является целевой функцией... понятно, что можно в цикле RMSE считать, но это не реально долго будет...

@Alexis · 05.09.2020, 21:51

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
function x =Iter_Solve(A, b, x, n_iter)
% уравнение Ax=b;
% n_iter - число итераций алгоритма
 
% проверка сходимости и определение оптимального шага
s=eig(A);
if (isreal(s) && min(s)>0) 
     step=2/(max(s)+min(s));   % шаг 
     N=size(A,1);
 
     g=b*step;
     B=eye(N)-A*step;
     nb=norm(B,inf);
 
    % итеративное уточнение
    for ii=1:n_iter;
            xold=x;
            x=B*x+g;
         err=(nb*norm(x-xold,inf)+norm(x,inf)*eps)/(1-nb);  % максимальная ошибка x
        disp(err);
    end
else
     error('Алгоритм не сходится. Матрица системы вырождена.');
end
end

Добавлено через 2 минуты

Сообщение от vittiv

err = -5e-13

она точно отрицательная?

@vittiv · 05.09.2020, 21:53

Сообщение от Alexis333

она точно отрицательная?

Да... все (для каждого столбца К)

@Alexis · 05.09.2020, 22:04

Но это невозможно, ведь норма всегда положительна, может какое нибудь переполнение.
Посмотрите в динамике, когда err становится отрицательной, и какое перед этим имеет значение

@vittiv · 05.09.2020, 22:06

Matlab M
1
err=(nb*norm(x-xold,inf)+norm(x,inf)*eps)/(1-nb)

если норма nb больше 1 то легко....

@Alexis · 05.09.2020, 22:12

Сообщение от vittiv

если норма nb больше 1 то легко....

нет, norm(B)<1 это условие сходимости, и оно выполняется, иначе там вообще ни о какой точности говорить бы не пришлось

@vittiv · 05.09.2020, 22:28

-4.521857593879361e-13

-4.521857593879361e-13

-4.521857593879361e-13

-4.521857593879361e-13

-4.521857593879361e-13

-4.521857593879360e-13

-4.521857593879360e-13

-4.521857593879360e-13

-4.521857593879360e-13

-4.521857593879360e-13

что-то не сильно меняется

@Alexis · 05.09.2020, 22:49

Что толку от изменений, если она отрицательная

Добавлено через 10 минут
Я разобрался. Действительно, максимальная норма B больше единицы, но евклидова меньше единицы. Поэтому алгоритм сходится, а ошибка получается отрицательной.
Посчитайте ещё максимальные нормы столбцов K

Matlab M
1
norm(K(:,i),inf)

@vittiv · 05.09.2020, 23:06

[nb norm(x-xold,inf) norm(x,inf) eps]

[4.9415227247783 0.0000000000004 207.2060849426498 0.0000000000000]

err=-5.283449373702147e-13

положительные компоненты но сумма меньше 0

Добавлено через 6 минут
ой, да 4.94>1

@Alexis · 05.09.2020, 23:13

нужно посчитать норму столбцов K, чтобы оценить относительную погрешность

Matlab M
1
norm(K(:,i),inf)

@vittiv · 05.09.2020, 23:15

norm(K(:,i),inf)=~2е7

@Alexis · 05.09.2020, 23:21

а норма самого x какая?

Matlab M
1
norm(x,inf)

@vittiv · 05.09.2020, 23:22

207.2060849426498

@Alexis · 05.09.2020, 23:28

Верхний предел относительной погрешности. В евклидовой норме она будет ещё меньше:
5.283449373702147e-13/207.2=2.5e-15
точнее для double просто уже некуда.

Чтобы проследить его работу, задайте начальное приближение x нулями, тогда будет видно как меняется ошибка. Видимо, начальное решение и так довольно точное, поэтому изменений и нет.

Новые блоги и статьи Все статьи Все блоги /
Был там один разговор по поводу свободы в материальном мире. kumehtar 19.08.2026 Суть: рассматривается живое существо, оказавшееся внутри довольно странной системы (этого мира) и пытающееся обустроить в ней свой кусок пространства. Жизнь действительно предъявляет каждому. . .	Когда логика программы не спасает от человеческих ошибок Maks 18.08.2026 В последнее время всё чаще и чаще сталкиваюсь с таким явлением, как абсолютная невнимательность (или глупость) пользователей. Проявляется это чаще всего на работе в коллективе. Допустим, человек с. . .	Лето уходит kumehtar 17.08.2026	Мысли в слух kumehtar 17.08.2026 Забавно, насколько сейчас стала доступна информация. Например о магии, духовном развитии, медитациях, и других подобных направлениях, ранее зачастую тайных, передаваемых от учителя к ученику. Хотя. . .
Перемещение строк из ТЧ в другой документ с учетом текущего пробега Maks 17.08.2026 Реализация из решения ниже выполнена на примере нетипового документа "Автозапчасти", с ТЧ "Шины". За основу взят алгоритм отсюда: https:/ / www. cyberforum. ru/ blogs/ 359708/ 10838. html Задача: . . .	Саморегулирующийся социальный контракт для сервера cross-section. Hrethgir 14.08.2026 С кодом конечно таких глубоких размышлений пока не было, впрочем я уже привык к алгоритмизации. Суть предмета записи: снова в диалоге с нейросетью (я взял пока себе ник для учётки админа - Rector). . . .	Часы электронные Uhbif79 12.08.2026 Выкладываю программу часов. Программа позволяет: 1. Использовать системное время и дату, 2. Есть возможность вводить время и дату вручную. 3. Реализованы 2 будильника: начало и конец рабочего дня. . . .	Часы с будильником на основе класса QLCDNumber Uhbif79 12.08.2026 Всем добрый день, выкладываю программу часов с будильником на основе класса QLCDNumber. Здесь я пробовал самостоятельно создавал классы, впервые столкнулся с видимостью переменной одного класса из. . .

@vittiv 9 / 1 / 0 Регистрация: 08.09.2019 Сообщений: 343
	05.09.2020, 19:51
	Я для себя это понимаю так, что 25 мерная фигура очень вытянута вниз и в месте, где находится ее минимум она очень острая, как игла... и угол на кончике этой иглы очень не стабильный. Но при этом это все-равно около минимума... игла то острая... а значение угла касательной к вершине иглы не так иж и важно. 0

@Alexis 495 / 542 / 131 Регистрация: 31.10.2016 Сообщений: 4,388
	05.09.2020, 22:04
	Но это невозможно, ведь норма всегда положительна, может какое нибудь переполнение. Посмотрите в динамике, когда err становится отрицательной, и какое перед этим имеет значение 0

@vittiv 9 / 1 / 0 Регистрация: 08.09.2019 Сообщений: 343
	05.09.2020, 22:28
	-4.521857593879361e-13 -4.521857593879361e-13 -4.521857593879361e-13 -4.521857593879361e-13 -4.521857593879361e-13 -4.521857593879360e-13 -4.521857593879360e-13 -4.521857593879360e-13 -4.521857593879360e-13 -4.521857593879360e-13 что-то не сильно меняется 0

@vittiv 9 / 1 / 0 Регистрация: 08.09.2019 Сообщений: 343
	05.09.2020, 23:06
	[nb norm(x-xold,inf) norm(x,inf) eps] [4.9415227247783 0.0000000000004 207.2060849426498 0.0000000000000] err=-5.283449373702147e-13 положительные компоненты но сумма меньше 0 Добавлено через 6 минут ой, да 4.94>1 0

@vittiv 9 / 1 / 0 Регистрация: 08.09.2019 Сообщений: 343
	05.09.2020, 23:15
	norm(K(:,i),inf)=~2е7 0

@vittiv 9 / 1 / 0 Регистрация: 08.09.2019 Сообщений: 343
	05.09.2020, 23:22
	207.2060849426498 0

@Alexis 495 / 542 / 131 Регистрация: 31.10.2016 Сообщений: 4,388
	05.09.2020, 23:28
	Верхний предел относительной погрешности. В евклидовой норме она будет ещё меньше: 5.283449373702147e-13/207.2=2.5e-15 точнее для double просто уже некуда. Чтобы проследить его работу, задайте начальное приближение x нулями, тогда будет видно как меняется ошибка. Видимо, начальное решение и так довольно точное, поэтому изменений и нет. 0