Проект по физике: ряды Фурье, метод Эйлера и параллельное программирование

В третьем семестре мы разработали проект по физики - физический маятник. Я занимался вычислениями. Если в общих чертах, то там решалось дифференциальное уравнение методом Эйлера. И функция заданная эти уравнением раскладывалась в ряд Фурье.
Подсчёт коэффициентов Фурье занимал достаточно много времени (на компах в терминалках на это уходило до двух минут), поэтому решил попробовать ускорить это всё, используя OpenMP.
Сама задача не хитрая - распараллелить цикл подсчёта интегральных сумм(метод трапеций). Однако, в силу неявного задания функции очень неудобно считать её значения в точках, т.к. высчитываются они отдельно от предыдущих. Сказано - сделано.
Теперь просто - запускаем параллельный цикл с помощью #pragma parallel for и использованием reduction и получаем уже расспареллеленую версию интегрирования.
Что меня поразило, так это ускорение: если последовательный вариант считался 26.6 с, то параллельный 3.3 с.
Приведу схематичный пример кода (на С), который был и который стал(пояснения после)

double Ssin = 0, Scos = 0;
double fp = f(x1,fimax,vl,x1,&penenq,vl,w);
double fp1,fc,fc1,fs,fs1;

while(absol(x1-l)>=h){

 fp1 = f(x1, fp, vl, x1+h,&penenq,vl,w);
 fc = fp * cos(2*n*(x1)*pi/l)*h;
 fc1 = fp1 * cos(2*n*(x1+h)*pi/l)*h;
 fs = fp * sin(2*n*(x1)*pi/l)*h;
 fs1 = fp1 * sin(2*n*(x1+h)*pi/l)*h;
 Scos += (fc+fc1)/2;
 Ssin += (fs+fs1)/2;
 x1 += h1;
 fp = fp1;

};
Scos = 2*Scos/l;
Ssin = 2*Ssin/l;

Параллельный вариант:

double Ssin = 0, Scos = 0;
double fp1,fc,fc1,fs,fs1;
int it_numb = l/h;
double *fp = new double [it_numb];
fp[0] = f(x1,fimax,vl,x1,&penenq,vl,w);
for(int i = 1; i < it_numb; i++)
 fp[i] = f(x1+i*h);
int i;

#pragma omp parallel for reduction (+:Scos) reduction (+:Ssin) private (fc,fc1,fs,fs1,i) schedule(dynamic,4)

for(i = 1; i < it_numb; i++){
 fc = fp[i-1] * cos(2*n*(x1+i*h)*pi/l)*h;
 fc1 = fp[i] * cos(2*n*(x1+(i+1)*h)*pi/l)*h;
 fs = fp[i-1] * sin(2*n*(x1+i*h)*pi/l)*h;
 fs1 = fp[i] * sin(2*n*(x1+(i+1)*h)*pi/l)*h;
 Scos += (fc+fc1)/2;
 Ssin1 += (fs+fs1)/2;
}

Scos = 2*Scos/l;
Ssin = 2*Ssin/l; 
delete fp[];

x1 - текущая точка(в последовательном варианте) или начальная точка(в параллельном), l - промежуток интегрирования (если быть более точным [0,l]), h - шаг разбиения.

Я считаю, что сами изменения кода не очень большие, а вот результат очень хороший, так что при возможности - распараллеливание вычисления.