Численные методы +Qt, Под Qt Опции

[QMainWindow]

Добрый день!
Подскажите, пожалуйста, какие-нибудь библиотеки (или классы) для работы с численными методами (в данном случае требуется решить дифф. уравнение). Может быть и на Qt чего уже сделали. Можно, конечно, и самому алгоритм сделать, но хотелось бы уже готовый и оптимизированный.
Заранее спасибо!))

[DEADHUNT]

вот как то проводил эксперимент, какой метод решает задачу Коши первого порядка с наименьшей погрешностью:

Раскрывающийся текст

#include <iostream>
#include <iomanip>
#include <cmath>
#include <algorithm>

/* 
  y' = y * (cos(x) + sin(x))
  y(0) = 1
  y(x) = e(sin(x) - cos(x)) * e
*/

double exact_solution(double x)
{
    return exp(sin(x)) / exp(cos(x)) * exp(1);
}

double f(double x, double y)
{
    return y * (cos(x) + sin(x));
}

// y(x1) = y(x0) + h*f(x0)
double euler(int n)
{
    double h = 1.0 / n;
    double x = 0.0;
    double y = 1;
    double err = 0;
    for (int i = 1; i <= n; ++i, x += h)
    {
        y = y + h * f(x, y);
        err = std::max(err, fabs(y - exact_solution(x)));        
    }
    return err;
}

double rectangle(int n)
{
    double h = 1.0 / n;
    double x = 0.0;
    double y = 1;
    double err = 0;
    for (int i = 1; i <= n; ++i, x += h)
    {
        y = y + h * f(x + h / 2, y);
        err = std::max(err, fabs(y - exact_solution(x)));        
    }
    return err;
}

double simpson(int n)
{
    double h = 1.0 / n;
    double x = 0.0;
    double y = 1;
    double err = 0;
    for (int i = 1; i <= n; ++i, x += h)
    {
        y = y + (h / 6) * (f(x, y) + 4 * f(x + h / 2, y) + f(x + h, y));
        err = std::max(err, fabs(y - exact_solution(x)));        
    }
    return err;
}

double runge(int n)
{
    double h = 1.0 / n;
    double x = 0.0;
    double y = 1;
    double err = 0;
    for (int i = 1; i <= n; ++i, x += h)
    {
        double k1 = f(x, y),
            k2 = f(x + 0.5 * h, y + h / 2 * k1),
            k3 = f(x + 0.5 * h, y + h / 2 * k2),
            k4 = f(x + h, y + h * k3);
        y = y + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4);
        err = std::max(err, fabs(y - exact_solution(x)));        
    }
    return err;
}

int main()
{
    const int n = 80;
    std::cout << std::fixed << std::setprecision(30);
    std::cout << "euler\t" << euler(n) << std::endl;
    std::cout << "rect\t" << rectangle(n) << std::endl;
    std::cout << "simpson\t" << simpson(n) << std::endl;
    std::cout << "runge\t" << simpson(n) << std::endl;
    return 0;
}

в исходнике реализованы следующие методы:
метод эйлера(фактически интегрирование уравнения методом левых прямоугольников) - погрешность O(h)
интегрирование уравнения методов центральных прямоугольников - погрешность O(h^2)
интегрирование уравнения методов симпсона - погрешность O(h^4)
метод рунге-кутта 4 порядка - погрешность O(h^4)
странно но почему то на моих примерах точнее всего считает метод эйлера.