Вычисление длинны дуги Эллипса

[vsrmis]

Брякнуло мне в голову вычислить длинну дуги нашей планеты. Лезу в Википедию. Параметры эллипса орбиты есть. Смотрю вычисление длинны дуги сводится к неберущемуся интегралу.

Ну думаю не проблема. Интеграл - площадь между функцией и Ox( простите за столь вульгарное упрощение. Как выразиться более грамотно - я не знаю).
Составил программу на Си

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

alex@alex-Lenovo-B570 /media/Files/Programming/Cpp/C $ cat length.c 
#include "math.h"
#include <stdio.h>


int main(){
        double a = 149597887.5;
        double e = 0.016710219;
        double b =a*sqrt(1-pow(e,2));
        double t;
        double dt = 0.0005;
        double Summ =0;
        for(t = 0;t<360;t+=dt){
                Summ = Summ+dt*sqrt(1-pow(e,2)*pow(cos(t*M_PI/180),2));
                //Summ = Summ+dt*sqrt(a*a*sin(t*M_PI/180)*sin(t*M_PI/180)
                //                  +b*b*cos(t*M_PI/180)*cos(t*M_PI/180));
        }
        Summ = Summ*a;
         printf("%f \n",Summ);
        double L= 4*(M_PI*a*a+ pow(a-b,2))/(a+b);
        printf("%f \n",L);

}
alex@alex-Lenovo-B570 /media/Files/Programming/Cpp/C $ ./a.out 
53851554572.905037 
940016879.701284 
 

Считает полный бред. Скажите, где я ошибся?

И это ещё не всё. Я решил сам выразить длинну эллипса по этой формуле

Для этого я выразил y из канонического уравнения эллипса.
Продифиринцировал его и возвёл в квадрат. Получился простой интеграл. Который в ответе даёт полный бред. Снова, где я ошибся?


Я понимаю, что после взятия интегралла от выраженного Y. Нужно будет умножить на 2, но и это не помогает. Да и Я попенцией чую, что я сделал что-то не так. Потому, что мой получившийся интеграл проще палёной репы, а в ВИКИ написано, что длинна эллипса не выражается в элементарных функциях.

[gris]

Брякнуло в голову, а чуете попенцией. Парадокс, однако.

Вполне достаточное приближение для длины дуги Земли даёт формула: $L=30\cdot (t_2-t_1)$. Время в секундах, длина в километрах. Я сам всегда так считаю и репа не воспаляется.

[vsrmis]

gris в сообщении #603892 писал(а):

Брякнуло в голову, а чуете попенцией. Парадокс, однако.

Вполне достаточное приближение для длины дуги Земли даёт формула: $L=30\cdot (t_2-t_1)$. Время в секундах, длина в километрах. Я сам всегда так считаю и репа не воспаляется.

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

[Mathusic]

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

Цитата:

pow(cos(t*M_PI/180),2)

Вот это что такое? Зачем pow для второй степени?
Вы вообще интеграл каким способом вычисляете?

-- Вт авг 07, 2012 21:11:52 --

Вообще для интеграла должна быть задана либо точность с которой считать, либо кол-во итераций. У вас 360 берется откуда-то.

[Sonic86]

(Оффтоп)

Mathusic в сообщении #603894 писал(а):

Вот это что такое? Зачем pow для второй степени?

Любитель оптимизации детектед

[gris]

А разве в Википедии есть длина орбиты Земли

[Mathusic]

(Sonic86)

Sonic86 в сообщении #603896 писал(а):

(Оффтоп)

Mathusic в сообщении #603894 писал(а):

Вот это что такое? Зачем pow для второй степени?

Любитель оптимизации детектед

Да не. Оптимизация -- это когда думать надо А тут даже этого не надо. Просто прямо руки чешутся взять, да исправить. Или заменить на какую-то самописную square.

[vsrmis]

Mathusic в сообщении #603894 писал(а):

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

Цитата:

pow(cos(t*M_PI/180),2)

Вот это что такое? Зачем pow для второй степени?
Вы вообще интеграл каким способом вычисляете?

-- Вт авг 07, 2012 21:11:52 --

Вообще для интеграла должна быть задана либо точность с которой считать, либо кол-во итераций. У вас 360 берется откуда-то.

Подправил код

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

#include "math.h"
#include <stdio.h>

#define a 149597887.5 //большая полуось
#define e 0.016710219 //Эксцентринитет

#define eps 0.0001;//точность

double GradToRadians(double val) 
//функция пересчёта градусов в радианы
{
        return val*M_PI/180;
}

