Зачем программисту микроконтроллеров численные методы?29.01.2024 17:15

Вот на 12 м году работы программистом микроконтроллеров мне наконец-то пригодились численные методы из ВУЗ (овского) курса по математике.

Причем не сколько для работы сколько для pet-проекта на котором я отлаживаю накопленную кодовую базу.

Я делаю солнечный навигатор. Про это у меня есть отдельный подробный текст. https://habr.com/ru/articles/687640/ Там надо, в частности, по измеренной секундомером кварцевого RTC продолжительности светового дня вычислить географическую широту наблюдения. Вычисления следует производить с учетом рефракции Солнца.

Для этого существует формула (1). Её вывод можно найти в любом советском учебнике астрономии для моряков.

$Day_{Duaration}^{hour}(\phi) = \frac{24}{\pi} \ acos( \frac{-sin \ {51}'-sin \ \phi \ sin \ \delta }{cos \ \phi \ cos \ \delta} ) \qquad \qquad \qquad (1)$

тут

Обозначение	Назначение переменной	Ед. измерения
Day	Продолжительность светового дня	часы
phi	Географическая широта	радианы
delta	Солнечное склонение	радианы
pi	число пи	радианы

Вот график этой функции для солнечного склонения -18.74 градусов. Это 26 декабря каждого года.

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

Вот простой Си-код для вычисления формулы (1) на микроконтроллере.

double calc_day_length_duration_refraction_h(double phi_deg, double delta_deg) {
    double day_length_duration_h = 0.0;
    double phi_rad = DEG_2_RAD(phi_deg);
    double delta_rad = DEG_2_RAD(delta_deg);
    double refrac_rad = DEG_2_RAD(51.0 / 60.0);

    double numerator = -sin(refrac_rad) - sin(phi_rad) * sin(delta_rad);
    double denominator = cos(phi_rad) * cos(delta_rad);
    double cos_arg = numerator / denominator;
    if(cos_arg<=1.0){
    	if(-1.0<=cos_arg){
           day_length_duration_h = (24.0 * acos(cos_arg)) / M_PI;
           LOG_DEBUG(PLANETARIUM, "Phi:%7.2f Deg,Delta:%7.2f Deg,Dur:%5.2f h", phi_deg, delta_deg,day_length_duration_h);
    	}else{
        	day_length_duration_h = 25.0;
            LOG_ERROR(PLANETARIUM, "AcosSmallErr %f", cos_arg);
        }
    }else{
    	day_length_duration_h = 0;
        LOG_ERROR(PLANETARIUM, "AcosBigErr %f", cos_arg);
    }
    return day_length_duration_h;
}

Однако для задачи навигации надо по известной продолжительности светового дня вычислить широту

$\phi(Day_{Duaration}^{hour}) = f(Day_{Duaration}^{hour})$

Фоторезистор измерил мне долготу дня 8 часов 14 минут. Какая же широта соответствует этой продолжительности? Эту задачу можно решить графически. Взять рейсшину, угольник и карандашом провести две линии.

Однако как решить эту задачу на микроконтроллере? Внутри микроконтроллера нет бумаги, рейсшины и угольника. Надо написать код на языке Си.

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

Вот тут-то нам как раз и помогут численные методы. Заметьте, что функция (1) монотонная. При отрицательном солнечном склонении функция убывает. Это хорошая новость. Можно воспользоваться бинарным поиском, чтобы найти нужную нам широту (X) по долготе для (Y).

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

double planetarium_day_dur_to_phi(double day_length_duration_h, double delta_deg){
    double phi_deg = -180.0;
    LOG_INFO(PLANETARIUM, "CalcPhiRefr: DayDur: %f H, SunDec: %f Deg", day_length_duration_h, delta_deg);
    if(0.0 < day_length_duration_h) {
        bool res = is_double_equal_absolute(0.0, delta_deg, 0.001);
        if(false == res) {
            double cur_phi_deg = 0.0;
            double max_phi_deg = 88.0;//1.2*ARCTIC_CIRCLE_DEG;
            double min_phi_deg = -88.0; //1.2*ANTARCTIC_CIRCLE_DEG;
            double cur_day_length_duration_diff_abs_h=0.0;
            double prev_day_length_duration_diff_abs_h=0.0;
            double cur_day_length_duration_diff_h=0.0;
             double cur_day_length_duration_h=0.0;
            uint32_t step = 0;

            for(;;) {
                step++;
                cur_phi_deg = AVERAGE_2(min_phi_deg,max_phi_deg);
                cur_day_length_duration_h = calc_day_length_duration_refraction_h(cur_phi_deg, delta_deg);

                cur_day_length_duration_diff_h = day_length_duration_h -cur_day_length_duration_h;
                LOG_DEBUG(PLANETARIUM, "Step:%u,Phi:[%7.3f..%7.3f..%7.3f] Deg->CurDur:%f h,Diff:%f h,",
                		step,
						min_phi_deg, cur_phi_deg,max_phi_deg,
						cur_day_length_duration_h,
						cur_day_length_duration_diff_h
						);
                cur_day_length_duration_diff_abs_h = fabs(cur_day_length_duration_diff_h);

                if(cur_day_length_duration_diff_abs_h



Функция planetarium_day_dur_to_phi перестает работать как только приближение к решению меньше приемлемой для нас погрешности. В данном случае нас устраивает, если мы найдем решение с погрешностью +/- 1на угловая минута.

Отладка

Для исходных данных возьмем продолжительность светового дня 8.23 часов и солнечное склонение -18.742 deg, алгоритм нашел значение широты 55.907…+/-0.016 градусов. Буквально 15 итераций потребовалось на то, чтобы решить уравнение (1).



Плюсы численных методов

++Можно достигнуть любую наперед заданную точность.

++Меньше вычислений, чем в методе перебора.

Минусы численных методов

--Больше вычислений чем в вычислении формулы.

Итоги
Численные методы позволяют получать решения уравнений, для которых нет аналитического решения с любой точностью. 

Links

Фоторезистор = Навигатор https://habr.com/ru/articles/687640/
    
            © Habrahabr.ru