Метод Симпсона. Вычислить определённый интеграл [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]).