int main(){

        double b =a*sqrt(1-pow(e,2)); // малая полуось


        double Summ =0; //в конце тут будет результат
        
        for(double t = 0;t<360;t+=eps)//считаю по всей окружности ( если я правильно понял формулу t в градусах 
        {
                //я вычисляю значение функции в каждой точке с погрешностью eps и умножаю её на eps
                //Я так получаю площать маленького прямоугольника
                //И все эти прямоугольники я суммирую
                Summ = Summ+eps*sqrt(1-e*e*cos(GradToRadians(t))*cos(GradToRadians(t)));
        }
        Summ = Summ*a
         printf("%f \n",Summ);
}
 

gris в сообщении #603897 писал(а):

А разве в Википедии есть длина орбиты Земли

К стати нет. Но гуглится она легко.

[Mathusic]

vsrmis в сообщении #603902 писал(а):

Подправил код

Почему объявление переменной double t внутри цикла? (это же C)
Точность не так используют вообще. Для гладкой функции типа эллипса тот же метод прямоугольников хорошую точность даст. А у вас 3.600.000 итераций
Summ = Summ*a
Странно, хотя в семантику вообще не вникал.

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

[Shtorm]

vsrmis в сообщении #603885 писал(а):

....
Для этого я выразил y из канонического уравнения эллипса.
Продифиринцировал его и возвёл в квадрат. Получился простой интеграл. Который в ответе даёт полный бред. Снова, где я ошибся?

.....

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

Ну-ка ещё раз повнимательнее, как Вы разбиваете дробь:

$$\dfrac {x^2}{1-x^2}$$

[vsrmis]

Mathusic в сообщении #603903 писал(а):

vsrmis в сообщении #603902 писал(а):

Подправил код

Почему объявление переменной double t внутри цикла? (это же C)
Точность не так используют вообще. Для гладкой функции типа эллипса тот же метод прямоугольников хорошую точность даст. А у вас 3.600.000 итераций
Summ = Summ*a
Странно, хотя в семантику вообще не вникал.

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

Вы меня погонять по коду решили? Ваша правда - код не особо вылизан. Ну как знаете. Код поправил.

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

#include "math.h"
#include <stdio.h>
#define a 149597887.5 //большая полуось
#define e 0.016710219 //Эксцентринитет

double GtoRs(double val) 
//функция пересчёта градусов в радианы
{
        return val*M_PI/180;
}

double square(double val){
        // возведение в квадрат
        return val*val;
}

double Ellipse(double x){       
        //подинтегральная функция 
        // вот формула http://upload.wikimedia.org/wikipedia/r ... 963f66.png 
        return a*sqrt(1-square(e)*square(cos(GtoRs(x))));
}


double Integral(double right,double left,int pointcount, double(*func)(double) ){
        if(right<left){
                //идём от меньшего к большему
                double temp = right;
                right = left;
                left = temp;
        }
        double dx = (right -left)/pointcount;// вычисляем шаг(ширину нижнего края прямоугольника)
        double Summ=0;// изначально площать нулевая.
        
        while(left<right){
                Summ = Summ + dx*func(left);// умножаю длинну нижнего края прямоугольника на его высоту для получения площади.
                left+=dx;
        }
        
        return Summ;
}

int main(){
double res = Integral(0,360,360,&Ellipse);//интеграл в полярных координатах, в качестве парпметра - угол, от 0 до 360. точность до 1 градуса
printf("%f",res);
return 1 ;
}
 

Shtorm в сообщении #603937 писал(а):

vsrmis в сообщении #603885 писал(а):

....
Для этого я выразил y из канонического уравнения эллипса.
Продифиринцировал его и возвёл в квадрат. Получился простой интеграл. Который в ответе даёт полный бред. Снова, где я ошибся?

.....

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

Ну-ка ещё раз повнимательнее, как Вы разбиваете дробь:

$$\dfrac {x^2}{1-x^2}$$

Да, тут я лопухнулся по чёрному.

[Xaositect]

Так вы же в своей программе считаете не $\int_0^{2\pi} \sqrt{1 - e^2\cos^2 x} dx$, а $\int_0^{360} \sqrt{1 - e^2\cos^2 \frac{\pi x}{180}} dx$.

[vsrmis]

Xaositect в сообщении #603966 писал(а):

Так вы же в своей программе считаете не $\int \sqrt{1 - e^2\cos^2 x} dx$, а $\int \sqrt{1 - e^2\cos^2 \frac{\pi x}{180}} dx$.

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

[Xaositect]

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

[Shtorm]

Xaositect в сообщении #603968 писал(а):

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

Да. Нужно тогда вместо $dx$ писать $d(\frac {\pi x}{180})$, а перед интегралом тогда появится дополнительный множитель $\frac {180}{\pi}$