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

[Kevsh]

svv
Да, они равноотстоящие. Это натуральные числа от 1 до 74 включительно.

[Amw]

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

что значит "отсканировать"

Это значит задать область изменения переменных, выбрать шаг их изменения и вычислить функцию в каждой точке. Если переменных четыре и каждая принимает 100 значений (например от 1 до 100 с шагом 1) то всего будет $100^4$ значений функции. Выберите из них минимальное - это и будет первым приближением решения Вашей задачи.

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

что за локальные минимумы

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

[Amw]

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

значения $d,p,q$, обеспечивающие минимум, (ну, и сам минимум) находятся методом наименьших квадратов.

Вот уж "метод наименьших квадратов" тут совсем ни при чем.

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

Amw в сообщении #1592748 писал(а):

Вот уж "метод наименьших квадратов" тут совсем ни при чем.

0. Это именно задача метода наименьших квадратов. И все предложенные подходы - его варианты.
1. Метод Марквардта - нелинейный МНК, полезный при "умеренной нелинейности", итеративно приближающий решение линейными МНК, в которых матрица регрессоров - значения производных нелинейной функции по параметрам при текущем значении вектора коэффициентов, регрессанд - разность между наблюдаемыми значениями зависимой переменной и её оценкой при данном значении вектора коэффициентов, а полученная оценка коэффициентов данной линейной регрессии - поправка к вектору коэффициентов оцениваемой нелинейной модели.
2. В данной модели есть упрощение, здесь лишь два параметра существенно нелинейны, и можно оставить лишь один, воспользовавшись $A\sin(\omega t+\varphi)=A_1\sin(\omega t)+A_2\cos(\omega t)$
В частности, можно организовать перебор по значениям омеги, оценивая прочие параметры, как линейную регрессию.
3. Для уточнения можно рассматривать задачу, как одномерную оптимизацию по омеге, используя в качестве целевой функции сумму квадратов линейной регрессии при заданной омеге. Ну, скажем, "золотым сечением".
4. Оценку частотного параметра омега можно получить, рассмотрев ряд из разностей. При этом первые разности избавятся от постоянного слагаемого, но синус превратится в косинус, вторые вернут синус, третьи опять косинус и т.п. Причём коэффициент будет меняться в зависимости от омеги, выражения получите сами. Вероятно, удобно сопоставлять первые и третьи разности.
5. Простой, но работающий метод оценки - считать переходы через 0 (поскольку есть свободный член - через среднее значение). Расстояние между переходами - полпериода.

[Amw]

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

Это именно задача метода наименьших квадратов.

С чего бы - здесь уже задана целевая функция в явном виде.
А если бы она была не

$S(a,b,c,d)=\sum\limits_n(y_n-a\sin{(2\pi bx_n+c)-d)^2}$,
а
$S(a,b,c,d)=\sum\limits_n(y_n-a\sin{(2\pi bx_n+c)-d)^3}$

тогда бы это был "метод наименьших кубов?

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

Метод Марквардта - нелинейный МНК

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

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

[svv]

Amw в сообщении #1592748 писал(а):

Вот уж "метод наименьших квадратов" тут совсем ни при чем.

Для пущей интриги: посмотрите, я в своём сообщении утверждаю, что применим не просто МНК, а ещё и в линейном варианте. При условии, что $b$ известно, он "убивает" остальные три искомых параметра.

[Amw]

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

посмотрите, я в своём сообщении
утверждаю, что применим не просто МНК, а ещё и в линейном варианте

Вот если бы у ТС задача была сформулирована как:
есть модель $y=a\sin(2\pi bx+c)+d$ и $n$ пар экспериментальных значений $x$ и $y$ и следует найти параметры $a,b,c,d$
тогда бы следовало применив МНК, составить функцию
$S(a,b,c,d)=\sum\limits_n(y_n-a\sin{(2\pi bx_n+c)-d)^2}$
и решить задачу нелинейного программирования - найти её минимум.
Или воспользоваться Вашим советом (который годится не всегда) и свести задачу к линейной - получить уравнение регрессии.
Но задача была другой

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

