Skip to content

Instantly share code, notes, and snippets.

@shegeley
Created November 20, 2018 14:36
Show Gist options
  • Save shegeley/e707a0755617370a2fc5d9a15c203d69 to your computer and use it in GitHub Desktop.
Save shegeley/e707a0755617370a2fc5d9a15c203d69 to your computer and use it in GitHub Desktop.
Лабораторная работа по предмету "Введение в специальность" №1
% Шаблон (версия от 15.02.2016) предназначен
% для использования студентами каф. ПМиИ СамГТУ
% при оформлении отчетов по лабораторным работам.
% Copyright (c) 2016 by Mikhail Saushkin (msaushkin@gmail.com)
% All rights reserved except the rights granted by the
% Creative Commons Attribution 4.0 International Licence
% <https://creativecommons.org/licenses/by/4.0/>
% Свежая версия шаблона здесь <https://www.overleaf.com/read/sqvxbnhgxxdm>
\documentclass[14pt,a4paper,report]{ncc}
\usepackage[a4paper, left=2cm, right=2cm, top=2cm, bottom=2cm]{geometry}
\usepackage[utf8]{inputenc}
\usepackage[english,russian]{babel}
\usepackage{indentfirst}
\usepackage[dvipsnames]{xcolor}
\usepackage[colorlinks]{hyperref}
\usepackage{listings}
\usepackage{enumitem}
\usepackage{multirow}
\usepackage{multicol}
\usepackage{amssymb}
\usepackage{color}
\usepackage{csvsimple}
\newcommand{\ignore}[1]{}
\lstset{
belowcaptionskip=1\baselineskip,
breaklines=true,
xleftmargin=\parindent,
language=Haskell,
identifierstyle=\color{blue},
stringstyle=\color{orange},
keywordstyle=\color{violet}
}
\begin{document}
% Оформление титульного листа
\begin{titlepage}
\begin{center}
{\scriptsize{\textbf{ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ АВТОНОМНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧЕРЕЖДЕНИЕ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ}}} \\
\textsc{ «БЕЛГОРОДСКИЙ ГОСУДАРСТВЕННЫЙ НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ УНИВЕРСИТЕТ» \\ (НИУ «БелГУ») \\[5mm] Институт инженерных и цифровых технологий \\[2mm] Кафедра общей математики}
\vfill
\textbf{ОТЧЁТ ПО ЛАБОРАТОРНОЙ РАБОТЕ №1 \\[3mm] по дисциплине «Введение в специальность»\\[6mm] на тему: «Определение площади геометрических фигур методом Монте-Карло»}
\end{center}
\hfill
\begin{minipage}{.5\textwidth}
Выполнил: студент \\[2mm] Шепелев Григорий Михайлович\\ группы 12001804\\[5mm]
Проверил:\\[2mm] заведующий кафедрой, доктор технических наук, профессор\\ Аверин Геннадий Викторович
\end{minipage}
\vfill
\begin{center}
Белгород, \theyear\ г.
\end{center}
\end{titlepage}
\newpage
\section*{Задача}
Используя метод Монте-Карло определить площадь геометрической фигуры согласно варианту задания. Найти среднее и дисперсию экспериментальных данных по значениям площади. Сравнить значения полученной площади фигуры с её точным значением, используя для этого математические формулы. \\
\subsection*{Вариант 3}
Геометрические фигуры для исследования:
\begin{enumerate}[label=(\alph*)]
\item Фигура, образованная первой полуволной синусоиды $ y = \sin x $, функцией $ y = \cos x $ и прямой $ x = 0 $
\item Кольцо с центром в начале координат и радиусом 2 и 1.
\item Эллипс с центром в начале координат с полуосями $ a = 1 $, $ b = 2 $ и вырезанный круг с радиусом 1
\item Круг радиуса 1 с вырезанным треугольником с вершинами $ A(-0.5; 0) $, $ B(0.5; 0) $, $ C(0; 1) $
\end{enumerate}
\newpage
\section*{Выполнение}
\subsection*{Введение}
\subsubsection*{Техническая информация}
Для выполнения работы я выбрал функциональный язык программирования Haskell (версия GHC (Glasgow Haskell Compiler) - 1.9.1). Среда разработки - Visual Studio Code, операционная система - Ubuntu 18.04 LTS.
\subsubsection*{Краткое описание метода Монте-Карло}
Метод Монте-Карло используется для моделирования случайных процессов и заключается в том, что для анализа процесса используется генератор случайных величин. На основе полученных данных вычисляют вероятностные характеристики задачи.
\subsection*{Часть a}
\textit{Фигура, образованная первой полуволной синусоиды $ y = \sin x $, функцией $ y = \cos x $ и прямой $ x = 0 $}
Площадь фигуры можно вычилить заранее:
$$ S = (\int_0^{\frac{\pi}{4}} \sin x) + (\int_{\frac{\pi}{4}}^{\frac{\pi}{2}} \cos x) \approx 0.59 $$
Пусть $ x \in [0; \frac{\pi}{2}] $. Аналитическое выражение для пренадлежности точки образованной фигуре можно задать так:
\[
f (x, y) =
\begin{cases}
\text{True}: \text{ if } (y \leq \cos x) \text{ and } (y \leq \sin x) \\
\text{False}: \text{ otherwise }
\end{cases}
\]
На Haskell'е его выглядит так:
\begin{lstlisting}
isInRange :: Double -> Double -> Bool
isInRange x y
| (y <= sin x) && (y <= cos x) = True
| otherwise = False
\end{lstlisting}
Реализация метода Монте-Карло для анализа площади заданной фигуры выглядит так:
\begin{lstlisting}
monteCarlo :: Int -> IO Double
monteCarlo n = do
xs <- sequence $ replicate n $ randomRIO (0 :: Double, pi/2 :: Double)
ys <- sequence $ replicate n $ randomRIO (0 :: Double, 1 :: Double)
return $ (pi/2) * (fromIntegral (strikes xs ys) / fromIntegral n)
where
strikes xs ys = (length . filter (==True)) $ map (\el -> isInRange (fst el) (snd el)) (zip xs ys)
\end{lstlisting}
Где \texttt{n} - число прогонов, и \texttt{sequence \$ replicate n \$ randomRIO (0 :: Double, 1 :: Double)} генерирует спиок из n случайных чисел в промежутке 0, 1.
Вызывая функцию несколько раз с заданным \texttt{n}, получаю:
\begin{multicols}{2}
\begin{lstlisting}[language=bash]
*Main> monteCarlo 100
0.6283185307179586
*Main> monteCarlo 200
0.48694686130641796
*Main> monteCarlo 1000
0.578053048260522
*Main> monteCarlo 5000
0.5906194188748811
*Main> monteCarlo 10000
0.5918760559363171
\end{lstlisting}
\columnbreak
{\small Далее в отчете вывод из терминала будет пропущен}
\end{multicols}
Повторяя процесс несколько раз составляю таблицу: \\
\begin{table}[h!]
\centering
\begin{tabular}{|r|l|l|l|l|l|}
\hline
\multicolumn{1}{|c|}{\multirow{2}{*}{\begin{tabular}[c]{@{}c@{}}Номер прогона\end{tabular}}} & \multicolumn{5}{c|}{Число испытаний (n)} \\ \cline{2-6}
\multicolumn{1}{|c|}{} & 100 & 200 & 1 000 & 5 000 & 10 000 \\ \hline
1 & 0,534 & 0,596 & 0,5482 & 0,5976 & 0,5874 \\ \hline
2 & 0,659 & 0,565 & 0,6173 & 0,577 & 0,585 \\ \hline
3 & 0,502 & 0,589 & 0,625 & 0,584 & 0,586 \\ \hline
4 & 0,565 & 0,534 & 0,6 & 0,577 & 0,586 \\ \hline
5 & 0,518 & 0,565 & 0,611 & 0,586 & 0,589 \\ \hline
\multicolumn{1}{|l|}{Среднее} & 0,5556 & 0,5698 & 0,6003 & 0,58432 & 0,58668 \\ \hline
\multicolumn{1}{|l|}{Дисперсия} & 0,2496 & 0,1562 & 0,1747 & 0,092 & 0,0394 \\ \hline
\end{tabular}
\end{table}
\subsection*{Часть b}
\textit{Кольцо с центром в начале координат и радиусом 2 и 1} \\
Пусть $R_1$ и $R_2$ - радиусы внутренней и внешней окружностей кольца соответственно.
Площадь фигуры (для первой четверти):
$$ S_1 = (\int_0^{R_2} \sqrt{R_2^2 - x^2} dx) - (\int_0^{R_1} \sqrt{R_1^2 - x^2} dx) $$
При заданных $ R_1 = 1; R_2 = 2 $ : $ S_1 = \pi - \frac{\pi}{4} \approx 2.356 $ \\
Соответственно для всей всей фигуры: $ S = 2.356 \cdot 4 = 9.424 $ \\
Пренадлежность точки кольцу можно задать аналитическим выражением:
\[
f(x, y)=
\begin{cases}
\text{ True }: & \text{ if } (R_1^2 - x^2) \leq y^2 \leq (R_2^2 - x^2) \\
\text{ False }: & \text{ otherwise }
\end{cases}
\]
На Haskell'е это выражение выглядит так:
\begin{lstlisting}
isInRing :: Double -> Double -> Bool
isInRing x y
| (y^2 >= r1^2 - x^2 ) && (y^2 <= r2^2 - x^2) = True
| otherwise = False
\end{lstlisting}
Реализация метода Монте-Карло:
\begin{lstlisting}
monteCarlo :: Int -> IO Double
monteCarlo n = do
xs <- sequence $ replicate n $ randomRIO (- r2 :: Double, r2 :: Double)
ys <- sequence $ replicate n $ randomRIO (- r2 :: Double, r2 :: Double)
return $ 4 * r2^2 * (fromIntegral (strikes xs ys) / fromIntegral n)
where
strikes xs ys = (length . filter (==True)) $ map (\el -> isInRing (fst el) (snd el)) (zip xs ys)
\end{lstlisting}
Таблица с результатами: \\
\begin{table}[h!]
\centering
\begin{tabular}{|r|l|l|l|l|l|}
\hline
\multicolumn{1}{|c|}{\multirow{2}{*}{Номер прогона}} & \multicolumn{5}{c|}{Число испытаний (n)} \\ \cline{2-6}
\multicolumn{1}{|c|}{} & 200 & 1 000 & 5 000 & 10 000 & 10 000 \\ \hline
1 & 8,16 & 9,52 & 9,696 & 9,44 & 9,4448 \\ \hline
2 & 6,88 & 10,16 & 9,472 & 9,328 & 9,2784 \\ \hline
3 & 10,08 & 9,92 & 9,488 & 9,472 & 9,52 \\ \hline
4 & 9,28 & 8,96 & 8,896 & 9,448 & 9,4064 \\ \hline
5 & 10,24 & 7,68 & 9,648 & 9,264 & 9,2292 \\ \hline
\multicolumn{1}{|l|}{Среднее} & 8,928 & 9,248 & 9,44 & 9,3904 & 9,37576 \\ \hline
\multicolumn{1}{|l|}{Дисперсия} & 1,1876 & 0,9936 & 0,5652 & 0,2998 & 0,3462 \\ \hline
\end{tabular}
\end{table}
\subsection*{Часть c}
\textit{Эллипс с центром в начале координат с полуосями $ a = 1 $, $ b = 2 $ и вырезанный круг с радиусом 1} \\
Каноническое уравнение эллипса $$ \frac{x^2}{a^2} + \frac{y^2}{b^2} = 1 $$
Площадь фигуры (в первой четветри):
$$ S_1 = (\int_0^1 2 \sqrt{1 - x^2}) - (\int_0^1 \sqrt{1 - x^2}) = \frac{\pi}{4} $$
Тогда площадь всей фигуры: $ S = 4 S_1 = \pi $ \\
Пренадлежность точки эллипсу можно задать выражением:
\[
f(x, y)=
\begin{cases}
\text{ True }: & \text{ if } (1 - x^2) \leq y^2 \leq 4 (1 - x^2) \\
\text{ False }: & \text{ otherwise }
\end{cases}
\]
На Haskell'е выглядит так: \\
\begin{lstlisting}
isInEllipse :: Double -> Double -> Bool
isInEllipse x y
| (y^2 >= (1 - x^2)) && (y^2 <= 4 * (1 - x^2)) = True
| otherwise = False
\end{lstlisting}
Реализация метода Монте-Карло:
\begin{lstlisting}
monteCarlo :: Int -> IO Double
monteCarlo n = do
xs <- sequence $ replicate n $ randomRIO (- a :: Double, a :: Double)
ys <- sequence $ replicate n $ randomRIO (- b :: Double, b :: Double)
return $ 2 * a * 2 * b * (fromIntegral (strikes xs ys) / fromIntegral n)
where
strikes xs ys = (length . filter (==True)) $ map (\el -> isInEllipse (fst el) (snd el)) (zip xs ys)
\end{lstlisting}
Таблица с результатами: \\
\begin{table}[h!]
\centering
\begin{tabular}{|r|l|l|l|l|l|}
\hline
\multicolumn{1}{|c|}{\multirow{2}{*}{Номер прогона}} & \multicolumn{5}{c|}{Число испытаний (n)} \\ \cline{2-6}
\multicolumn{1}{|c|}{} & 200 & 1 000 & 5 000 & 10 000 & 10 000 \\ \hline
1 & 3.63 & 2,92 & 3,032 & 3,2 & 3,0456 \\ \hline
2 & 3,6 & 3,48 & 3,28 & 3,172 & 3,1184 \\ \hline
3 & 2,88 & 3,4 & 1,144 & 3,148 & 3,1024 \\ \hline
4 & 2,72 & 3,42 & 3,224 & 3,204 & 3,1472 \\ \hline
5 & 3,44 & 3,52 & 3,08 & 3,192 & 3,1568 \\ \hline
\multicolumn{1}{|l|}{Среднее} & 3,16 & 3,348 & 2,752 & 3,1832 & 3,11408 \\ \hline
\multicolumn{1}{|l|}{Дисперсия} & 0,6526 & 0,4939 & 0,9511 & 0,1524 & 0,2099 \\ \hline
\end{tabular}
\end{table}
\subsection*{Часть d}
\textit{Круг радиуса 1 с вырезанным треугольником с вершинами $ A(-0.5; 0) $, $ B(0.5; 0) $, $ C(0; 1) $} \\
Площадь фигуры: $ S = \pi - 0.5 \approx 2.64 $ \\
Выражение для вычисления пренадлежности точки образованной фигуре:
\[
f(x, y)=
\begin{cases}
\text{ True }: & \text{ if } (y^2 \leq 1 - x^2) \text{ and } ( (y < 0) \text { or } ( y > - 2 |x| + 1 )) \\
\text{ False }: & \text{ otherwise }
\end{cases}
\]
На Haskell выглядит так:
\begin{lstlisting}
isInFigure :: Double -> Double -> Bool
isInFigure x y
| (y^2 <= 1 - x^2) && ( (y < 0) || (y > -2 * abs x + 1)) = True
| otherwise = False
\end{lstlisting}
Релизация метода Монте-Карло:
\begin{lstlisting}
monteCarlo :: Int -> IO Double
monteCarlo n = do
xs <- sequence $ replicate n $ randomRIO (- 1 :: Double, 1 :: Double)
ys <- sequence $ replicate n $ randomRIO (- 1 :: Double, 1 :: Double)
return $ 4 * (fromIntegral (strikes xs ys) / fromIntegral n)
where
strikes xs ys = (length . filter (==True)) $ map (\el -> isInFigure (fst el) (snd el)) (zip xs ys)
\end{lstlisting}
Таблица с результатами: \\
\begin{table}[h!]
\centering
\begin{tabular}{|r|l|l|l|l|l|}
\hline
\multicolumn{1}{|c|}{\multirow{2}{*}{Номер прогона}} & \multicolumn{5}{c|}{Число испытаний (n)} \\ \cline{2-6}
\multicolumn{1}{|c|}{} & 200 & 1 000 & 5 000 & 10 000 & 10 000 \\ \hline
1 & 2,8 & 2,5 & 2,592 & 2,644 & 2,6196 \\ \hline
2 & 2,2 & 2,7 & 2,64 & 2,666 & 2,6384 \\ \hline
3 & 2,76 & 2,68 & 2,66 & 2,632 & 2,6328 \\ \hline
4 & 2,44 & 2,86 & 2,744 & 2,602 & 2,6163 \\ \hline
5 & 2,6 & 2,66 & 2,624 & 2,612 & 2,6236 \\ \hline
\multicolumn{1}{|l|}{Среднее} & 2,56 & 2,68 & 2,652 & 2,6312 & 2,62614 \\ \hline
\multicolumn{1}{|l|}{Дисперсия} & 0,4966 & 0,3579 & 0,239 & 0,1596 & 0,0961 \\ \hline
\end{tabular}
\end{table}
\section*{Выводы}
\begin{itemize}
\item Каждый прогон модели можно рассматривать как одно наблюдение в проводимом эксперименте на модели.
\item C увеличением продолжительности наблюдения отклонение измеряемой величины от ее точного значения уменьшается, поскольку наблюдаемая система переходит в стационарное состояние.
\item Влияние переходных условий можно уменьшить, если увеличить количество экспериментов.
\item Существует предел, за которым увеличение продолжительности прогона модели уже не дает существенного повышения точности результата, измеряемой дисперсией.
\item Имитационное моделирование не ограничивается разработкой модели и написанием соответствующей программы, а требует подготовки и проведения статистического эксперимента. В связи с этим результаты имитационного моделирования следует рассматривать как экспериментальные данные, требующие специальной обработки и анализа.
\end{itemize}
\pagebreak
\section*{Листинги программ}
\subsection*{Часть a}
\begin{lstlisting}[basicstyle=\small]
import System.Random
import Numeric
isInRange :: Double -> Double -> Bool
isInRange x y
| (y <= sin x) && (y <= cos x) = True
| otherwise = False
monteCarlo :: Int -> IO Double
monteCarlo n = do
xs <- sequence $ replicate n $ randomRIO (0 :: Double, pi/2 :: Double)
ys <- sequence $ replicate n $ randomRIO (0 :: Double, 1 :: Double)
return $ (pi/2) * (fromIntegral (strikes xs ys) / fromIntegral n)
where
strikes xs ys = (length . filter (==True)) $ map (\el -> isInRange (fst el) (snd el)) (zip xs ys)
main :: IO ()
main = do
result <- mapM (\x -> monteCarlo x) [100, 200, 1000, 2000, 10000]
print result
\end{lstlisting}
\subsection*{Часть c}
\begin{lstlisting}[basicstyle=\small]
import System.Random
import Numeric
r1 :: Double
r1 = 1.0
r2 :: Double
r2 = 2.0
isInRing :: Double -> Double -> Bool
isInRing x y
| (y^2 >= r1^2 - x^2 ) && (y^2 <= r2^2 - x^2) = True
| otherwise = False
monteCarlo :: Int -> IO Double
monteCarlo n = do
xs <- sequence $ replicate n $ randomRIO (- r2 :: Double, r2 :: Double)
ys <- sequence $ replicate n $ randomRIO (- r2 :: Double, r2 :: Double)
return $ 4 * r2^2 * (fromIntegral (strikes xs ys) / fromIntegral n)
where
strikes xs ys = (length . filter (==True))
$ map (\el -> isInRing (fst el) (snd el)) (zip xs ys)
main :: IO()
main = do
result <- mapM (\x -> monteCarlo x) [100, 200, 1000, 2000, 10000]
print result
\end{lstlisting}
\subsection*{Часть b}
\begin{lstlisting}[basicstyle=\small]
import System.Random
import Numeric
a :: Double
a = 1.0
b :: Double
b = 2.0
isInEllipse :: Double -> Double -> Bool
isInEllipse x y
| (y^2 >= (1 - x^2)) && (y^2 <= 4 * (1 - x^2)) = True
| otherwise = False
monteCarlo :: Int -> IO Double
monteCarlo n = do
xs <- sequence $ replicate n $ randomRIO (- a :: Double, a :: Double)
ys <- sequence $ replicate n $ randomRIO (- b :: Double, b :: Double)
return $ 2* a * 2 * b * (fromIntegral (strikes xs ys) / fromIntegral n)
where
strikes xs ys = (length . filter (==True)) $ map (\el -> isInEllipse (fst el) (snd el)) (zip xs ys)
main :: IO()
main = do
result <- mapM (\x -> monteCarlo x) [100, 200, 1000, 2000, 10000]
print result
\end{lstlisting}
\subsection*{Часть d}
\begin{lstlisting}[basicstyle=\small]
import System.Random
import Numeric
isInFigure :: Double -> Double -> Bool
isInFigure x y
| (y^2 <= 1 - x^2) && ( (y < 0) || (y > -2 * abs x + 1)) = True
| otherwise = False
monteCarlo :: Int -> IO Double
monteCarlo n = do
xs <- sequence $ replicate n $ randomRIO (- 1 :: Double, 1 :: Double)
ys <- sequence $ replicate n $ randomRIO (- 1 :: Double, 1 :: Double)
return $ 4 * (fromIntegral (strikes xs ys) / fromIntegral n)
where
strikes xs ys = (length . filter (==True)) $ map (\el -> isInFigure (fst el) (snd el)) (zip xs ys)
main :: IO()
main = do
result <- mapM (\x -> monteCarlo x) [100, 200, 1000, 2000, 10000]
print result
\end{lstlisting}
\end{document}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment