Метод Симпсона. Большая погрешность при большом шаге.

[dimanet]

Здравствуйте!
Стоит задача посчитать методом Симпсона интеграл $\int\limits_{0}^{1}\frac{\ln(1+x)}{x}dx$.
В нуле я доопределил значение подынтегральной функции до 1 (пролопиталил).
Формула, по которой делал:

$$\frac{{b - a}}{{6N}}\left( {{f_0} + {f_{2n}} + 2\left( {{f_2} + {f_4} + \cdots + {f_{2N - 2}}} \right) + 4\left( {{f_1} + {f_3} + \cdots + {f_{2N - 1}}} \right)} \right)$ $

Вот код на си++:

double simp(double lim_a, double lim_b, double h){
    int N = static_cast<int>((lim_b - lim_a)/h);
    double delta = (lim_b - lim_a)/(6*N);
    double result = 0;

    result = (func(lim_a+0.000001) + func(lim_a + h*N) + 1)*delta;

    int i;
    for(i = 2; i <= 2*N-2; i+=2){
        result += func(lim_a + 0.5*h*i)*delta*2;
    }

    for(i = 1; i <= 2*N-1; i+=2){
        result += 4*func(lim_a + 0.5*h*i)*delta;
    }

    return result;
}

(Понимаю, что можно назвать говнокодом.)

На малых значениях шага $h$ погрешность у Симпсона самая маленькая по сравнению с другими методами, как и должно быть. Но преподаватель вводил шаги 0,1 и 0,01 (количество узлов соответственно 10 и 100) и погрешность у метода Симпсона оказывалась самая большая. Сказал, чтоб я разобрался.
Кроме того, я переделал полностью эту функцию по другой формуле и погрешность получается еще большей на больших и даже малых шагах.
Ткните меня пожалуйста носом, где и что я пропустил в программе матанализа и численных методов.
Спасибо.

[AKM]

После этой трагедии решение таких интегралов на форуме запрещено.

Запишите () подынтегральное выражение корректно.

[dimanet]

Да, я сразу понял о чем речь. Спасибо. У нас один профессор из-за этого на одной научной работе написал "Ересь".

[мат-ламер]

dimanet в сообщении #636667 писал(а):

(Понимаю, что можно назвать говнокодом.)

А почему? Если прорамма работает (надеюсь Вы её тестировали), то всё хорошо, и зря Вы тут текст программы привели.

dimanet в сообщении #636667 писал(а):

Ткните меня пожалуйста носом, где и что я пропустил в программе матанализа и численных методов.

У Вас метод Симпсона был, наверное, с оценками погрешности. Почитайте этот вопрос.

[_Ivana]

мат-ламер в сообщении #636679 писал(а):

А почему? Если прорамма работает (надеюсь Вы её тестировали), то всё хорошо, и зря Вы тут текст программы привели.

Категорически не согласен. Ничего не хорошо и код приведен не зря. Можно указать много неоптимальных мест, ошибка может крыться в одном из них.

Вот этот ужас

result = (func(lim_a+0.000001) + func(lim_a + h*N) + 1)*delta;

перепишите по нормальному - берите прямо func(lim_a) (лопитальте прямо внутри функции по условию if (lim_a == 0) else), func(lim_a + h*N) заменяйте на func(lim_b), на delta можно умножить весь результат потом (но тут зависит от величины шага, чтобы не переполнился тип), и главное - откуда там в скобках +1?

[dimanet]

Да, с оценками. Сейчас почитаю отдельно про оценки, но сразу скажу, что ни я, ни Mathematica максимум второй и четвертой производной не смогли найти. Следовательно, таким способом оценить у меня не получилось (это отдельный вопрос). Преподаватель сказал взять значение, которое дает Mathematica --

$$\frac{{{\pi ^2}}}{{12}}$$

, и вычитать от него значение посчитанного программой интеграла.
Спасибо.

-- Сб окт 27, 2012 22:57:42 --

Уже_Ivana, уже переписал

/* Моя функция */
double func(double x){
   if(x == 0) return 1;
   else return log(1+x)/x;
}

Метод Симпсона:

Единица взялась как значение функции в нуле.

Вот новая функция

double simp(double lim_a, double lim_b, double h){
    int N = static_cast<int>((lim_b - lim_a)/h);
    double delta = (lim_b - lim_a)/(6*N);
    double result = 0;

    result = (func(lim_a + h*N) + func(0))*delta;

    int i;
    for(i = 2; i <= 2*N-2; i+=2){
        result += func(lim_a + 0.5*h*i)*delta*2;
    }

    for(i = 1; i <= 2*N-1; i+=2){
        result += 4*func(lim_a + 0.5*h*i)*delta;
    }

    return result;
}

