Вычислить определённый интеграл методом Монте-Карло

@EARLBKV · Регистрация: 19.01.2021

Студворк — интернет-сервис помощи студентам

Нужно решение в Maple для определённого интеграла через метод Монте-Карло

Интеграл:
int(1/(x^(2)+16),x=0..2);

@Royal_X · 24.12.2021, 15:48

EARLBKV, Монте-Карло довольно неточный метод, если сравнивать с другими методами численного интегрирования. Этот метод используется только для сложных интегралов, когда другие более точные не могут найти или неэффективно (тратят много времени) находят интеграл.
Что касается Maple, то в ней можно выбрать этот метод только для кратных интегралов. Даже если вы напишите:

Haskell
1
evalf(Int(1/(x^2 + 16), x = 0 .. 2, _MonteCarlo))

программа выдаст ошибку Error, (in evalf/int) one-dimensional integral, method=_MonteCarlo is invalid.
Кстати, в справке Maple этот метод упоминается в разделе Multiple Integration Methods и предупреждается о ее плохой точности: Monte Carlo method; for low accuracy only (less than 5 digits of accuracy).
Но вот Maple способна вычислить ваш интеграл очень точно другими методами численного интегрирования.
Если нужен именно метод Монте-Карло, то можно это сделать в Wolfram Mathematica, функция NIntegrate которой разрешает выбрать методы MonteCarlo, QuasiMonteCarlo, AdaptiveMonteCarlo, AdaptiveQuasiMonteCarlo даже для 1-D интегралов (интегралы Римана), как вот ваш интеграл.

@murom2013 · 13.01.2022, 16:22

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

1. В системе maple интегрирование методом Монте-Карло реализовано только для кратных интегралов. Ничто не мешает при этом ввести новую переменную и представить область интегрирования как цилиндр с соответствующим основанием:

Haskell
1
2
3
4
5
restart; 
y := x-> x*sin(x); 
int1 := evalf(int(y(x), x = 0 .. Pi)); 
int2 := evalf(Int(y(x), x = 0 .. Pi, t = 0 .. 1, method = _MonteCarlo, epsilon = 0.001)); 
Delta := evalf[5](100*abs(int2-int1)/int1)

Здесь Delta -- относительная погрешность расчета:
Название: mc_1.JPG
Просмотров: 521

Размер: 12.0 Кб

2. Большей наглядности работы метода мы добьёмся при непосредственном программировании:

Haskell
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
restart; 
y := x-> x*sin(x); 
a := 0; b := Pi; 
intrange := a .. b; 
 
N := 10000; 
 
Digits := 20; 
A := Optimization[Maximize](y(x), x = intrange)[1]; 
 
Pv := []; Pn := []; 
st := time[real](); 
randomize(); 
xx := rand(0 .. 1.); 
 
for i to N do 
ax := b*xx(); 
ay := A*xx(); 
flag := evalb(y(ax) < ay); P := [ax, ay]; 
if flag then 
Pv := [op(Pv), P] 
else Pn := [op(Pn), P] 
end if 
end do; 
 
if nops(Pv) = 0 then Pv := [Pn[1]] end if; #если f(x)=Const`
 
plot([Pn, Pv, y(x)], x = intrange, style = [point, point, line], symbol = solidcircle, symbolsize = 2, color = [red, green, black], axes = boxed);
 
nn := nops(Pn); 
cat("количество точек: ", N); 
Digits := 5; 
intA := evalf(int(y(x), x = intrange)); cat("интеграл(точно): ", intA); 
intB := evalf(A*b*nn/N); cat("интеграл(расчет): ", intB); 
dlta := evalf(100*abs(intA-intB)/intA); cat("отн.погрешность(%): ", dlta); 
stme := evalf(time[real]()-st); cat("время расчета(сек): ", stme);

Вычислить определённый интеграл методом Монте-Карло

3. Эту же программу можно оформить в виде несложного маплета:

Кликните здесь для просмотра всего текста

Haskell
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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
##
##kernelopts(version)
##Maple 2017.0, X86 64 WINDOWS, May 17 2017, Build ID 1231047
##
##Маплет работает только в верхней полуплоскости!
##
restart; 
with(Maplets); 
with(Maplets[Elements]); 
with(Maplets[Tools]); 
 
