Численно решить уравнение с 4-мя неизвестными

[svv]

Спасибо за информацию.

[Amw]

Вообще не очень понятно, чего хочет добиться ТС.
1. Решить задачу и сообщить преподавателю решение
2. Запрограммировать микроконтроллер для системы АСУТП, чтобы он в течение десятилетий подхватывал в реальном времени незатухающий синус и идентифицировал его параметры
3. Просто научиться (раз и навсегда) идентифицировать незатухающий синус...
4. Научиться идентифицировать нелинейные модели по экспериментальным данным.

В зависимости от этого и надо выбрать метод решения

О, ещё забыл:
5. Поразить всех нетрадиционным подходом к решению тривиальной задачи.

[Kevsh]

Amw
Ну... вообще этот синус – погрешность, которую я таким способом решил обработать. Эффективность такого метода обработки, честно говоря, не очень высокая, так как есть цифровая фильтрация, которую я уже для этой задачи применил (и она за куда меньшее время выдаёт результат не хуже). Просто было интересно, что можно для такой задачи получить, используя численные методы.

Также интерес у меня тут есть чисто учебный – я не на математика учусь, поэтому задача для меня эта далеко не элементарная . Я до недавнего момента думал, что МНК – это когда формулу из лекции берешь (уж не буду искать, какую) и рассчитываешь какие-то коэффициенты. Про предложенные выше методы (очень удобные, кстати) через псевдообратную матрицу (о ней я тоже не знал) я и не подозревал. Про то, что есть задачи минимизации функции, я тоже не знал. Чтобы вы понимали, наша программа по численным методам выглядела так: решение СЛАУ парой методов, МНК (два типовых расчёта, где приближать всё нужно было полиномами), интерполяция (методы Ньютона и Лагранжа), численное дифференцирование и интегрирование, численное решение дифференциальных уравнений. То есть я что-то плюс-минус знаю, но я уже задумываюсь о том, чтобы в эту тему углубиться, потому что уже не первый раз нужно что-то подобное, а я туплю.

Кстати, вот вы писали про АСУТП и идентификацию параметров незатухающего синуса. А как бы вы такую задачу решали? Через спектральный анализ?

[Amw]

Kevsh в сообщении #1593228 писал(а):

МНК (два типовых расчёта, где приближать всё нужно было полиномами)

Это как раз тот метод, что здесь Вам предложил svv (Линейное программирование, регрессионный анализ)

Kevsh в сообщении #1593228 писал(а):

...про АСУТП и идентификацию параметров незатухающего синуса. А как бы вы такую задачу решали? Через спектральный анализ?

Опять же всё зависит от конечной цели. Можно спектр использовать для выбора начального приближения, а затем "шлифануть" МНК. АСУТП я упомянул как пример того, что задача непрерывно повторяется и её реализация требует особого подхода. Одноразовая же решается "на коленке".

[Евгений Машеров]

Соображения применительно к "промышленным" (а не "разовым" научным) расчётам, в особенности в рассуждении микроконтроллера. Здесь быстрота вычисления может быть критична, особенно на архитектуре микроконтроллера (в котором может вообще не быть расчёта тригонометрии, и синусы/косинусы надо считать подпрограммой).
1. Перебор по сетке ("сканирование") для единичного расчёта в одномерном случае может быть удовлетворительным решением, но уже в многомерном нарываемся на "проклятие размерности", и он годится разве что для получения начальных приближений, уточняемых оптимизационными алгоритмами (градиентным, сопряжённых направлений и т.п.), конкурируя здесь со "случайным набросом". При многократном же расчёте есть резон использовать какой-либо метод одномерной оптимизации ("золотого сечения", квадратичной аппроксимации и др.), Для ЗС, скажем, число вычислений линейной модели при выбранном b для точности 0.001 начального интервала поисков сократится с 1000 до примерно 15. Это может оправдать усложнение (весьма умеренное) программы (напомню, что метод "золотого сечения" в одномерной оптимизации предполагает, что мы указали отрезок, в котором лежит единственный оптимум, вычисляем значения функции на концах его и в двух точках внутри него, разбивая на три отрезка, и отбрасывая тот, где оптимума заведомо нет, затем повторяя процедуру для сокращённого отрезка; "золотое сечение" здесь задаёт положения двух внутренних точек, что связано не с "матемистикой", а с тем, что после отбрасывания части отрезка у нас уже есть вычисленное значение функции в новом отрезке, так что нужно считать не в двух, а в одной точке). В ЗС отрезок, внутри которого оптимум, на каждое вычисление функции сокращается в 1.618... раза. Несложно реализовать и квадратичную аппроксимацию.
2. Вычисление тригонометрических функций затратно даже и на процессорах общей архитектуры, требуя задержки 50-160 тактов (от модели)
https://www.agner.org/optimize/instruction_tables.pdf
И считать в каждой точке синусы и косинусы может быть весьма накладно. Если же микроконтроллер без специальных функций, то вычисление по подпрограмме может затребовать и тысяч тактов.
Однако в данной задаче шаг по времени постоянен, и можно воспользоваться обычными тригонометрическими соотношениями:
$ \sin\left( \alpha \pm \beta \right)= \sin\alpha \, \cos\beta \pm \cos\alpha \, \sin\beta$
$ \cos\left( \alpha \pm \beta \right)= \cos\alpha \, \cos\beta \mp \sin\alpha \, \sin\beta$
Где $\beta$ соответствует шагу аргумента по времени, а $\alpha$ текущему значению аргумента.
То есть собственно расчёт тригонометрических функций делается однократно для заданного значения параметра b, а затем для каждой точки считается не синус и косинус, а 4 умножения и два сложения.

