Тема №4. Численное дифференцирование и интегрирование (2 часа)
4.1 Численное дифференцирование функции одной переменной
Пусть функция g(x) задана на отрезке [a, b] в некоторой системе точек:
x |
x0 |
x1 |
... |
xn-1 |
xn |
g(x) |
y0 |
y1 |
... |
yn-1 |
yn |
Известно что у этой функции есть производные всех порядков. Необходимо найти в точке x = c, производную функции g(x).
Первый способ: данный способ заключается в замене функции её многочленом Лагранжа, построенным по таблице значений данной функции, после этого требуемая производная берется от многочлена, пользуясь простотой строения многочлена в смысле правил дифференцирования.
Второй способ: данный способ предполагает, что точка в которой надо найти производную — один из узлов таблицы, например xn. Тогда в качестве искомой первой производной берется число:
(4.1)
По этому принципу можно вычислить у' в точках х1, х2, ..., хn. Затем по первым производным и прежнему принципу можно найти у'' в точках х2, х3, ...,хn, затем у''' в точках х3, х4, ..., хn и так далее.
Пример программы дифференцирования:
Дано: y = sin(x); Найти y' на отрезке [0, π].
Решение:
#include "math.h" // Включаемый файлы с прототипами
#include "stdio.h" // используемых функций и
#include "conio.h" // математических выражений
int main()
{
float x, y, dx, y1, y_1; // Объявление необходимых переменных
clrscr(); // Очистка экрана перед выводом таблицы
dx = 0.1; // Установка шага приращения
y_1 = 0; // Переменная под предыдущее значение у
printf("| x | y | y' |\n"); // Заголовок таблицы
for (x = 0; x <= M_PI; x+=dx) // Цикл обсчета отрезка [0, π]
{
y = sin(x); // Текущее значение у
y1 = (y – y_1)/dx; // Находим y'
printf("| %.2f | %.2f | %.2f |\n", x, y, y1); // Вывод
y_1 = y; // Сохраняем y на следующий цикл
}
return 0; // Конец программы
}
4.2 Численное интегрирование функции одной переменной
Требуется с той или иной степенью точности найти:
(4.2)
Способ №1: Метод прямоугольников.
Отрезок [а, b] на котором ищется интеграл разбивается на n равных частей длиной h.
, xi = a + i·h, i = 1,2, ..., n; (4.3)
Затем на каждом участке [xi-1, xi], функция f(x) заменяется на константу yi-1 = f(xi-1), после чего искомый интеграл заменяется на интеграл от новой ступенчатой функции, т.е. на число S = h·(y0 + y1 + y2 + y3+...+yn-1).
Можно доказать, что справедлива следующая оценка точности:
|S – I| ≤ h·(b – a)·M(1),
где M(1) — максимум модуля первой производной функции f(x) на отрезке [a, b].
//==============================
Интегрирование методом прямоугольников с оценкой точности.
Считается интеграл функции F на отрезке [a,b] с погрешностью
порядка Epsilon.
//==============================
#include "math.h"
double integralrect (const double& a, const double& b, const double& epsilon)
{
double result, h, s1, s2;
int i, n;
n = 1;
h = b-a;
s2 = h*f((a+b)/2);
do
{
n = 2*n;
s1 = s2;
h = h/2;
s2 = 0;
i = 1;
do
{
s2 = s2+f(a+h/2+h*(i-1));
i = i+1;
}
while(i<=n);
s2 = s2*h;
}
while(fabs(s2-s1)>3*epsilon);
result = s2;
return result;
}
Способ № 2: Метод трапеций.
При использовании данного метода искомый интеграл после разбиения отрезка на равные части заменяется на следующую сумму (суммируются площади трапеций, а не прямоугольников, как в предыдущем случае):
(4.4)
где yi = f(xi), i = 0, 1, 2, ..., n.
Можно доказать, что если I — исходный обсуждаемый интеграл, то
(4.5)
где |f''(x)| ≤ M(2) на отрезке [a, b] .
//===============================
Интегрирование методом трапеций с оценкой точности.
Считается интеграл функции F на отрезке [a,b] с погрешностью порядка Epsilon.
//================================
#include "math.h"
double integraltrapezium(const double& a, const double& b, const double& epsilon)
{
int i, n;
double result, h, s, s1, s2;
n = 1;
h = b-a;
s2 = h*(f(a)+f(b))/2;
do
{
s1 = s2;
s2 = 0;
i = 1;
do
{
s2 = s2+f(a-h/2+h*i);
i = i+1;
}
while(i<=n);
s2 = s1/2+s2*h/2;
n = 2*n;
h = h/2;
}
while(fabs(s2-s1)>3*epsilon);
result = s2;
return result;
}
Способ № 3: метод парабол (Симпсона)
В этой ситуации отрезок разбивается на 2n равных частей, где
(4.6)
На участках [x2i-2, x2i], функцию f(x) заменяют на параболу, проходящую через точки (x2i-2, y2i-2), (x2i-1, y2i-1), (x2i, y2i) и интегралом от этой параболы на участке (x2i-2, x2i) заменяют интеграл от функции f(x) на этом же участке, после чего все эти интегралы суммируют и результаты принимают за интеграл от f(x) по всему отрезку [a, b]. Полученная приближенная формула называется формулой парабол или формулой Симпсона:
(4.7)
Если левую часть этого приближенного равенства обозначить через I, а правую через S, то окажется выполненной следующая формула для оценки погрешности:
,
где M(IV) — максимум на интервале (a, b) четвертой производной функции f(x).
#include "math.h"
//================================
Интегрирование методом Симпсона с оценкой точности.
Считается интеграл функции F на отрезке [a,b] с погрешностью
порядка Epsilon.
//================================
double integralsimpson(const double& a, const double& b, const double& epsilon)
{
int i, n;
double result, h, s, s1, s2, s3, x;
s2 = 1;
h = b-a;
s = f(a)+f(b);
do
{
s3 = s2;
h = h/2;
s1 = 0;
x = a+h;
do
{
s1 = s1+2*f(x);
x = x+2*h;
}
while(x<b);
s = s+s1;
s2 = (s+s1)*h/3;
x = fabs(s3-s2)/15;
}
while(x>epsilon);
result = s2;
return result;
}