tff := TextField['tf1'](28, "x*sin(x)"); 
tfa := TextField['tf2'](7, "0"); 
tfb := TextField['tf3'](7, "Pi"); 
tfn := TextField['tf4'](5, "1000"); 
chb := CheckBox['cb']("Вывод графика", true); 
 
plotopt := labels = ["", ""], axes = boxed; 
plt := Plotter['pl'](height = 300, plot(x*sin(x), x = 0 .. Pi, plotopt)); 
tbx := TextBox['tb'](editable = false, height = 6, font = Font("Times", 12)); 
 
MCPlot := proc () 
local f, a, b, minf, maxf, Pv, Pn, st, xx, ax, ay, 
flag, i, p, P, m, N, IntA, IntB, dlta, stme; 
 
f := Get('tf1'::algebraic); 
a := Get('tf2'::algebraic); 
b := Get('tf3'::algebraic); 
 
Digits := 20; 
minf := Optimization[Minimize](f, x = a .. b)[1]; 
maxf := Optimization[Maximize](f, x = a .. b)[1]; 
 
Set('pl' = plot(f, x = a .. b, 0 .. maxf, plotopt)); 
Set('tb' = ""); 
 
m := maximize(f, x = a .. b); 
if m = infinity then 
error "Функция разрывна на интервале!" 
end if; 
 
if minf < 0 then 
error "Минимум функции на интервале должен быть неотрицательным!" 
end if; 
 
Set('tb' = "     ... производится расчет ..."); 
N := Get('tf4'::algebraic); 
 
Pv := []; Pn := []; 
st := time[real](); 
randomize(); 
xx := rand(0 .. 1.); 
for i to N do 
ax := (b-a)*xx()+a; 
ay := maxf*xx(); 
flag := evalb(evalf(subs(x = ax, f)) < ay); P := [ax, ay]; 
if flag then Pv := [op(Pv), P] 
else Pn := [op(Pn), P] 
end if 
end do; 
if nops(Pv) = 0 then Pv := [Pn[1]] end if; # если f(x)=Const`
 
if Get('cb') then 
p := plot([Pv, Pn, f], x = a .. b, style = [point, point, line], 
symbol = solidcircle, symbolsize = 1, color = [green, red, black], plotopt); 
Set('pl' = p) 
end if; 
 
Digits := 5:
IntA := evalf(int(f, x = a .. b)); IntB := evalf(maxf*(b-a)*nops(Pn)/N); 
dlta := evalf(100*abs(IntA-IntB)/IntA); stme := evalf(time[real]()-st); 
 
Set('tb' = "     функция y = "); Set('tb'('append') = f); Set('tb'('append') = ",   "); 
Set('tb'('append') = a); Set('tb'('append') = "..."); Set('tb'('appendline') = b); 
Set('tb'('append') = "     количество точек:  "); Set('tb'('appendline') = N); 
Set('tb'('append') = "     интеграл(точно) = "); Set('tb'('appendline') = IntA); 
Set('tb'('append') = "     интеграл(расчет) = "); Set('tb'('appendline') = IntB); 
Set('tb'('append') = "     отн. погрешность (%) = "); Set('tb'('appendline') = dlta); 
Set('tb'('append') = "     время расчета (сек)  = "); Set('tb'('append') = stme) 
end proc; 
 
btnGo := Button("Выполнить", Evaluate(function = MCPlot)); 
btnEx := Button("Закрыть", Shutdown(['tb'])); 
 
hg := HorizontalGlue(); 
Lfunc := [hg, "Функция (от x): ", tff, hg, halign = none]; 
Lpred := [hg, " Пределы интегрирования: от", tfa, "до", tfb, hg, halign = none]; 
Lnast := [hg, "  Количество точек: ", tfn, hg, chb, hg, halign = none]; 
Lbtns := [btnGo, btnEx]; 
Ltbx := [hg, tbx, hg, halign = none]; 
 
L := [Lfunc, Lpred, Lnast, plt, Ltbx, Lbtns]; 
 
s := "Вычисление определенного интеграла методом Монте-Карло"; 
w := Window['w1'](s, L, resizable = false); 
m := Maplet(w); 
Display(m)

@Royal_X · 13.01.2022, 20:52

Сообщение от EARLBKV