[Евгений Машеров]

Kevsh в сообщении #1593228 писал(а):

Ну... вообще этот синус – погрешность, которую я таким способом решил обработать.

Вот для погрешности способ едва ли не пессимальный. Он предполагает, что у нас есть синусоида неизвестной, но постоянной частоты, а помеха, искажающая её значения, мала. Для оценки спектра погрешности лучше бы спектральные методы, того же Фурье, использовать. Вот оценивать какие-нибудь периодические процессы (экономические циклы, или там что-то хронобиологическое) может подойти.

Kevsh в сообщении #1593228 писал(а):

про АСУТП и идентификацию параметров незатухающего синуса. А как бы вы такую задачу решали? Через спектральный анализ?

Если частота неизвестна - какие-нибудь авторегрессионные методы оценивания. По два коэффициента на частотную составляющую.
Если известны частоты, которые могут быть (но не все одновременно), скажем, это какой-нибудь протокол связи с частотной модуляцией, наподобие DTMF в телефонии: 1209 Hz, 1336 Hz, 1477 Hz , 1633 Hz, 697 Hz, 770 Hz, 852 Hz, 941 Hz, может быть полезен алгоритм Герцеля или что-то в этом роде.
http://dsp-book.narod.ru/dtmf.pdf

[Kevsh]

Евгений Машеров
Я посмотрел, что такое метод золотого сечения, на метод бисекции похоже. Только мне не очень понятно, как его тут применить. На данный момент мы делали так: заготавливаем кучу вариантов $b$ с каким-то шагом, для каждого из этих вариантов применяем МНК, находим вариант с наименьшей невязкой – готово. Как я понял, вы предлагаете не просто разбивать отрезок, в котором лежит $b$, с каким-то шагом, а находить $b$ как-то быстрее и умнее. Для применения ЗС должна быть возможность вычислить функцию в каждой точке, но ведь при старте программы у нас вообще никаких параметров нет (а не только $b$), то есть применить ЗС не выйдет. Может, сделать начальное приближение всех остальных параметров, для них через ЗС найти $b$, а потом для этого $b$ через МНК найти приближение функции?

Также не очень понятно, как использовать предложенные вами тригонометрические соотношения. А первой точке у меня угол $\alpha$, а потом $\alpha + \beta$, $\alpha + 2\beta$, $\alpha + 3\beta$ и т. д., разве не так? Получается, что при таком подходе можно избежать расчёта двух синусов и двух косинусов, но ещё два синуса и два косинуса придётся рассчитывать (а ещё 2 сложения/вычитания и 4 произведения). Хотя, это, конечно, всё равно быстрее.

[Евгений Машеров]

1. Ну, я бы поступал так: получаем некоторое приближение к b (например, посчитав пересечения нулей), исходя из него задаём отрезок поиска $(b_0-\Delta;b_0+\Delta)$(в случае перебора по сетке нужно то же самое). Обращаемся к процедуре расчёта параметров при заданном b (которая вычисляет последовательности синусов и косинусов и строит линейную регрессию, возвращая меру приближения, например, сумму квадратов отклонений), на стартовом шаге четырежды (для значений b на концах выбранного отрезка и в двух точках "золотого сечения"), на последующих отбрасывается часть отрезка, у нас есть три точки, вычислять надо только в одной точке на каждом шаге, минимизируя меру приближения, пока отрезок не укоротится до нужной точности.
2. Рекуррентно. На первом шаге есть синусы-косинусы альфы и беты, получаются синусы-косинусы альфы плюс беты. На следующем используем полученные $\sin(\alpha+\beta)$ и $\cos(\alpha+\beta)$ так же, как на предыдущем использовались синусы-косинусы альфы, получая $\sin(\alpha+2\beta)$ и $\cos(\alpha+2\beta)$ и так далее.