Там как раз наверное этот эпсилон как избавление от особенностей давал большую погрешность. Вроде сейчас погрешность гораздо меньше. Спасибо.

-- Сб окт 27, 2012 23:03:44 --

мат-ламер, _Ivana,
большое спасибо! Теперь все точно работает. Только что проверил.

[_Ivana]

dimanet в сообщении #636688 писал(а):

Единица взялась как значение функции в нуле.

Тогда первое слагаемое func(lim_a+0.000001) было лишнее. И вообще плохой тон задавать значение когда передаются пределы явно - препод скажет взять от 0.5 до 1 и вы что будете делать - все переписывать?
func(0) меняйте на func(lim_a).

Раз уж привели код, то делайте и это

_Ivana в сообщении #636687 писал(а):

func(lim_a + h*N) заменяйте на func(lim_b), на delta можно умножить весь результат потом (но тут зависит от величины шага, чтобы не переполнился тип - это под вопросом)

Вдобавок посчитайте 0.5*h единожды а не в каждом цикле. Это если не считать того, что цикл можно было сделать и один.

-- 27.10.2012, 23:08 --

dimanet в сообщении #636688 писал(а):

Там как раз наверное этот эпсилон как избавление от особенностей давал большую погрешность.

Погрешность вам давал не этот эпсилон а лишнее слагаемое - или единица, или func(lim_a) - надо было оставить что-то одно.

[dimanet]

_Ivana, я его и ранье убирал, но почему-то вообще ересь какую-то программа выдавала, поэтому решил не трогать. Сейчас вот вышло.
Про замену func(0) на func(lim_a) и func(lim_a + h*N) func(lim_b) понял. Ага, и про 0,5h понял. Про единый цикл -- пусть пока так, я попозже разберусь более глубоко в синтаксисе. Я не очень хорошо программирую, мягко говоря.
Спасибо.

-- Сб окт 27, 2012 23:14:23 --

_Ivana в сообщении #636693 писал(а):

Погрешность вам давал не этот эпсилон а лишнее слагаемое - или единица, или func(lim_a) - надо было оставить что-то одно.

Ага, понял. Это ведь одно и то же.

[ewert]

_Ivana в сообщении #636687 писал(а):

на delta можно умножить весь результат потом (но тут зависит от величины шага, чтобы не переполнился тип),

Тип переполниться практически не может -- грубо говоря, для этого нужно брать шаг порядка машинного нуля, а этого никто в здравом рассудке делать не будет. Развлекаться оптимизацией кода можно, конечно, сколько угодно, но какого-то принципиального выигрыша это не даст. А вот что принципиально надо сознавать, так это неудачность самого интерфейса (хотя при предлагавшихся входных данных она и не проявляется). В качестве входного аргумента лучше задавать не шаг, а количество отрезков. Если же хочется задавать именно шаг, то внутри процедуры следует предусмотреть его коррекцию: после определения $N$ надо пересчитать шаг обратным делением длины всего промежутка на $N$.

[_Ivana]

ewert в сообщении #636762 писал(а):

Тип переполниться практически не может

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

ewert в сообщении #636762 писал(а):

Развлекаться оптимизацией кода можно, конечно, сколько угодно, но какого-то принципиального выигрыша это не даст.

Согласен с первым но не согласен со вторым. Смотря что считать принципиальным выигрышем - если существенное уменьшение количества операций и времени выполнения кода не считать выигрышем.... Ответ то не изменится. А лишние 1000 раз в циклах умножения на delta да на 2 и 4 делать - имхо сильно замедлит код.
Про задание шага согласен.

[ewert]

_Ivana в сообщении #636814 писал(а):

я имел в виду что может ухудшиться его точность из-за плавающей запятой,

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

_Ivana в сообщении #636814 писал(а):

лишние 1000 раз в циклах умножения на delta да на 2 и 4 делать - имхо сильно замедлит код.

Замедлит -- но не сильно. Там ведь только для вычисления аргумента функции нужны минимум две операции (можно обойтись и одной, но это немного ухудшит точность). И ещё минимум пять операций для вычисления функции (включая собственно её вызов). И ещё инкремент, и ещё какие-то затраты на собственно организацию цикла.

Вообще исходный код ТС как раз и был оптимизирован -- по краткости кода в сочетании с наглядностью. Любые попытки оптимизировать либо по скорости, либо по точности существенной пользы не принесут, но читабельность ухудшат.