Нужно решение в Maple для определённого интеграла через метод Монте-Карло
Интеграл: int(1/(x^(2)+16),x=0..2);

Быстрее будет, если реализовать на языке С++, а не в Maple. Можно написать такой простой код:

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
#include <iostream>
#include <cstdlib>
#include <ctime>
using namespace std;
 
double f(double);
double MonteCarlo(double, double, int);
 
 
int main()
{
    cout << "N = ";
    int n; cin >> n;
    cout.precision(16);
    cout << MonteCarlo(0., 2., n);
}
 
double f(double x)
{
    return 1. / (x * x + 16.);
}
 
double MonteCarlo(double lowerBound, double upperBound, int n)
{
    srand(time(nullptr));
    double sum = 0.;
    for (int i = 0; i < n; ++i)
        sum += f(rand() / static_cast<double>(RAND_MAX) * (upperBound - lowerBound) + lowerBound);
    return (upperBound - lowerBound) * sum / n;
}

Новые блоги и статьи Все статьи Все блоги /
SDL3 для Web (WebAssembly): Идентификация объектов на Box2D v3 - использование userData и событий коллизий 8Observer8 02.03.2026 Содержание блога Финальная демка в браузере. Итоговый код: finish-collision-events-sdl3-c. zip https:/ / www. cyberforum. ru/ blog_attachment. php?attachmentid=11680&d=1772460536 Одним из. . .	Реалии Hrethgir 01.03.2026 Нет, я не закончил до сих пор симулятор. Эта задача сложнее. Не получилось уйти в плавсостав, но оно и к лучшему, возможно. Точнее получалось - но сварщиком в палубную команду, а это значит, в моём. . .	Ритм жизни kumehtar 27.02.2026 Иногда приходится жить в ритме, где дел становится всё больше, а вовлечения в происходящее — всё меньше. Плотный график не даёт вниманию закрепиться ни на одном событии. Утро начинается с быстрых,. . .	SDL3 для Web (WebAssembly): Сборка библиотек: SDL3, Box2D, FreeType, SDL3_ttf, SDL3_mixer и SDL3_image из исходников с помощью CMake и Emscripten 8Observer8 27.02.2026 Недавно вышла версия 3. 4. 2 библиотеки SDL3. На странице официальной релиза доступны исходники, готовые DLL (для x86, x64, arm64), а также библиотеки для разработки под Android, MinGW и Visual Studio. . . .
SDL3 для Web (WebAssembly): Реализация движения на Box2D v3 - трение и коллизии с повёрнутыми стенами 8Observer8 20.02.2026 Содержание блога Box2D позволяет легко создать главного героя, который не проходит сквозь стены и перемещается с заданным трением о препятствия, которые можно располагать под углом, как верхнее. . .	Конвертировать закладки radiotray-ng в m3u-плейлист damix 19.02.2026 Это можно сделать скриптом для PowerShell. Использование . \СonvertRadiotrayToM3U. ps1 <path_to_bookmarks. json> Рядом с файлом bookmarks. json появится файл bookmarks. m3u с результатом. # Check if. . .	Семь CDC на одном интерфейсе: 5 U[S]ARTов, 1 CAN и 1 SSI Eddy_Em 18.02.2026 Постепенно допиливаю свою "многоинтерфейсную плату". Выглядит вот так: https:/ / www. cyberforum. ru/ blog_attachment. php?attachmentid=11617&stc=1&d=1771445347 Основана на STM32F303RBT6. На борту пять. . .	Камера Toupcam IUA500KMA Eddy_Em 12.02.2026 Т. к. у всяких "хикроботов" слишком уж мелкий пиксель, для подсмотра в ESPriF они вообще плохо годятся: уже 14 величину можно рассмотреть еле-еле лишь на экспозициях под 3 секунды (а то и больше),. . .

@EARLBKV 0 / 0 / 0 Регистрация: 19.01.2021 Сообщений: 25

	Maple 2017 Вычислить определённый интеграл методом Монте-Карло 02.12.2021, 11:39. Показов 3638. Ответов 3 Метки нет (Все метки) Нужно решение в Maple для определённого интеграла через метод Монте-Карло Интеграл: int(1/(x^(2)+16),x=0..2); 0

Вычислить определённый интеграл методом Монте-Карло

Решение