[Alex Krylov]

Kevsh в сообщении #1593228 писал(а):

Amw
Ну... вообще этот синус – погрешность, которую я таким способом решил обработать. Эффективность такого метода обработки, честно говоря, не очень высокая, так как есть цифровая фильтрация, которую я уже для этой задачи применил (и она за куда меньшее время выдаёт результат не хуже). Просто было интересно, что можно для такой задачи получить, используя численные методы.

Также интерес у меня тут есть чисто учебный – я не на математика учусь, поэтому задача для меня эта далеко не элементарная . Я до недавнего момента думал, что МНК – это когда формулу из лекции берешь (уж не буду искать, какую) и рассчитываешь какие-то коэффициенты. Про предложенные выше методы (очень удобные, кстати) через псевдообратную матрицу (о ней я тоже не знал) я и не подозревал. Про то, что есть задачи минимизации функции, я тоже не знал. Чтобы вы понимали, наша программа по численным методам выглядела так: решение СЛАУ парой методов, МНК (два типовых расчёта, где приближать всё нужно было полиномами), интерполяция (методы Ньютона и Лагранжа), численное дифференцирование и интегрирование, численное решение дифференциальных уравнений. То есть я что-то плюс-минус знаю, но я уже задумываюсь о том, чтобы в эту тему углубиться, потому что уже не первый раз нужно что-то подобное, а я туплю.

Кстати, вот вы писали про АСУТП и идентификацию параметров незатухающего синуса. А как бы вы такую задачу решали? Через спектральный анализ?

Цитата:

Кстати, вот вы писали про АСУТП и идентификацию параметров незатухающего синуса. А как бы вы такую задачу решали? Через спектральный анализ?

Если речь о синусе, имеющем постоянную частоту (или вариация частоты мала на интервале накопления), то БПФ - это ОПТИМАЛЬНЫЙ алгоритм ОБНАРУЖЕНИЯ/ОЦЕНИВАНИЯ. Можно оценить амплитуду, частоту, фазу. Уточнить измерения можно как sinc-интерполяцией (заполнив часть БПФ-выборки нулями справа и выполнив БПФ), так и используя для уточнения частоты локальные алгоритмы (нелинейный МНК, дискриминатор итд), инициализировав их данными (включая первичную оценку частоты) от БПФ.

[Евгений Машеров]

Проблема оценок через БПФ (вообще ДПФ, не только "быстрое") в том, что одновременно оцениваются мощности по всем частотам "гребёнки". То есть оцениваем много, но с ошибкой ("оценки несостоятельны").

[Geen]

и уж лучше не Фурье, а Хартли - хоть памяти вдвое меньше жрёт...

[Евгений Машеров]

Нет пар комплексных чисел - но вещественных вдвое больше. "Что совой об пень, что пнём об сову".

[Kevsh]

Евгений Машеров
svv
Вот, решил использовать метод золотого сечения. У нас уже есть график зависимости $S$ от $b$:

Тут, в принципе, видно, что при более-менее нормальном начальном приближении минимум должен искаться хорошо – главное начальным отрезком попасть в "большую впадину". Начальное приближение возьмём, считая пересечения через ноль, оно у нас достаточно близко к искомому $b$: 0,0942 (при искомом 0,0934). Начальный отрезок я беру как начальное приближение плюс-минус 2 процента от него, в нашем случае это от 0,0923 до 0,0961. Если отметить эти точки на графике, то видно, что в " большую впадину" мы попали, причём с запасом (ну и искомый минимум лежит внутри нашего отрезка).

Вот код, который должен искать набор коэффициентов, обеспечивающих минимальную невязку методом золотого сечения:

Используется синтаксис Matlab M

