Метод Симпсона. Вычислить определённый интеграл [latex]I=\int_{a}^{b}{f\left(x \right)}dx[/latex] по формуле Симпсона: [latex]I\approx \frac{b-a}{6n}\left(y_{0}+4y_{1}+2y_{2}+\cdot \cdot \cdot +4y_{2n-1}+y_{2n} \right)[/latex], где [latex]2n[/latex] — количество отрезков разбиения, [latex]y_{0}[/latex], [latex]y_{1}[/latex], …, [latex]y_{2n}[/latex] — значение функции [latex]f\left(x \right)[/latex] на концах отрезков.
В задачах на численное интегрирование определённый интеграл требуется найти с заданной точностью, для чего вычисление по формуле метода рекомендуется проводить многократно, каждый раз уменьшая шаг интегрирования в два раза, пока разница между соседними приближениями не станет меньше заданной погрешности.
| Функция | [latex]a[/latex] | [latex]b[/latex] | [latex]eps[/latex] | Интеграл | Комментарий |
| [latex]f\left(x \right)=\sin \left(x^{2}+2x \right)[/latex] | 1 | 3 | 0.0001 | -0.143058 | Тест пройден. |
| [latex]\ln \left(1+x \right)[/latex] | 1 | 3 | 0.0001 | 2.15888 | Тест пройден. |
| [latex]\tan \left(3x^{3} \right)[/latex] | 2 | 15 | 0.01 | 0.0256033 | Тест пройден. |
| [latex]x\left( x^{2}-1\right)\left(x+1 \right)[/latex] | -1 | 1 | 0.3 | -0.265625 | Тест пройден. |
C++:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 |
#include <iostream> #include <math.h> using namespace std; double f(double x) { return sin(x*x+2*x); } int main() { double a, b, eps;//Нижний и верхний пределы интегрирования (a, b), погрешность (eps). cin >> a >> b >> eps; double I=eps+1, I1=0;//I-предыдущее вычисленное значение интеграла, I1-новое, с большим N. for (int N=2; (N<=4)||(fabs(I1-I)>eps); N*=2) { double h, sum2=0, sum4=0, sum=0; h=(b-a)/(2*N);//Шаг интегрирования. for (int i=1; i<=2*N-1; i+=2) { sum4+=f(a+h*i);//Значения с нечётными индексами, которые нужно умножить на 4. sum2+=f(a+h*(i+1));//Значения с чётными индексами, которые нужно умножить на 2. } sum=f(a)+4*sum4+2*sum2-f(b);//Отнимаем значение f(b) так как ранее прибавили его дважды. I=I1; I1=(h/3)*sum; } cout << I1 << endl; return 0; } |
Java:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 |
import java.util.*; import java.lang.*; import java.io.*; class Ideone { public static double f (double x) { return Math.sin(x*x+2*x); } public static void main (String[] args) { Scanner in = new Scanner(System.in); double a = in.nextDouble();//Нижний и верхний пределы интегрирования (a, b), погрешность (eps). double b = in.nextDouble(); double eps = in.nextDouble(); double I=eps+1, I1=0;//I-предыдущее вычисленное значение интеграла, I1-новое, с большим N. for (int N=2; (N<=4)||(Math.abs(I1-I)>eps); N*=2){ double h, sum2=0, sum4=0, sum=0; h=(b-a)/(2*N);//Шаг интегрирования. for (int i=1; i<=2*N-1; i+=2){ sum4+=f(a+h*i);//Значения с нечётными индексами, которые нужно умножить на 4. sum2+=f(a+h*(i+1));//Значения с чётными индексами, которые нужно умножить на 2. } sum=f(a)+4*sum4+2*sum2-f(b);//Отнимаем значение f(b) так как ранее прибавили его дважды. I=I1; I1=(h/3)*sum; } System.out.println(I1); } } |
С помощью программы можно вычислить интеграл любой непрерывной на [latex]\left[a;b \right][/latex] функции, для этого нужно изменить 7 строку.
В условии главного цикла [latex]N\leq 4[/latex], так как важно, чтобы цикл выполнился хотя бы дважды. Потому что первый раз мы сравниваем новое значение интеграла не с предыдущим вычисленным, а с нулём, и если новое значение интеграла будет меньше погрешности, то цикл прекратится после первого же выполнения (без условия [latex]N\leq 4[/latex]).
Долго смотрел на Ваше решение — чувствуя какой-то скрытый подвох. И вот он — если Вы вычислите приближенное значение интеграла при N=2, и окажется что оно равно нулю или очень близко к нулю, то выполнится условие fabs(I1-I)>eps и будет выдан ответ (ведь в конце первой итерации в переменной I будет «фантомный» ноль, не полученный по формуле Симпсона). Этот ответ будет выдан даже, если точное значение интеграла существенно отличается от нуля. Нужно обеспечить, чтобы выполнилось как минимум две итерации (при N=2 и N=4). Понимаю, что такая ошибка будет очень редко проявлять себя. Но тем более она коварна — в реальной программе такая скрытая ошибка может проявить себя в самый неожиданный момент. Если это просто обычная программа расчета чего-либо, то будет выдан неверный ответ, а если эта программа будет управлять атомной станцией или космическим кораблем может произойти катастрофа. Так что постарайтесь исправить этот мелкий недочет. Достаточно дополнить условие продолжения цикла.
Исправила. Теперь цикл точно выполнится два раза.
Засчитано, отлично, 10 баллов, только три небольшие просьбы:
1) Вместо условия N==4 (правильного в Вашем случае) надежнее и логичнее поставить N<=4.
2) Опишите в отчете этот нетривиальный момент с этим условием.
3) и, наконец, раз это условие уже поставлено, то давайте я придумаю тест, где оно будет существенно: давайте возьмем функцию x*(x*x-1)*(x+1) на интервале от -1 до 1 . Проверьте свою старую версию (без дополнительного условия) и новую версию — в отчет впишите результат правильной программы.
Исправила код и дополнила отчёт.
Отлично!