Форум программистов, компьютерный форум, киберфорум
Наши страницы
Численные методы
Войти
Регистрация
Восстановить пароль
 
Рейтинг 4.80/5: Рейтинг темы: голосов - 5, средняя оценка - 4.80
Quadra
10 / 10 / 0
Регистрация: 29.04.2013
Сообщений: 141
#1

Численное интегрирование кратного интеграла

25.03.2016, 19:55. Просмотров 851. Ответов 10
Метки нет (Все метки)

Доброго времени суток.
Сижу уже пару часов и не понимаю в чём проблема.
Для примера, задача состоит в интегрировании интеграла:
http://www.cyberforum.ru/cgi-bin/latex.cgi?\int_{a}^{b}\int_{c}^{d}(x*y)dxdy

Пробую взять Гауссом.

Квадратурную формулу вывел при помощи Википедии, получилось что-то такое:
http://www.cyberforum.ru/cgi-bin/latex.cgi?\sum_{i=1}^{N}\left[\left(f\left( \frac{({x}_{1}-{x}_{0})(1-\frac{1}{3})}{2} + {x}_{0},\frac{({y}_{1}-{y}_{0})(1-\frac{1}{3})}{2} + {y}_{0}\right) + f\left( \frac{({x}_{1}-{x}_{0})(1-\frac{1}{3})}{2} + {x}_{0},\frac{({y}_{1}-{y}_{0})(1+\frac{1}{3})}{2} + {y}_{0}\right) + f\left( \frac{({x}_{1}-{x}_{0})(1+\frac{1}{3})}{2} + {x}_{0},\frac{({y}_{1}-{y}_{0})(1-\frac{1}{3})}{2} + {y}_{0}\right) + f\left( \frac{({x}_{1}-{x}_{0})(1+\frac{1}{3})}{2} + {x}_{0},\frac{({y}_{1}-{y}_{0})(1+\frac{1}{3})}{2} + {y}_{0}\right)\right) * \frac{({x}_{1}-{x}_{0})({y}_{1}-{y}_{0})}{4}\right]

http://www.cyberforum.ru/cgi-bin/latex.cgi?{x}_{0} = {x}_{i-1},{x}_{1} = {x}_{i}

Ну и реализация:
Кликните здесь для просмотра всего текста
C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
#include <stdio.h>
#include <conio.h>
 
double f(double x, double y)
{
    return x*y;         // [0;1]*[0;1]= 1/4 = 0.25
}
 
double gauss_square(double a, double b, double c, double d, int N)
{
    int i;
    double res, hx, hy, xi_1, xi, yi_1, yi;
 
    hx = (b - a) / N;
    hy = (d - c) / N;
 
    res = 0;
    xi_1 = a;
    xi = xi_1 + hx;
    yi_1 = c;
    yi = yi_1 + hy;
 
    for (i = 0; i < N; i++)
    {
        res += (f((xi - xi_1) * (2.0 / 3.0) / 2.0 + xi_1, (yi - yi_1) * (2.0 / 3.0) / 2.0 + yi_1) +
            f((xi - xi_1) * (2.0 / 3.0) / 2.0 + xi_1, (yi - yi_1) * (4.0 / 3.0) / 2.0 + yi_1) +
            f((xi - xi_1) * (4.0 / 3.0) / 2.0 + xi_1, (yi - yi_1) * (2.0 / 3.0) / 2.0 + yi_1) +
            f((xi - xi_1) * (4.0 / 3.0) / 2.0 + xi_1, (yi - yi_1) * (4.0 / 3.0) / 2.0 + yi_1)) * hx * hy / 4.0;
 
        xi_1 = xi;
        xi = xi_1 + hx;
        yi_1 = yi;
        yi = yi_1 + hy;
    }
 
    res *= N;
 
    return res;
}
 