y=[61.8345252303583; 64.2186685630058; 68.5780907223303; 73.8351943469310; ...
   78.1178746381488; 79.4447599032848; 77.1209584326325; 72.3384843492439; ...
   67.1780874592201; 63.3083342155163; 61.5848686426315; 62.2811015529467; ...
   65.2944676251564; 70.0570903999628; 75.2539157928855; 78.8647097424086; ...
   79.1332309482167; 75.9175990711235; 70.8145100998371; 65.8890541188997; ...
   62.5767534515232; 61.5348068381206; 62.9195242022718; 66.5179434492526; ...
   71.5768723085928; 76.5420755341518; 79.3273502810222; 78.5249447120497; ...
   74.5575769012328; 69.3099594259272; 64.7368823911916; 62.0334551154671; ...
   61.6848623894562; 63.7418912006921; 67.8660696583989; 73.0933449460809; ...
   77.6481344114599; 79.4840538567364; 77.6481344114599; 73.0933449460811; ...
   67.8660696583989; 63.7418912006923; 61.6848623894562; 62.0334551154669; ...
   64.7368823911916; 69.3099594259272; 74.5575769012327; 78.5249447120497; ...
   79.3273502810222; 76.5420755341517; 71.5768723085929; 66.5179434492527; ...
   62.9195242022719; 61.5348068381207; 62.5767534515231; 65.8890541188995; ...
   70.8145100998372; 75.9175990711235; 79.1332309482166; 78.8647097424087; ...
   75.2539157928858; 70.0570903999627; 65.2944676251561; 62.2811015529469; ...
   61.5848686426318; 63.3083342155164; 67.1780874592199; 72.3384843492438; ...
   77.1209584326324; 79.4447599032847; 78.1178746381488; 73.8351943469311; ...
   68.5780907223305; 64.2186685630060 ];
n=size(y,1);
x=(1:n)';

init_approx = find_frequency(x, y);
init_approx_delta = init_approx * 0.02;
a_n = init_approx - init_approx_delta;
b_n = init_approx + init_approx_delta;
c_n = (b_n - a_n) * 0.312 + a_n;
d_n = (b_n - a_n) * 0.618 + a_n;
[answ_c, S_c] = LS(init_approx, x, y);
[answ_d, S_d] = LS(init_approx, x, y);

for counter = (1:10)
    if S_c > S_d
        a_n = c_n;
        c_n = d_n;
        d_n = (b_n - a_n) * 0.618 + a_n;
        %disp(["Выбран d: " string(d_n) ". При a=" string(a_n) ", c=" string(c_n) ", b=" string(b_n)])
        [answ_d, S_d] = LS(d_n, x, y);
    else
        b_n = d_n;
        d_n = c_n;
        c_n = (b_n - a_n) * 0.312 + a_n;
        %disp(["Выбран c: " string(c_n) ". При a=" string(a_n) ", d=" string(d_n) ", b=" string(b_n)])
        [answ_c, S_c] = LS(c_n, x, y);
    end
end

[c, S] = LS(a_n + (b_n - a_n) / 2, x, y);
disp(c)
disp(a_n + (b_n - a_n) / 2)
disp(S)

function [c, S] = LS(b, x, y)
    A=[ones(size(y,1),1), cos(2*pi*b*x), sin(2*pi*b*x)];
    c=pinv(A)*y;
    S=norm(y-A*c);
end

function frequency = find_frequency(x, y)
    counter = 0;
    first_i = 0;
    last_i = 0;
    zero = (min(y) + max(y)) / 2;
    for i = 2:size(x, 1)
        if y(i - 1) > zero && y(i) < zero || y(i - 1) < zero && y(i) > zero
            counter = counter + 1;
            if counter == 1
                first_i = i;
            else
                last_i = i;
            end
        end 
    end
    frequency = (counter - 1) / (last_i - first_i) / 2;
end
 

Данный код на выдаёт $b=0,0927$ при невязке $S=5,9165$. Интересно не только то, что это не искомый коэффициент $b$, но и то, что невязка не та . Если ткнуть в такую точку на приведённом выше графике, то невязка там равна 5,76103, отсюда и проблема. Видимо, я неправильно преобразовал функцию поиска коэффициентов. Вот так было:

Используется синтаксис Matlab M

