Вычисление интеграла функции на отрезке по методу Симпсона (Ясенков) (Лабораторная работа 1)

Московский Энергетический Институт

(Технический Университет)

Кафедра Прикладной Математики

Отчет по лабораторной работе 1

«MPI »

Выполнил: Ясенков Е.М.

Группа: А-13-03

Москва 2008

Задание

Задана некоторая произвольная функция f, непрерывная на отрезке [a, b]. Написать программу параллельного вычисления интеграла функции f на [a, b] по методу Симпсона с точностью eps.

Алгоритм

З начение определенного интеграла

н аходится методом Симпсона (парабол). Отрезок [a, b] разбивается на n=2m частей x₀ =a, x₁ =a+h, ..., x_n =b с шагом h=(b-a)/n. Вычисляются значения y_i  = F(x_i ) функции в точках x_i  и находится значение интеграла по формуле Симпсона:

S = S_n +R_n , где

Затем количество точек разбиения удваивается и производится оценка точности вычислений:

R_n  = |S_2n -S_n |/15

Если R_n  > e, то количество точек разбиения удваивается. Значение суммы 2(y₁ +y₂ +...+y_2m-1 ) сохраняется, поэтому для вычисления интеграла при удвоении количества точек разбиения требуется вычислять значения y_i  лишь в новых точках разбиения.

Реализация

#include "math.h"

#include

#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;

}

Тестирование

Зависимость времени вычислений от количества процессов

Зависимость ускорения вычислений от количества процессов

Вывод

Применение MPI для распараллеливания циклов в данной задаче очень эффективно. С увеличением числа процессов, вычисление интеграла заметно ускоряется.

4