Форум программистов [Powered by Invision Power Board]

Версия для печати темы
Нажмите сюда для просмотра этой темы в оригинальном формате

Форум программистов > C/C++: Для новичков > Метод последовательного интегрирования

Автор: АлексейX86 29.8.2015, 18:50

Нужно вычислить двойной интеграл методом последовательного интегрирования.
Делаю так : разбиваю тело интегрирования сечениями, затем считаю площадь каждого сечения как интеграл (по методу Симпсона), и суммирую площади всех сечений.
Кол-во шагов разбиений на сечения и в методе Симпсона задаются вводом.

Проблема в том, что при разбиении на одно сечение, метод Симпсона при увелечении в нем кол-ва разбиений все увеличивает и увеличивает площадь данного сечения.
Можете указать на мою ошибку.

Метод Симпсона взят из другого кода, там он считал свой интеграл верно. Из изменений в коде этого метода - добавлено x0 в f.

Код


double  X1 = 1;
double  X2 = 2;

double  Y1 = 2;
double  Y2 = 3;

double f(double x, double y){
    return 1/pow((1+x+y),2);
}

double solveBy_Simpson(double n, double a, double b, double x0){
    double h   = (b-a)/(2*n);
    double summ = f(x0,a)+f(x0,b)+4*f(x0,a+h);

    for(int i = 1; i<=n-1; i++){
        double x = a+2*h*i;
        summ += 2*f(x0,x)+4*f(x0,x+h);
    }
    summ = h*summ/3;

    return summ;
}

double solveBy_StepStrip(double nOuter, double nInner){
    double stepX = (X2-X1)/nOuter;
    double summ = 0;

    for(int i = 0; i <= nOuter-1; i++){
        double x0 = X1+(stepX*((double)i));
        summ += solveBy_Simpson(nInner,X1,X2,x0);
    } 
    return summ;
}

Автор: feodorv 30.8.2015, 00:09

У Вас и на первом, и на втором интеграле интервал интегрирования одинаков?

Что сразу обращает на себя внимание. Здесь есть умножение на ширину сечения:

Цитата(АлексейX86 @ 29.8.2015, 18:50

)

summ = h*summ/3;

А здесь нет:

Цитата(АлексейX86 @ 29.8.2015, 18:50

)

return summ;

Автор: АлексейX86 31.8.2015, 13:18

Цитата(feodorv @ 30.8.2015, 00:09

)

У Вас и на первом, и на втором интеграле интервал интегрирования одинаков?

Да, но должно быть так :

Код


 summ += solveBy_Simpson(nInner,Y1,Y2,x0);

(иксы уже при отладке появились, с игреками соответственно - не работает тоже)

Цитата(feodorv @ 30.8.2015, 00:09

)

А здесь нет:
Цитата(АлексейX86 @ 29.8.2015, 18:50 Найти цитируемый пост)
return summ;

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

О алгоритме :
user posted image

Делаю nOuter сечений тела интегрирования. Затем площадь каждого сечения вычисляю как интеграл от f(x0,y), где x0 - данное сечение, в пределах от y1 до y2. Эта площадь считается по формуле Симпсона, где кол-во шагов есть nInner. Затем суммирование площадей всех сечений.
Провел тест : nOuter = 1, при увеличении nInner площадь сечения увеличивается в бесконечность (на 1-ом шаге площадь сечения БОЛЬШЕ чем значение двойного интеграла).
(Может функция где-то на бесконечность уходит, но вроде бы на заданных промежутках такого нету)

График подынтегральной функции :
user posted image

Полный код

Код


#include <conio.h>
#include <iostream>
#include <cmath>
#include <windows.h>
using namespace std;

int nIn  = 0;
int nOut = 0;

double  X1 = 1;
double  X2 = 2;

double  Y1 = 2;
double  Y2 = 3;

void readN(){
cout << "Print count of cross-sections (nOutter):" << endl << ">";
scanf("%d",&nOut);

cout << "\nPrint count of cross-sections (nInner):" << endl << ">";
scanf("%d",&nIn);
}

double f(double x, double y){
    return (1.0)/(((x+y+1.0)*(x+y+1.0)));
}

double solveBy_Simpson(double n, double a, double b, double x0){
    double h = (b-a)/n;
    double summ = 0;
    for(int i = 0; i<n;i++){
        double pLeft = a+(h*(double)i);
        double pRight = a+(h*(double)(i+1));
        double fun = (pLeft+pRight)/2;
        double c = f(x0,fun)*h;
        summ += c;
    }
    return summ;
}

double solveBy_StepStrip(double nOuter, int nInner){
    double stepX = (X2-X1)/nOuter;
    double summ = 0.0;

    double x0 = X1;
    for(int i = 0; i <= nOuter-1; i++){
        summ += solveBy_Simpson(1,Y1,Y2,x0);
        x0 += stepX;
    } 
    return summ;
}