y=[61.8345252303583; 64.2186685630058; 68.5780907223303; 73.8351943469310; ...
   78.1178746381488; 79.4447599032848; 77.1209584326325; 72.3384843492439; ...
   67.1780874592201; 63.3083342155163; 61.5848686426315; 62.2811015529467; ...
   65.2944676251564; 70.0570903999628; 75.2539157928855; 78.8647097424086; ...
   79.1332309482167; 75.9175990711235; 70.8145100998371; 65.8890541188997; ...
   62.5767534515232; 61.5348068381206; 62.9195242022718; 66.5179434492526; ...
   71.5768723085928; 76.5420755341518; 79.3273502810222; 78.5249447120497; ...
   74.5575769012328; 69.3099594259272; 64.7368823911916; 62.0334551154671; ...
   61.6848623894562; 63.7418912006921; 67.8660696583989; 73.0933449460809; ...
   77.6481344114599; 79.4840538567364; 77.6481344114599; 73.0933449460811; ...
   67.8660696583989; 63.7418912006923; 61.6848623894562; 62.0334551154669; ...
   64.7368823911916; 69.3099594259272; 74.5575769012327; 78.5249447120497; ...
   79.3273502810222; 76.5420755341517; 71.5768723085929; 66.5179434492527; ...
   62.9195242022719; 61.5348068381207; 62.5767534515231; 65.8890541188995; ...
   70.8145100998372; 75.9175990711235; 79.1332309482166; 78.8647097424087; ...
   75.2539157928858; 70.0570903999627; 65.2944676251561; 62.2811015529469; ...
   61.5848686426318; 63.3083342155164; 67.1780874592199; 72.3384843492438; ...
   77.1209584326324; 79.4447599032847; 78.1178746381488; 73.8351943469311; ...
   68.5780907223305; 64.2186685630060 ];
n=size(y,1);
x=(1:n)';

b=0.02:0.0001:0.16;
N=size(b,2);
cc=zeros(3,N);
S=zeros(N);
for j=1:N
   A=[ones(n,1), cos(2*pi*b(j)*x), sin(2*pi*b(j)*x)];
   c=pinv(A)*y;
   S(j)=norm(y-A*c);
   cc(:,j)=c;
end

[Smin,Ind]=min(S);
b_=b(Ind(1))
c_=cc(:,Ind(1))
 

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

Используется синтаксис Matlab M

function [c, S] = LS(b, x, y)
    A=[ones(size(y,1),1), cos(2*pi*b*x), sin(2*pi*b*x)];
    c=pinv(A)*y;
    S=norm(y-A*c);
end
 

Ответ такая функция выдаёт не совсем верный. Хотя я ничего не менял – x тот же, y тот же. Не понимаю

-- 11.05.2023, 00:53 --

Ещё заметил – если в "изначальном" коде в 10 раз уменьшить шаг поиска $b$, то найденный коэффициент $b$ не меняется, невязка не меняется, а остальные коэффициенты меняются

[Alex Krylov]

Евгений Машеров в сообщении #1593359 писал(а):

Проблема оценок через БПФ (вообще ДПФ, не только "быстрое") в том, что одновременно оцениваются мощности по всем частотам "гребёнки". То есть оцениваем много, но с ошибкой ("оценки несостоятельны").

Вы, как я понимаю, под оценкой амплитуды посредством БПФ/ДПФ имеете в виду вот это: $\sqrt{\operatorname{Re}(...)^2+\operatorname{Im}(...)^2}$ и вас смущает неравенство матожидания этой величины амплитуде... Или вас смущают боковые лепестки фильтров гребенки?

Касательно 1-го момента... Мы ж не знаем реальных условий (SNR, требований по точности). Мож там высокий SNR, не предъявляются высокие требования по точности...

Касательно боковых лепестков... В радиолокации это существенно, и заранее заложена возможность выполнения аподизации. А в дешевых связных или навигационных приемниках с 3-битным сигналом такая возможность в принципе не заложена. Т.е. и в этом случае мы не знаем реалий....

Касательно измерения мощности сигнала/шума/SNR. Вот тут есть обзорная статья (можно скачать pdf), где сравниваются 5 классических алгоритмов, используемых в индустрии для этих целей:
GNSS Solutions. Are Carrier-to-Noise Algorithms Equivalent in All Situations?
https://insidegnss.com/are-carrier-to-noise-algorithms-equivalent-in-all-situations/

[svv]

Kevsh в сообщении #1593393 писал(а):

Ещё заметил – если в "изначальном" коде в 10 раз уменьшить шаг поиска $b$, то найденный коэффициент $b$ не меняется, невязка не меняется, а остальные коэффициенты меняются

Мы нашли с некоторой точностью точку экстремума $b=b_0$ невязки $S(b)$. В этой точке $\frac{dS}{db}=0$. Поэтому в малой окрестности точки $b_0$ справедливо приближение
$S(b)\approx S(b_0)+\frac 1 2 S''(b_0)\;(b-b_0)^2$
То есть $S(b)$ в окрестности $b_0$ меняется квадратично. То есть медленно.

А у коэффициентов в точке $b_0$ нет экстремума. Их первая производная там не равна нулю, и они меняются относительно быстро.

MATLAB по умолчанию выводит вещественные числа с четырьмя знаками после точки. (Кстати, кто знает, как это изменить?) Для коэффициентов этой точности хватает, чтобы заметить разницу при большем и меньшем шаге, для медленно меняющейся невязки — нет.