void main()
{
    double res;
 
    res = gauss_square(0, 1, 0, 1, 100);
    printf("gauss res = %lf", res);
 
    _getch();
}


На единичном квадрате ответ 0.25, но Гаусс выдаёт 0.33

Видимо принципиально что-то не так, ибо интегралы вида http://www.cyberforum.ru/cgi-bin/latex.cgi?\int_{a}^{b}\int_{c}^{d}\alpha{x}^{n}+\beta{y}^{m} dxdy берутся без проблем, даже методом прямоугольников.

Заказываю контрольные, курсовые, дипломные и любые другие студенческие работы здесь.

0
Лучшие ответы (1)
Надоела реклама? Зарегистрируйтесь и она исчезнет полностью.
Similar
Эксперт
41792 / 34177 / 6122
Регистрация: 12.04.2006
Сообщений: 57,940
25.03.2016, 19:55
Ответы с готовыми решениями:

Численное интегрирование
Мои помощники-модераторы и просвещенные в вычислительной математике люди,...

Численное интегрирование
http://pers.narod.ru/study/methods/04.html Внизу этой странице есть таблица с...

Численное интегрирование
Делаю лабу по численному интегрированию. Столкнулся с такой проблемой:...

Численное интегрирование
Квадратурные формулы интерполяционного типа. Нужно написать программку, пытаюсь...

Численное интегрирование
Задан курсовик на тему &quot;Численное интегрирование(Графическая интерполяция)&quot;...

10
mathidiot
Эксперт по математике/физике
2674 / 2378 / 1027
Регистрация: 14.01.2014
Сообщений: 5,056
25.03.2016, 20:50 #2
У Вас стоит одномерная сумма по х, а должна стоять двухмерная по х и у. Фактически берется одномерный интеграл, а не двухмерный!
1
jogano
Модератор
Эксперт по математике/физике
4099 / 2611 / 881
Регистрация: 09.10.2009
Сообщений: 4,587
Записей в блоге: 4
25.03.2016, 20:56 #3
Лучший ответ Сообщение было отмечено Quadra как решение

Решение

Ковыряться в чужом коде - вещь не благодарная...
А у меня всё получается. Для объёма над прямоугольным участком http://www.cyberforum.ru/cgi-bin/latex.cgi?\left[x; \: x+\Delta x \right]\times \left[y; \: y+\Delta y \right] можно, при некотором напряжении мысли, получить такую приближённую формулу:
http://www.cyberforum.ru/cgi-bin/latex.cgi?\Delta V=\left(f\left(x;y \right)+f\left(x+\Delta x,y \right)+f\left(x,y+\Delta y \right)+f\left(x+\Delta x,y+\Delta y \right) \right)\frac{\Delta x \Delta y}{4}
Если отрезок по х разбивается на n одинаковый кусков, а отрезок по y - на m одинаковых кусков, то делается двойное суммирование
http://www.cyberforum.ru/cgi-bin/latex.cgi?\sum_{i=0}^{n-1}\sum_{j=0}^{m-1}\Delta V_{\left[ a+i \Delta x;a+\left(i+1 \right)\Delta x\right]\times \left[ b+j \Delta y;b+\left(j+1 \right)\Delta y\right]}
Вот на единичном квадрате разбил отрезки по х и y на 10 кусков по каждой оси, и выходит ровно 0,25, даже без компьютера.
А у вас почему-то одинарное суммирование.
1
Quadra
10 / 10 / 0
Регистрация: 29.04.2013
Сообщений: 141
25.03.2016, 21:32  [ТС] #4
Ну вот с двумя суммами, только толку...