Мне нужно найти минимум данной функции:

-- 06.05.2023, 17:20 --

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

применим не просто МНК, а ещё и в линейном варианте.

Причем тысячу раз при разных заданных значениях $b$

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

Предложенная ТС постановка задачи - найти решение методом наименьших квадратов. поскольку:
а. Невязки квадрируются.
б. Их сумма минимизируется.
Если бы у него были бы кубы - можно было бы рассуждать о "методе наименьших кубов", если он кому-то нужен. Но там квадраты.

-- 06 май 2023, 18:27 --

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

Мне нужно найти минимум данной функции:
$S(a,b,c,d)=\sum\limits_n(y_n-a\sin{(2\pi bx_n+c)-d)^2}$

[Amw]

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

Невязки квадрируются.

Я думаю, что преподаватель хотел научить ТС находить минимум нелинейной функции, а взял её в таком виде потому, что просто знает ответ.

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

Я так не думаю, но не вижу надобности угадывать. Если ТС захочет - расскажет обстоятельства постановки задачи.

[Kevsh]

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

-- 06.05.2023, 21:07 --

Евгений Машеров
Ещё стоит заметить, что функция точно имеет указанный мной вид, это доказано аналитически.

[Amw]

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

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

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

[svv]

Amw в сообщении #1592767 писал(а):

Причем тысячу раз при разных заданных значениях $b$

Да, именно это и предлагается.

Давайте приблизительно оценим трудоёмкость процесса, заодно опишу его подробнее для ТС. Вот выбрали $b$. У нас есть $m=3$ базисных функции: $u_1(x)=1, u_2(x)=\cos 2\pi bx, u_3=\sin 2\pi bx$. Нам нужно подобрать коэффициенты $c_1,c_2,c_3$ так, чтобы минимизировать невязку (мне самому не нравится этот термин, но он самый лаконичный)
$\sum\limits_{i=1}^n (c_1u_1(x_i)+c_2u_2(x_i)+c_3u_3(x_i)-y_i)^2$
Здесь $n=74$ (ну, круглым счётом $100$).

Пусть $A$ — матрица $n\times m$ с элементами $a_{ik}=u_k(x_i)$, где $i=1...n, k=1...m$.
Вектор $c$ имеет размер $m$, вектор $y$ — размер $n$.
Согласно теории, минимальную невязку обеспечивает вектор $c=(A^\top A)^{-1}A^\top y$. Сама невязка равна
$(Ac-y)^\top (Ac-y)=y^\top y-y^\top Ac$
Но поскольку вектор $y^\top y$ не зависит от $b$, для каждого $b$ достаточно вычислять
$y^\top Ac=(A^\top y)^\top(A^\top A)^{-1}(A^\top y)$

Что сюда входит:
$\bullet$ вычисление $n$ синусов и $n$ косинусов, составляющих матрицу $A$;
$\bullet$ вычисление симметричной матрицы $A^\top A$ размера $3\times 3$ — это $6n$ умножений;
$\bullet$ вычисление вектора $A^\top y$ — ещё $3n$ умножений;
$\bullet$ обращение матрицы $A^\top A$ — чуть больше $20$ умножений;
$\bullet$ перемножение уже вычисленных $(A^\top y)^\top,(A^\top A)^{-1},(A^\top y)$ размеров $1\times 3,\,3\times 3,\,3\times 1$ — всего $12$ умножений.
Итого: $150$ триг.функций и не более $1000$ умножений.
Для всех $1000$ значений $b$ потребуется $150000$ триг.функций и миллион умножений.

Вы при сканировании берёте всего 100 различных значений $b$, то есть (почему-то) в 10 раз меньше, чем предложили мне. Ну да ладно. Интересно сделать аналогичный подсчёт для Вашего предложения

