Лабораторные работы (2008) / Лабраторная работа 1 / Вычисление интеграла функции на отрезке по методу Симпсона (Ясенков)
.docМосковский Энергетический Институт
(Технический Университет)
Кафедра Прикладной Математики
Отчет по лабораторной работе 1
«MPI»
Выполнил: Ясенков Е.М.
Группа: А-13-03
Москва 2008
Задание
Задана некоторая произвольная функция f, непрерывная на отрезке [a, b]. Написать программу параллельного вычисления интеграла функции f на [a, b] по методу Симпсона с точностью eps.
Алгоритм
Значение определенного интеграла
находится методом Симпсона (парабол). Отрезок [a, b] разбивается на n=2m частей x0 =a, x1 =a+h, ..., xn =b с шагом h=(b-a)/n. Вычисляются значения yi = F(xi ) функции в точках xi и находится значение интеграла по формуле Симпсона:
S = Sn +Rn , где
Затем количество точек разбиения удваивается и производится оценка точности вычислений:
Rn = |S2n -Sn |/15
Если Rn > e, то количество точек разбиения удваивается. Значение суммы 2(y1 +y2 +...+y2m-1 ) сохраняется, поэтому для вычисления интеграла при удвоении количества точек разбиения требуется вычислять значения yi лишь в новых точках разбиения.
Реализация
#include "math.h"
#include <stdio.h>
#include "time.h"
#include "mpi.h"
#undef SEEK_SET
#undef SEEK_CUR
#undef SEEK_END
double f(double x) { return (x*x*x*x*x*sin(x))+x; }
int main(int argc, char * argv[])
{
int nProcNum; // число процессов в коммуникаторе
int nProcRank; // ранг данного процесса
double a = 0.0, b = 10.0;// границы отрезка
int n = 0; // число отрезков
double h = 0; // шаг вычисления интеграла
double fIntegr1=0, fIntegr2=0;
double fPartIntegr1, fPartIntegr2, fPartSum1, fPartSum2;
double eps=0.00000000001, x, fa,fb;
int i, nFlag, nIter = 0;
double timebegin, timeend, time;
MPI_Init(&argc, &argv); //инициализация
MPI_Comm_size(MPI_COMM_WORLD, &nProcNum); // число процессов в коммуникаторе
MPI_Comm_rank(MPI_COMM_WORLD, &nProcRank);
if(nProcRank == 0) timebegin = MPI_Wtime();
n = nProcNum;
h = (b-a)/(double)n;// считаем интеграл первый раз
fPartSum1 = 0;
for(i = nProcRank + 1; i <= n; i += nProcNum)
{
x = a + h * ((double)i - 0.5);
fPartSum1 += f(x);
}
fPartSum1 = fPartSum1 * 4;
fPartSum2 = 0;
for(i = nProcRank + 1; i < n; i += nProcNum)
{
x = a + h * i;
fPartSum2 += f(x);
}
fPartSum2 = fPartSum2 * 2;
fPartIntegr1 = fPartSum1 + fPartSum2;
MPI_Reduce(&fPartIntegr1, &fIntegr1, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
fIntegr1 = (f(a) + fIntegr1 +f(b)) * h / 6.0;
MPI_Bcast(&fIntegr1, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
do
{
if(nIter != 0) fIntegr1 = fIntegr2;
n = n * 2;
h = (b-a)/(double)n;
fPartSum1 = 0;
for(i = nProcRank + 1; i <= n; i += nProcNum)
{
x = a + h * ((double)i - 0.5);
fPartSum1 += f(x);
}
fPartSum1 = fPartSum1 * 4;
fPartSum2 = 0;
for(i = nProcRank + 1; i < n; i += nProcNum)
{
x = a + h * i;
fPartSum2 += f(x);
}
fPartSum2 = fPartSum2 * 2;
fPartIntegr2 = fPartSum1 + fPartSum2;
MPI_Reduce(&fPartIntegr2, &fIntegr2, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
if(nProcRank == 0) fIntegr2 = (f(a) + fIntegr2 +f(b)) * h / 6.0;
nIter++;
if(nProcRank == 0) printf("%.15f\n", fIntegr2);
if(nProcRank == 0) nFlag = (fabs(fIntegr1-fIntegr2)/15.0 > eps);
MPI_Bcast(&nFlag, 1, MPI_INT, 0, MPI_COMM_WORLD);
}
while(nFlag!=0);
if(nProcRank == 0)
{
FILE * f = fopen(".\\res.txt", "w");
printf("\nIntegral = %.15f; \nIterations: %d; \nNumber of Intervals: %d \n", fIntegr2, nIter, n);
fprintf(f, "\nIntegral = %.15f; \nIterations: %d; \nNumber of Intervals: %d \n", fIntegr2, nIter, n);
timeend = MPI_Wtime();
time = timeend - timebegin;
printf("time = %f\n", time);
fprintf(f, "time = %f\n", time);
}
MPI_Finalize();
return 0;
}
Тестирование
Зависимость времени вычислений от количества процессов
1 |
2 |
3 |
4 |
6 |
8 |
9 |
12 |
15 |
0,75 |
0,43 |
0,39 |
0,3 |
0,28 |
0,24 |
0,23 |
0,2 |
0,18 |
Зависимость ускорения вычислений от количества процессов
1 |
2 |
3 |
4 |
6 |
8 |
9 |
12 |
15 |
1 |
1,74 |
1,92 |
2,5 |
2,68 |
3,13 |
3,26 |
3,75 |
4,17 |
Вывод
Применение MPI для распараллеливания циклов в данной задаче очень эффективно. С увеличением числа процессов, вычисление интеграла заметно ускоряется.