Кликните здесь для просмотра всего текста
C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
double gauss_square(double a, double b, double c, double d, int N)
{
    int i, j;
    double res, hx, hy, xi_1, xi, yj_1, yj, V;
 
    hx = (b - a) / N;
    hy = (d - c) / N;
 
    res = 0;
    xi_1 = a;
    xi = xi_1 + hx;
    yj_1 = c;
    yj = yj_1 + hy;
 
    for (i = 0; i < N; i++)
    {
        for (j = 0; j < N; j++)
        {
            res += (f((xi - xi_1) * (2.0 / 3.0) / 2.0 + xi_1, (yj - yj_1) * (2.0 / 3.0) / 2.0 + yj_1) +
                f((xi - xi_1) * (2.0 / 3.0) / 2.0 + xi_1, (yj - yj_1) * (4.0 / 3.0) / 2.0 + yj_1) +
                f((xi - xi_1) * (4.0 / 3.0) / 2.0 + xi_1, (yj - yj_1) * (2.0 / 3.0) / 2.0 + yj_1) +
                f((xi - xi_1) * (4.0 / 3.0) / 2.0 + xi_1, (yj - yj_1) * (4.0 / 3.0) / 2.0 + yj_1)) * hx * hy / 4.0;
 
            yj_1 = yj;
            yj = yj_1 + hy;
        }
 
        xi_1 = xi;
        xi = xi_1 + hx;         
    }
 
    res *= N;
 
    return res;
}
0
jogano
Модератор
Эксперт по математике/физике
4099 / 2611 / 881
Регистрация: 09.10.2009
Сообщений: 4,587
Записей в блоге: 4
25.03.2016, 21:44 #5
Обратно код... Ладно.
Перенесите строчку 13 после строчки 16 - вы же при очередном проходе цикла по i должны заново пробегать значения y от с до d, а у вас после первого прохода (для i=0) yj продолжает увеличиваться за пределы отрезка [c;d]
1
Quadra
10 / 10 / 0
Регистрация: 29.04.2013
Сообщений: 141
25.03.2016, 21:44  [ТС] #6
Цитата Сообщение от jogano Посмотреть сообщение
Вот на единичном квадрате разбил отрезки по х и y на 10 кусков по каждой оси, и выходит ровно 0,25, даже без компьютера.
Ваш вариант, 0.75, что не так?

Кликните здесь для просмотра всего текста
C
1
2
3
4
5
6
7
8
9
10
for (i = 0; i < N; i++)
    {
        for (j = 0; j < N; j++)
        {
            res += (f(a + i*hx, b + j*hy) + 
                f(a + (i + 1)*hx, b + j*hy) + 
                f(a + i*hx, b + (j + 1)*hy) + 
                f(a + (i + 1)*hx, b + (j + 1)*hy)) * hx * hy / 4.0;
        }       
    }
0
jogano
Модератор
Эксперт по математике/физике
4099 / 2611 / 881
Регистрация: 09.10.2009
Сообщений: 4,587
Записей в блоге: 4
25.03.2016, 21:51 #7
Цитата Сообщение от Quadra Посмотреть сообщение
Ваш вариант, 0.75, что не так?
Всё так. Только я написал по переменной y значения от b , а у вас же отрезок от с до d. Поменяйте b на c.