Amw в сообщении #1592520 писал(а):

Предварительно отсканировать $100^4$

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

В общем, мне дело представляется так:
1. Это МНК. Ipso facto того, что минимизируется сумма квадратов невязок. Было бы другое определение целевой функции - был бы не МНк.
2. В линейном случае МНК имеет грандиозные преимущества, в силу того, что производные квадратичной функции линейны. И если сама модель описывается линейной зависимостью, то решение сводится к решению системы линейных уравнений, что не только вычислительно просто, но и гарантирует единственность ответа.
3. Если модель описана нелинейной функцией - преимущества частично сохраняются, но лишь частично. Эксплуатирует эти преимущества метод Левенберга-Марквардта, сводя задачу к серии линейных, полученных разложением в ряд Тейлора относительно значений параметров, и выбором лишь линейных членов ряда. Единственность уже не гарантирована, нужно начальное приближение. Сложность возрастает в разы, сравнительно с линейной моделью, в число раз, равное числу расчёта поправок, уточняющих коэффициенты.
4. Метод Левенберга-Марквардта не требует хотя бы частичной линейности функции, но в данном случае она есть. Задача линейна по двум параметрам из четырёх, а после очевидного преобразования (которое ещё и исключает появление бесконечного числа локальных минимумов, для борьбы с которыми надо вводить ограничения на фазу) - по трём из четырёх. Это позволяет рассматривать задачу, как композицию двух - одномерной оптимизации по нелинейно входящему параметру и линейной по прочим.
5. Перебор по сетке (тут именуемый "сканированием") чрезвычайно просто реализуется. Этим его достоинства исчерпываются. Если он даёт хорошую (тут уместнее говорить "удовлетворительную") точность - требует большой затраты ресурсов, особенно в многомерном случае ("экспоненциальный взрыв"). Однако им приходится пользоваться, если нет иного способа получить начальное приближение, уточняя грубое его приближение более тонким методом.
6. Поскольку модель хорошо определена, можно попытаться поискать лучшее начальное приближение, чем просто тыкать во все подряд точки сетки. Например, можно считать переходы через ноль, или, в данном случае, через среднее значение (если в одном направлении. то расстояние между ними равно периоду синусоиды, если в обоих - полупериоду), а также искать расстояние между максимумами или минимумами. Другой подход - отыскать разности разных порядков (например, первого. убирая свободный член, и третьего. чтобы синус перешёл в синус той же фазы, отличающийся от синуса в первых разностях амплитудой, причём коэффициент зависит от частоты и позволяет её оценить.
7. Затем, используя какой-либо метод одномерной оптимизации ("золотого сечения", квадратичный и т.п.), целевой функцией которого будет мера приближения линейной моделью при заданном значении нелинейно входящего параметра (сумма квадратов отклонений, например), уточняем этот параметр. Если, скажем, для ЗС точность на каждом шаге растёт в 1.618 раз, то для уточнения в сто раз сравнительно с начальным интервалом понадобится около 10 шагов вместо ста по сетке, а в тысячу 14 вместо тысячи.

[svv]

Kevsh
Код на MATLAB.

Используется синтаксис 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))
xx=1:0.1:n;
yy=c_(1)+c_(2)*cos(2*pi*b_*xx)+c_(3)*sin(2*pi*b_*xx);
plot(x,y,'ro', xx,yy,'b');
% plot(b,S);

Я только начал изучать язык MATLAB, возможно, что-то делается проще или эффективнее (хотя у меня на ноуте десятилетней давности всё и так вычисляется за одну секунду).

(Основной график)

Когда построите основной график, раскомментируйте последний оператор, чтобы посмотреть, как зависит от $b$ величина
$\sqrt{\sum\limits_{i=1}^n (c_1u_1(x_i)+c_2u_2(x_i)+c_3u_3(x_i)-y_i)^2}$

Задавайте вопросы!