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

[Alex Krylov]

svv в сообщении #1593439 писал(а):

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

Цитата:

Кстати, кто знает, как это изменить?

format('long');

[svv]

Alex Krylov, большое спасибо
Kevsh
Ну вот, с большей точностью значения $S,c(1),c(2),c(3)$ равны
для шага $0.0001$
0.093400000000000
70.000447437561220
-8.537710132040177
-2.723221091009695
для шага $0.00001$
0.093360000000000
69.996794398696181
-8.564958407287474
-2.642092797448871

[svv]

Kevsh
Хочу показать один интересный момент. Взгляните на картинку

Кружочки — заданные по условию опорные точки $(x_i,y_i)$.

Обратите внимание, что кружочки, близкие к максимумам и минимумам нашей аппроксимирующей синусоиды, находятся заметно выше "оптимальной" синусоиды. Возникает вопрос: почему же наша синусоида не подвинулась чуть выше (за счёт увеличения $d=c_1$), чтобы уменьшить невязку?

Ответ. В Вашем сигнале $y_i$ присутствует довольно заметная вторая гармоника. И её фаза такова, что её максимумы совпадают с максимумами и минимумами основной синусоиды, а её минимумы приходятся на "склоны" основной синусоиды. А на склонах разница мало заметна на графике. Вот и кажется, что синусоиду стоило бы подвинуть вверх.
Вот как это выглядит: красная кривая — исходный сигнал $y_i$, синяя — остаток (исходный сигнал минус наша аппроксимация, т.е. $y_i-f(x_i)$). Постоянные составляющие выброшены. Видно, что львиная доля остатка приходится на вторую гармонику.

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

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

[B@R5uk]

svv, матлаб имеет встроенные функции для поиска экстремума: одномерная fminbnd на отрезке и многомерная fminsearch, использующая симплекс-метод Нелдера-Мида. Работают гораздо быстрее, чем поиск в лоб перебором.

-- 12.05.2023, 14:14 --

Алсо, высшие гармоники в этом сигнале однозначно присутствуют. Отбрасывание их в модели сигнала приведёт к системной ошибке расчёта частоты.

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

>> coords
coords =
  70.000000022070964
  -8.571272756147744
   2.587819742489675
   0.420959973399330
  -0.279685155313799
  -0.015509248841234
   0.018747413635993
   0.000348554161859
  -0.000829269933344
>> w
w =
   0.586430628851210
>>

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

clc
clearvars
format compact
format long

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
];
y = y';
dsize = numel (y);
x = (1 : dsize)' - 1;

func_base = @ (p) [ones(numel (p), 1), cos(p), sin(p), cos(2 * p), sin(2 * p), cos(3 * p), sin(3 * p), cos(4 * p), sin(4 * p)];
func_coords = @ (p, y_) func_base (p) \ y_;
func_target = @ (base, y_) sum ((y_ - base * (base \ y_)) .^ 2);

w = fminbnd (@ (w) func_target (func_base (w * x), y), 0.5, 0.7, optimset ('TolX', 1e-12));
coords = func_coords (w * x, y);
xx = (0 : 0.02 : dsize)';
yy = func_base (w * xx) * coords;
y_err = y - func_base (w * x) * coords;

subplot (3, 1, 1 : 2)
plot (xx, yy, 'k')
hold on
plot (x, y, 'o')
hold off
grid on

subplot (3, 1, 3)
plot (x, y_err, 'o')
grid on
 

[svv]

B@R5uk, большое спасибо, но в задаче есть специфика:

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

решать это всё планируется на микроконтроллере

Поэтому я пытаюсь писать код, прозрачный с точки зрения того, что там делается:

svv в сообщении #1593176 писал(а):

Чтобы избавиться от слишком высокоуровневых функций и приблизить код MATLAB к тому, что будет исполняться на микроконтроллере, перепишем цикл так:

Кроме того, я изначально отводил для МНК только следующую роль:

svv в сообщении #1592659 писал(а):

Так Ваша задача сводится к задаче минимизации функции одной переменной $b$.

А эта минимизация при желании может быть реализована разными методами. Но это уже будет минимизация функции одной переменной, с вот таким минимумом:

С остальными тремя переменными расправляется линейный МНК (после фиксации переменной $b$ получается линейная комбинация трёх известных функций с неизвестными коэффициентами), и он делает это быстро и хорошо.

[B@R5uk]

svv, разумеется здесь надо минимизировать только по частоте. Это даже в моём коде выше реализовано (функция fminbnd, а не fminsearch).

[Kevsh]

B@R5uk
Я посмотрел ваш код и решил похожим образом "спасти" свой. Вроде сделал то же самое, но таких отличных результатов не получил. Вот код для одной гармоники:

B_COEF=fminbnd(@LS1, a_n, b_n); @a_n и b_n – начальные приближения, которые точно выбраны верно 
A=[ones(size(y,1),1), cos(2*pi*B_COEF*x), sin(2*pi*B_COEF*x)];
COEFS=pinv(A)*y

function S = LS1(b)
    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)';
    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

Получаем следующее:

COEFS =

  69.997493484633836
  -8.559819480269374
  -2.657659629898907

Хорошо, добавим ещё одну гармонику:

B_COEF2=fminbnd(@LS2, a_n, b_n)
A2=[ones(size(y,1),1), cos(2*pi*B_COEF2*x), sin(2*pi*B_COEF2*x),...
    cos(2*pi*B_COEF2*x*2), sin(2*pi*B_COEF2*x*2)];
COEFS2=pinv(A2)*y

function S = LS2(b)
    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)';
    A=[ones(size(y,1),1), cos(2*pi*b*x), sin(2*pi*b*x), cos(2*pi*b*x*2), sin(2*pi*b*x*2)];
    c=pinv(A)*y;
    S=norm(y-A*c); 
end

Получим следующее:

COEFS2 =

  69.999382251747747
  -8.577036036723200
  -2.568021866649669
   0.423315666248895
   0.278063775718108

Точность повысилась (напомню, что постоянная составляющая в идеальном случае должна равняться 70). Добавим третью гармонику:

B_COEF3=fminbnd(@LS3, a_n, b_n)
A3=[ones(size(y,1),1), cos(2*pi*B_COEF3*x), sin(2*pi*B_COEF3*x), ...
    cos(2*pi*B_COEF3*x*2), sin(2*pi*B_COEF3*x*2), cos(2*pi*B_COEF3*x*3), sin(2*pi*B_COEF3*x*3)];
COEFS3=pinv(A3)*y

function S = LS3(b)
    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)';
    A=[ones(size(y,1),1), cos(2*pi*b*x), sin(2*pi*b*x), cos(2*pi*b*x*2), sin(2*pi*b*x*2), ...
        cos(2*pi*b*x*3), sin(2*pi*b*x*3)];
    c=pinv(A)*y;
    S=norm(y-A*c); 
end

Получим вывод:

COEFS3 =

  69.999189367018261
  -8.577273863261093
  -2.568672462247298
   0.422812571305099
   0.278117683581578
  -0.015454583668819
  -0.018555204536330

Результат изменился, но не в нашу пользу. Не понимаю, почему такое происходит. Также не очень ясно, какой вы в своём графике (где синусоида проходит точно по точкам) выводите сигнал. Если полигармонический, то я не понимаю, как такое возможно – у меня выходит совсем не синусоида. Если просто используете основную гармонику, то тоже не очень ясно. Изначально отклонение от идеальной постоянной составляющей было 3 тысячных. У вас, конечно, в тысячи раз меньше получилось, но ведь на графике, который был приведён svv, видно, что синусоиду нужно сдвинуть где-то на 0.5

[B@R5uk]

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

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

Сумма найденных компонент изображена на верхнем графике (сумма синусов — поэтому он хорошо проходит по экспериментальным точкам). Отклонение экспериментальных точек от модели — на нижнем (обратите внимание на множитель десять в минус пятой для вертикальной шкалы).

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

Результат изменился, но не в нашу пользу. Не понимаю, почему такое происходит.

Вы забыли указать точность для функции fminbnd, с которой она ищет минимум (optimset ('TolX', 1e-12)). По умолчанию точность по аргументу вроде десять в минус четвёртой, что явно меньше, чем у меня в коде.

Вообще, вопрос точности в такой ситуации очень расплывчатый. Ясно, что чем больше гармоник в модели и чем больше периодов сигнала в выборке, тем точнее результат. Однако, каждый конкретный результат сильно зависит в том числе и от фазы сигнала (если повезёт можно получить идеальную точность даже с одной гармоникой). Тот факт, что коэффициент нелинейности синуса очень мал, гарантирует бессмысленность учёта большого числа старших гармоник. Просто потому, что точность вычисления коэффициентов в обратной матрице будет ограничена точностью арифметики с плавающей точкой. В матлабе точность 52 двоичных разряда или 16 десятичных. На МК такое будет тяжело вытянуть (хотя всё зависит от МК, конечно). Другое дело, конечно, если вы сможите вашу pinv(A2) выразить явной формулой. Во всяком случае для одной гармоники суммы квадратов синусов/косинусов очень хорошо сворачиваются.