Добавлено через 3 минуты
В вашем варианте #4 строчку 12 перенесите тоже ниже строчки 16 - вы ж меняете и значение yj_1.
1
Quadra
10 / 10 / 0
Регистрация: 29.04.2013
Сообщений: 141
25.03.2016, 21:55  [ТС] #8
Цитата Сообщение от jogano Посмотреть сообщение
В вашем варианте #4 строчку 12 перенесите тоже ниже строчки 16 - вы ж меняете и значение yj_1.
Уже поправил.
Кликните здесь для просмотра всего текста
C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
double gauss_square(double a, double b, double c, double d, int N)
{
    int i, j;
    double res, hx, hy, xi_1, xi, yj_1, yj, V;
 
    hx = (b - a) / N;
    hy = (d - c) / N;
    
    res = 0;
    xi_1 = a;
    xi = xi_1 + hx;
 
    for (i = 0; i < N; i++)
    {
        yj_1 = c;
        yj = yj_1 + hy;
 
        for (j = 0; j < N; j++)
        {
            res += (f((xi - xi_1) * (2.0 / 3.0) / 2.0 + xi_1, (yj - yj_1) * (2.0 / 3.0) / 2.0 + yj_1) +
                f((xi - xi_1) * (2.0 / 3.0) / 2.0 + xi_1, (yj - yj_1) * (4.0 / 3.0) / 2.0 + yj_1) +
                f((xi - xi_1) * (4.0 / 3.0) / 2.0 + xi_1, (yj - yj_1) * (2.0 / 3.0) / 2.0 + yj_1) +
                f((xi - xi_1) * (4.0 / 3.0) / 2.0 + xi_1, (yj - yj_1) * (4.0 / 3.0) / 2.0 + yj_1)) * hx * hy / 4.0;
 
            yj_1 = yj;
            yj = yj_1 + hy;
        }
 
        xi_1 = xi;
        xi = xi_1 + hx;
    }
 
    return res;
}

Невнимательность и позднее время делают своё дело. Всем спасибо. Тему можно закрывать.
0
mathidiot
Эксперт по математике/физике
2674 / 2378 / 1027
Регистрация: 14.01.2014
Сообщений: 5,056
25.03.2016, 22:04 #9
Обе квадратурные формулы отлично работают даже при n=1
1
Миниатюры
Численное интегрирование кратного интеграла  
jogano
Модератор
Эксперт по математике/физике
4099 / 2611 / 881
Регистрация: 09.10.2009
Сообщений: 4,587
Записей в блоге: 4
25.03.2016, 22:12 #10
Кстати, в догонку. Зачем каждый раз в цикле вычислять заново разницы xi-xi_1 и yi-yi_1, когда они уже посчитаны в начале программы? Выполняются лишние 8 операций вычитания. Поставив вместо этих разниц hx, hy , вы раза в 2 сокращаете длину формулы в цикле и убираете 6 строчек кода: 11, 16, 25, 26, 29 и 30. Только ж множители i,j поставьте в этой длинной формуле, как в вашем посте #6.

Добавлено через 2 минуты
И умножение на hx*hy/4 можно делать в самом-самом конце, а не в цикле. Машинное время сокращается....
1
Quadra
10 / 10 / 0
Регистрация: 29.04.2013
Сообщений: 141
25.03.2016, 23:13  [ТС] #11
Цитата Сообщение от jogano Посмотреть сообщение
Кстати, в догонку. Зачем каждый раз в цикле вычислять заново разницы xi-xi_1 и yi-yi_1, когда они уже посчитаны в начале программы? Выполняются лишние 8 операций вычитания. Поставив вместо этих разниц hx, hy , вы раза в 2 сокращаете длину формулы в цикле и убираете 6 строчек кода: 11, 16, 25, 26, 29 и 30. Только ж множители i,j поставьте в этой длинной формуле, как в вашем посте #6.
Добавлено через 2 минуты
И умножение на hx*hy/4 можно делать в самом-самом конце, а не в цикле. Машинное время сокращается....
Всё поправлено.
Задача была заставить работать, красивости навелись потом.
0
25.03.2016, 23:13
MoreAnswers
Эксперт
37091 / 29110 / 5898
Регистрация: 17.06.2006
Сообщений: 43,301
25.03.2016, 23:13

Численное интегрирование ОДУ
Решить (написать программу на паскале) задачу Коши методом Адамса-Моултона n=2 ...

Численное интегрирование методом Филона
Ребят, очень нужна помощь. Необходимо написать курсовую на эту тему, &quot;численное...

Численное интегрирование методом Симпсона
Посмотрите пожалуйста правильно ли посчитал интеграл в excel, а то когда делал...


Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
11
Ответ Создать тему
Опции темы

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2018, vBulletin Solutions, Inc.
Рейтинг@Mail.ru