Нужна помощь в исправлении алгоритма, который приближает пи - PullRequest
Новая
       42

Нужна помощь в исправлении алгоритма, который приближает пи

3 голосов
/ 10 апреля 2020

Я пытаюсь написать код C для алгоритма, который приближается к pi . Предполагается получить объем куба и объем сферы внутри этого куба (радиус сферы равен половине стороны куба). Затем я должен разделить объем куба на сферу и умножить на 6, чтобы получить число пи.

Это работает, но он делает что-то странное в той части, которая должна получить объемы. Я полагаю, что это связано с дельтой, которую я выбрал для приближений. С кубом стороны 4 вместо того, чтобы дать мне объем 64, это дает мне 6400. Сферой вместо 33 это дает мне 3334. что-то.

Может кто-нибудь понять это? Вот код (я прокомментировал соответствующие части): 

#include <stdio.h>      

int in_esfera(double x, double y, double z, double r_esfera){
    double dist = (x-r_esfera)*(x-r_esfera) + (y-r_esfera)*(y-r_esfera) + (z-r_esfera)*(z-r_esfera);

    return  dist <= (r_esfera)*(r_esfera) ? 1 : 0;   
}   

double get_pi(double l_cubo){   
    double r_esfera = l_cubo/2;   
    double total = 0;
    double esfera = 0;    
//this is delta, for the precision. If I set it to 1E anything less than -1 the program continues endlessly. Is this normal?
    double delta = (1E-1);   

    for(double x = 0; x < l_cubo; x+=delta){
        printf("x => %f; delta => %.6f\n",x,delta);
        for(double y = 0; y <l_cubo; y+=delta){
            printf("y => %f; delta => %.6f\n",y,delta);
            for(double z = 0; z < l_cubo; z+=delta){
                printf("z => %f; delta => %.6f\n",z,delta);
                total+=delta;
                if(in_esfera(x,y,z,r_esfera))
                    esfera+=delta;
            }
        }
    }

    //attempt at fixing this
        //esfera/=delta;
        //total/=delta;
    //

//This printf displays the volumes. Notice how the place of the point is off. If delta isn't a power of 10 the values are completely wrong.   
    printf("v_sphere = %.8f; v_cube = %.8f\n",esfera,total);   

    return (esfera)/(total)*6;
}   

void teste_pi(){        
    double l_cubo = 4;    
    double pi = get_pi(l_cubo);

    printf("%.8f\n",pi);
}   

int main(){   
    teste_pi();
}

Ответы [ 2 ]

2 голосов
/ 10 апреля 2020 
total+=delta;
if(in_esfera(x,y,z,r_esfera))
    esfera+=delta;

total и esfera - трехмерные объемы, тогда как delta - одномерная длина. Если бы вы отслеживали единицы измерения, у вас было бы м 3 слева и м справа. Единицы несовместимы.

Чтобы исправить это, куб delta, так что вы концептуально накапливаете крошечные кубики вместо крошечных линий. 

total+=delta*delta*delta;
if(in_esfera(x,y,z,r_esfera))
    esfera+=delta*delta*delta;

Выполнение, которое фиксирует вывод, а также работает для любого значения delta: 

v_sphere = 33.37400000; v_cube = 64.00000000
3.12881250

Обратите внимание, что этот алгоритм "работает" для произвольных delta значений, но имеет серьезные проблемы с точностью. Это невероятно склонно к проблемам округления. Он работает лучше всего, когда delta - это степень двух: 1/64.0 лучше, чем 1/100.0, например: 

v_sphere = 33.50365448; v_cube = 64.00000000
3.14096761

Кроме того, если вы хотите, чтобы ваша программа работала быстрее, избавьтесь от всех этих распечаток! Или, по крайней мере, те, что во внутренних циклах ...

0 голосов
/ 10 апреля 2020

Дело в том, что умножение на целые числа , как a * b * c, равнозначно добавлению 1 + 1 + 1 + 1 + ... + 1 a * b * c раз, верно?

Вы добавляете delta + delta + ... (x / delta) * (y / delta) * (z / delta) раз. Или, другими словами, (x * y * z) / (delta ** 3) раз.

Теперь эта сумма delta s такая же, как эта: 

delta * (1 + 1 + 1 + 1 + ...)
         ^^^^^^^^^^^^^^^^^^^^ (x * y * z) / (delta**3) times

Итак, если delta - это степень из 10, (x * y * z) / (delta**3) будет целым числом, и оно будет равно сумме 1 в скобках (потому что это то же самое, что и product x * y * (z / (delta**3)), где последний член является целым числом - см. самое первое предложение этого ответа). Таким образом, ваш результат будет следующим: 

delta * ( (x * y * z) / (delta ** 3) ) == (x * y * z) / (delta**2)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ the sum of ones

Вот так вы и рассчитали произведение , деленное на delta квадрат .

Чтобы решить эту задачу, умножьте все тома на delta * delta.

Тем не менее, я не думаю, что возможно использовать эту логику c для delta с, которые не являются степенью 10. И действительно, код будет go всех видов haywire для delta == 0.21 и l_cubo == 2, например: вы получите 9.261000000000061 вместо 8.

...