int main(){
readN();
printf("Value is : %f",solveBy_StepStrip(nOut,nIn));
_getch();
return 0;
}

Автор: feodorv 31.8.2015, 17:17

Цитата(АлексейX86 @ 31.8.2015, 13:18

)

домножать в конце не надо.

Вы интеграл заменяете на сумму

Код


sum( i: 0...nOuter-1) func( X1+i*deltaX ) * deltaX

что есть

Код


deltaX * { sum( i: 0...nOuter-1) func( X1+i*deltaX ) }

Сумму Вы вычисляете. На deltaX не умножаете. И не важно, что в func() Вы снова вычисляете интеграл.

И еще. Первый интеграл тоже надо бы посчитать по методу Симпсона smile

Цитата(АлексейX86 @ 31.8.2015, 13:18

)

Провел тест : nOuter = 1, при увеличении nInner площадь сечения увеличивается в бесконечность (на 1-ом шаге площадь сечения БОЛЬШЕ чем значение двойного интеграла).

Вообще не понял,о чем Вы. Если бы Вы делали nOuter = 1, то код выглядел бы так:

Цитата(АлексейX86 @ 31.8.2015, 13:18

)

    nOuter = 1;
    for(int i = 0; i <= nOuter-1; i++){
        summ += solveBy_Simpson(nInner,Y1,Y2,x0);
        x0 += stepX;
    }

Вы же фиксировали nInner в 1, а nOuter меняете, так? Так вот умножать результат (сумму площадей) на stepX надо, чтобы объём получить...

Автор: АлексейX86 31.8.2015, 17:37

Да, Вы правы, спасибо, исправил.

Цитата(feodorv @ 31.8.2015, 17:17

)

Вообще не понял,о чем Вы. Если бы Вы делали nOuter = 1, то код выглядел бы так:
Цитата(АлексейX86 @  31.8.2015,  13:18 Найти цитируемый пост)

    nOuter = 1;
    for(int i = 0; i <= nOuter-1; i++){
        summ += solveBy_Simpson(nInner,Y1,Y2,x0);
        x0 += stepX;
    }

Вы же фиксировали nInner в 1, а nOuter меняете, так? Так вот умножать результат (сумму площадей) на stepX надо, чтобы объём получить...

С учетом Вашей поправки :

Цитата(feodorv @ 31.8.2015, 17:17

)

Сумму Вы вычисляете. На deltaX не умножаете. И не важно, что в func() Вы снова вычисляете интеграл.

и пока без :

Цитата(feodorv @ 31.8.2015, 17:17

)

И еще. Первый интеграл тоже надо бы посчитать по методу Симпсона smile

программа работает.
Но странно, при малом кол-во сечений тела, выходит результат больший, чем сам объем тела. Но если увеличить число сечений (до 100), то интеграл вычилсляется верно.
И теперь при фиксированном кол-во сечений, интеграл, считаемый по формуле Симпсона сходится к числу, а не уходит в бесконечность.

Автор: feodorv 31.8.2015, 17:38

А что Вы сделали с solveBy_Simpson??? И почему n имеет тип double???

Автор: АлексейX86 31.8.2015, 17:43

И все же....
Не могу понять зачем там нужно было на stepX умножать.
Вот например есть колбаса у меня. Я ее взвешиваю (m0). Затем я ее нарезаю перпендикулярно на 10 кусков (m1-m10). При нарезке часть колбасы останется на ноже, часть испарится, и т.д. Если я сложу веса всех частей (m1-m10), то я получу m0-погрешность. Вообщем, по такой идее я и писал первую версию программы.
Может быть я не учитываю то, при бесконечно малом уменьшении размера ломтя (а соответственно и массы), будет большее кол-во разрезов, и более большая часть массы будет теряться.

Автор: feodorv 31.8.2015, 17:45

Цитата(АлексейX86 @ 31.8.2015, 17:37

)

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

Совсем не странно. У Вас интегрируемая функция падает что по x, что по y, соответственно, прямоугольники, по которым Вы считаете площадь, будут покрывать все пространство под функцией, да ещё немного выступать - вычисленная площадь получиться больше реальной, да ещё по объему та же ситуация.

Добавлено через 2 минуты и 35 секунд

Цитата(АлексейX86 @ 31.8.2015, 17:43

)

Не могу понять зачем там нужно было на stepX умножать.

Интересно, а то что нужно умножать на h во внутреннем интеграле, Вам понятно? Объем слайса есть его площадь умноженная на толщину этого объема.

Автор: АлексейX86 31.8.2015, 17:57

Все, понял.
Я что-то зациклился что толщина сплайса единица, и соответственно домножать не надо.... smile