Главная страница
Навигация по странице:

  • 4.2. Приближение функции многочленами Тейлора

  • 4.3. Интерполяция функции многочленами Лагранжа

  • Погрешность интерполяции.

  • 5. Численное интегрирование функций одной переменной 5.1. Постановка задачи численного интегрирования

  • 5.2. Метод прямоугольников

  • Оценка погрешности.

  • 5.4. Метод Симпсона (метод парабол)

  • 5.5. Правило Рунге практической оценки погрешности

  • Лабы Антонова. Тема Решение задач вычислительными методами. Основные понятия


    Скачать 1.35 Mb.
    НазваниеТема Решение задач вычислительными методами. Основные понятия
    Дата24.10.2018
    Размер1.35 Mb.
    Формат файлаpdf
    Имя файлаЛабы Антонова.pdf
    ТипРешение
    #54374
    страница3 из 5
    1   2   3   4   5
    Тема 4. Приближение функций
    4.1. Постановка задачи
    Задача приближения (аппроксимации) функций заключается в том, чтобы для данной функции построить другую, отличную от нее функцию, значения которой достаточно близки к зна- чениям данной функции. Такая задача возникает на практике достаточно часто. Укажем наиболее типичные случаи.
    1. Функция задана таблицей в конечном множестве точек, а вычисления нужно произвести в других точках.
    2. Функция задана аналитически, но ее вычисление по формуле затруднительно.
    При решении задачи поиска приближенной функции возникают следующие проблемы.
    1. Необходимо выбрать вид приближенной функции. Для приближения широко использу- ются многочлены, тригонометрические функции, показательные функции и т. д.
    2. Необходимо выбрать критерий близости исходной и приближенной функции. Это может быть требование совпадения обеих функций в узловых точках (задача интерполяции), минимизация среднеквадратического уклонения (метод наименьших квадратов) и др.
    3. Необходимо указать правило (алгоритм), позволяющее с заданной точностью найти при- ближение функции.
    4.2. Приближение функции многочленами Тейлора

    Пусть функция y = f(x)определена в окрестности точки a и имеет в этой окрестности n + 1 производную. Тогда в этой окрестности справедлива формула Тейлора:
    f(x)= c
    0
    + c
    1
    (x – a) + c
    2
    (x – a)
    2
    + … + c
    n
    (x – a )
    n
    + R
    n
    (x) = T
    n
    (x) + R
    n
    (x),
    где
    c
    k
    =
    !
    )
    (
    )
    (
    k
    a
    f
    k
    T
    n
    (x) – многочлен Тейлора:
    T
    n
    (x)= c
    0
    + c
    1
    (x – a) + c
    2
    (x – a)
    2
    + … + c
    n
    (x – a )
    n
    , (4.1)
    R
    n
    (x)остаточный член формулы Тейлора. Его можно записать различными способами, например, в форме Лагранжа:
    R
    n
    (x)=
    1
    )
    1
    (
    )
    (
    )!
    1
    (
    )
    (




    n
    n
    a
    x
    n
    f

    , a



    x.
    Многочлен Тейлора (4.1) обладает свойством, что в точке x = a все его производные до по- рядка n включительно совпадают с соответствующими производными функции f, т. е.
    T
    )
    (k
    n
    (a)= f
    (k)
    (a), k = 0, 1, …, n.
    В этом легко убедиться, дифференцируя T
    n
    (x). Благодаря этому свойству многочлен Тейло- ра хорошо приближает функцию f в окрестности точки a. Погрешность приближения составляет
    |f(x) – T
    n
    (x)| = |R
    n
    (x)|, т. е. задавая некоторую точность

    > 0, можно определить окрестность точки a и значение n из усло- вия:
    |R
    n
    (x)| =
    1
    )
    1
    (
    )
    (
    )!
    1
    (
    )
    (




    n
    n
    a
    x
    n
    f

    <

    . (4.2)
    Пример 4.1.
    Найдем приближение функции y = sinx многочленом Тейлора в окрестности точки a = 0.
    Воспользуемся известным выражением для k-ой производной функции sinx:
    (sinx)
    (k)
    = sin x + k
    2

    (4.3)
    Применяя последовательно формулу (4.3), получим:
    f(0) = sin0 = 0;
    f '(0) = cos(0) = 1;
    f"(0) = –sin0 = 0;
    ……………………
    f
    (2k-1)
    (0) = sin (2k – 1)
    2

    = (–1)
    k – 1
    ;
    f
    (2k)
    (0) = 0;
    f
    (2k+1)
    () = (–1)
    k
    cos

    .
    Следовательно, многочлен Тейлора для функции y = sinx для n = 2k имеет вид:
    sinx = x –
    !
    3 3
    x
    + … + (1)
    k – 1
    )!
    1 2
    (
    1 2


    k
    x
    k
    + R
    2k
    (x),
    R
    2k
    (x) = (1)
    k
    1 2
    )!
    1 2
    (
    cos


    k
    x
    k

    .
    Зададим

    = 10
    –4
    и отрезок [
    4

    ,
    4

    ]. Определим n =2k из неравенства:
    |R
    2k
    (x)| =
    1 2
    |
    |
    )!
    1 2
    (
    |
    cos
    |


    k
    x
    k

    <
    1 2
    4
    )!
    1 2
    (
    1








    k
    k

    <
    )!
    1 2
    (
    8 0

    k
    <

    = 10
    -4

    Таким образом, на отрезке
    4

    ,
    4

    функция y = sinx с точностью до

    = 10
    -4
    равна много- члену 5-ой степени:
    sinx = x –
    !
    3 3
    x
    +
    !
    5 5
    x
    = x – 0.1667x
    3
    + 0.0083x
    5
    Пример 4.2.
    Найдем приближение функции y = e
    x
    многочленом Тейлора на отрезке [0, 1] с точностью

    = 10
    –5
    Выберем a = ½, т. е в середине отрезка. При этом величина погрешности в левой части (4.2) принимает минимальное значение. Из математического анализа известно, что для k-ой производной от e
    x
    справедливо равенство:
    (e
    x
    )
    (k)
    = e
    x
    .
    Поэтому
    (e
    a
    )
    (k)
    = e
    a
    = e
    1/2
    ,
    Следовательно, многочлен Тейлора для функции y = e
    x
    имеет вид:
    e
    x
    = e
    1/2 +
    e
    1/2
    (x – ½) +
    !
    2 2
    /
    1
    e
    (x – ½)
    2
    + … +
    !
    2
    /
    1
    n
    e
    (x – ½)
    n
    + R
    n
    (x),
    При этом, учитывая, что x[0, 1], получим оценку погрешности:
    |R
    n
    (x)| <
    )!
    1
    (
    2 1


    n
    e
    n
    . (4.4)
    Составим таблицу погрешностей, вычисленных по формуле (4.4):
    n
    2 3
    4 5
    6
    R
    n
    0.057 0.0071 0.00071 0.000059 0.0000043
    Таким образом, следует взять n = 6.
    4.3. Интерполяция функции многочленами Лагранжа
    Рассмотрим другой подход к приближению функции многочленами. Пусть функция y = f(x)
    определена на отрезке [a, b] и известны значения этой функции в некоторой системе узлов x
    i
     [a,
    b], i = 0, 1, … , n. Например, эти значения получены в эксперименте при наблюдении некоторой величины в определенных точках или в определенные моменты времени x
    0
    , x
    1
    , … , x
    n
    .Обозначим эти значения следующим образом: y
    i
    = f(x
    i
    ), i = 0, 1, … , n. Требуется найти такой многочлен P(x) степени m,
    P(x) = a
    0
    + a
    1
    x + a
    2
    x
    2
    + … + a
    m
    x
    m
    , (4.5) который бы в узлах x
    i
    , i = 0, 1, … , n принимал те же значения, что и исходная функция y = f(x), т. е.
    P(x
    i
    ) = y
    i
    , i = 0, 1, … , n. (4.6)
    Многочлен (4.5), удовлетворяющий условию (4.6), называется интерполяционным много-
    членом.
    Другими словами, ставится задача построения функции y = P(x), график которой проходит через заданные точки (x
    i
    , y
    i
    ), i = 0, 1, … , n (рис. 4.1).

    Рис. 4.1
    Объединяя (4.5) и (4.6), получим:
    a
    0
    + a
    1
    x
    i
    + a
    2
    x
    2
    i
    + … + a
    m
    x
    m
    i
    = y
    i
    ,
    i = 0, 1, … , n. (4.7)
    В искомом многочлене P(x) неизвестными являются m +1 коэффициент a
    0
    , a
    1
    , a
    2
    , …, a
    m
    Поэтому систему (4.7) можно рассматривать как систему из n +1 уравнений с m +1 неизвестными.
    Известно, что для существования единственного решения такой системы необходимо , чтобы вы- полнялось условие: m = n. Таким образом, систему (4.7) можно переписать в развернутом виде:
    a
    0
    + a
    1
    x
    0
    + a
    2
    x
    2 0
    + … + a
    n
    x
    n
    0
    = y
    0
    a
    0
    + a
    1
    x
    1
    + a
    2
    x
    2 1
    + … + a
    n
    x
    n
    1
    = y
    1
    a
    0
    + a
    1
    x
    2
    + a
    2
    x
    2 2
    + … + a
    n
    x
    n
    2
    = y
    2
    (4.8)
    …………………………………………….
    a
    0
    + a
    1
    x
    n
    + a
    2
    x
    2
    n
    + … + a
    n
    x
    n
    n
    = y
    n
    Вопрос о существовании и единственности интерполяционного многочлена решает следу- ющая теорема:
    Теорема 4.1. Существует единственный интерполяционный многочлен степени n, удовле- творяющий условиям (4.6).
    Имеются различные формы записи интерполяционного многочлена. Широко распростра- ненной формой записи является многочлен Лагранжа
    L
    n
    (x) =






    i
    j
    j
    i
    j
    n
    i
    i
    x
    x
    x
    x
    y
    1
    =
    )
    )...(
    )(
    )...(
    (
    )
    )...(
    )(
    )...(
    (
    1 1
    0 1
    1 0
    1
    n
    i
    i
    i
    i
    i
    i
    n
    i
    i
    n
    i
    i
    x
    x
    x
    x
    x
    x
    x
    x
    x
    x
    x
    x
    x
    x
    x
    x
    y














    . (4.9)
    В частности, для линейной и квадратичной интерполяции по Лагранжу получим следующие интерполяционные многочлены:
    L
    1
    (x) = y
    0
    )
    (
    )
    (
    1 0
    1
    x
    x
    x
    x


    + y
    1
    )
    (
    )
    (
    0 1
    0
    x
    x
    x
    x


    ,
    L
    2
    (x) = y
    0
    )
    )(
    (
    )
    )(
    (
    2 0
    1 0
    2 1
    x
    x
    x
    x
    x
    x
    x
    x




    +
    y
    1
    )
    )(
    (
    )
    )(
    (
    2 1
    0 1
    2 0
    x
    x
    x
    x
    x
    x
    x
    x




    + y
    2
    )
    )(
    (
    )
    )(
    (
    1 2
    0 2
    1 0
    x
    x
    x
    x
    x
    x
    x
    x




    .
    Пример 4.3.
    Построим интерполяционный многочлен Лагранжа по следующим данным:
    x 0 2 3 5
    y 1 3 2 5

    Степень многочлена Лагранжа для n +1 узла равна n. Для нашего примера многочлен Ла- гранжа имеет третью степень. В соответствии с (4.9)
    L
    3
    (x) = 1
    )
    5 0
    )(
    3 0
    )(
    2 0
    (
    )
    5
    )(
    3
    )(
    2
    (







    x
    x
    x
    + 3
    )
    5 2
    )(
    3 2
    (
    2
    )
    5
    )(
    3
    (





    x
    x
    x
    + 2
    )
    5
    )(#
    2 3
    (
    3
    )
    5
    )(
    2
    (





    x
    x
    x
    +
    5
    )
    3 5
    )(
    2 5
    (
    5
    )
    3
    )(
    2
    (





    x
    x
    x
    = 1 +
    15 62
    x –
    6 13
    x
    2
    +
    10 3
    x
    3
    .
    Пример 4.4.
    Рассмотрим пример использования интерполяционного многочлена Лагранжа для вычисле- ния значения заданной функции в промежуточной точке. Эта задача возникает, например, когда заданы табличные значения функции с крупным шагом, а требуется составить таблицу значений с маленьким шагом.
    Для функции y = sinx известны следующие данные.
    x
    0
    /6
    /3
    /2
    y
    0
    ½
    2 3
    1
    Вычислим y(0.25).
    Найдем многочлен Лагранжа третьей степени:
    L
    3
    (x) = 0
    )
    6 0
    )(
    3 0
    )(
    6 0
    (
    )
    2
    )(
    3
    )(
    6
    (













    x
    x
    x
    +
    2 1
    )
    2 6
    )(
    3 6
    )(
    0 6
    (
    )
    2
    )(
    3
    )(
    0
    (














    x
    x
    x
    +
    2 3
    )
    2 3
    )(
    6 3
    )(
    0 3
    (
    )
    2
    )(
    6
    )(
    0
    (














    x
    x
    x
    + 1
    )
    3 2
    )(
    6 2
    )(
    0 2
    (
    )
    3
    )(
    6
    )(
    0
    (














    x
    x
    x
    При x = 0.25 получим y(0.25) = sin 0.25  0.249.
    Погрешность интерполяции. Пусть интерполяционный многочлен Лагранжа построен для известной функции f(x). Необходимо выяснить, насколько этот многочлен близок к функции в точ- ках отрезка [a, b], отличных от узлов. Погрешность интерполяции равна |f(x) – P
    n
    (x)|. Оценку по- грешности можно получить на основании следующей теоремы.
    Теорема 4.2. Пусть функция f(x)дифференцируема n +1 раз на отрезке [a, b], содержащем узлы интерполяции x
    i
     [a, b], i = 0, 1, … , n. Тогда для погрешности интерполяции в точке x  [a,
    b] справедлива оценка:
    |f(x) – L
    n
    (x)|
    )!
    1
    (
    1


    n
    M
    n
    |

    n+1
    (x)|, (4.10) где
    M
    n+1
    =
     
    b
    a,
    max
    |f
    (n+1)
    (x)|,

    n+1
    (x) = (x – x
    0
    )(x – x
    1
    )…. (x – x
    n
    ).
    Для максимальной погрешности интерполяции на всем отрезке [a, b] справедлива оценка:
     
    b
    a,
    max
    |f(x) – L
    n
    (x)| 
    )!
    1
    (
    1


    n
    M
    n
     
    b
    a,
    max
    |

    n
    (x)| (4.11)
    Пример 4.5.
    Оценим погрешность приближения функции f(x) =
    x
    в точке x = 116 и на всем отрезке
    [a, b], где a = 100, b = 144, с помощью интерполяционного много члена Лагранжа L
    2
    (x) второй сте- пени, построенного с узлами x
    0
    = 100, x
    2
    = 144.
    Найдем первую, вторую и третью производные функции f(x):

    f '(x)=
    2 1
    x
    – 1/2
    , f "(x)= –
    4 1
    x
    –3/2
    ,
    f'''(x)=
    8 3
    x
    –5/2
    M
    3
    =
     
    b
    a,
    max
    | f'''(x)| =
    8 3
    100
    –5/2
    =
    8 3
    10
    –5
    В соответствии с (4.9) получим оценку погрешности в точке x = 116:
    |
    116
    L
    2
    (116)| 
    !
    3 1
    |(116 – 100)(116 – 121)(116 – 144)| =
    16 1
    10
    –5
    16528 = 1.410
    – 3
    Оценим погрешность приближения функции f(x) =
    x
    на всем отрезке в соответствии с
    (4.11):
     
    b
    a,
    max
    |
    x
    – L
    2
    (x)| 
    16 10 5

     
    b
    a,
    max
    |(x – 100)(x – 121)(x –144)|  2.510
    –3 4.4. Аппроксимация функций. Метод наименьших квадратов
    В инженерной деятельности часто возникает необходимость описать в виде функциональ- ной зависимости связь между величинами, заданными таблично или в виде набора точек с коорди- натами (x
    i
    , y
    i
    ), i = 0, 1, 2,... , n, где n – общее количество точек. Как правило, эти табличные данные получены экспериментально и имеют погрешности (рис. 2.5)
    Рис.4.2
    При аппроксимации желательно получить относительно простую функциональную зависи- мость (например, многочлен), которая позволила бы "сгладить" экспериментальные погрешности, вычислять значения функции в точках, не содержащихся в исходной таблице.
    Эта функциональная зависимость должна с достаточной точностью соответствовать исход- ной табличной зависимости. В качестве критерия точности чаще всего используют критерий
    наименьших квадратов, т.е. определяют такую функциональную зависимость f(x), при которой
    S =



    n
    i
    i
    i
    f
    y
    0 2
    )
    (
    , (4.12) обращается в минимум.
    Погрешность приближения оценивается величиной среднеквадратического уклонения

    =
    S
    n 1 1

    . (4.13)
    В качестве функциональной зависимости рассмотрим многочлен
    P
    m
    (x)=a
    0
    +
    a
    1
    x
    +
    a
    2
    x
    2
    +...+a
    m
    x
    m
    . (4.14)
    Формула (4.12) примет вид

    S =



    n
    i
    i
    m
    i
    x
    P
    y
    0 2
    ))
    (
    (
    Условия минимума S можно записать, приравнивая нулю частные производные S по всем переменным a
    0
    ,
    a
    1
    ,
    a
    2
    , … , a
    m
    . Получим систему уравнений
    k
    дa
    дS
    = –






    n
    i
    k
    i
    m
    i
    m
    i
    x
    x
    a
    x
    a
    a
    y
    0 1
    1 0
    )
    (
    2
    = 0, или






    n
    i
    k
    i
    m
    i
    m
    i
    x
    x
    a
    x
    a
    a
    y
    0 1
    1 0
    )
    (
    = 0, k = 0, 1, … , m. (4.15)
    Систему уравнений (4.15) перепишем в следующем виде:
    a
    0


    n
    i
    k
    i
    x
    0
    + a
    1



    n
    i
    k
    i
    x
    0 1
    + … +a
    m



    n
    i
    m
    k
    i
    x
    0
    =


    n
    i
    k
    i
    i
    x
    y
    0
    , k = 0, 1, … , m (4.16)
    Введем обозначения:
    c
    k
    =


    n
    i
    k
    i
    x
    0
    , b
    k
    =


    n
    i
    k
    i
    i
    x
    y
    0
    Система (4.16) может быть записана так:
    a
    0
    c
    k
    + a
    1
    c
    k+1
    + … + c
    k+m
    a m
    = b
    k
    , k = 0, 1, … , m. (4.17)
    Перепишем систему (4.17) в развернутом виде:
    c
    0
    a
    0
    + c
    1
    a
    1
    + c
    2
    a
    2
    … + c
    m
    a
    m
    = b
    0
    c
    1
    a
    0
    + c
    2
    a
    1
    + c
    3
    a
    2
    … + c
    m+1
    a
    m
    = b
    1
    …………………………………….. (4.18)
    c
    m
    a
    0
    + c
    m+1
    a
    1
    + c
    m+2
    a
    2
    … + c
    2m
    a
    m
    = b
    m
    Матричная запись системы (4.18) имеет следующий вид:
    Ca = b. (4.19)
    Для определения коэффициентов a
    k
    , k = 0, 1, … , m, и, следовательно, искомого многочлена
    (4.14) необходимо вычислить суммы c
    k
    , b
    k
    и решить систему уравнений (4.18). Матрица C системы
    (4.19) называется матрицей Грама и является симметричной и положительно определенной. Эти полезные свойства используются при решении.
    Погрешность приближения в соответствии с формулой (4.13) составит

    =






    n
    i
    i
    m
    i
    x
    P
    y
    n
    0 2
    )
    (
    1 1
    . (4.20)
    Рассмотрим частные случаи m =1 и m = 2.
    1. Линейная аппроксимация (m = 1).
    P
    1
    (x) = a
    0
    + a
    1
    x.
    c
    0
    =


    n
    i
    i
    x
    0 0
    = n + 1; c
    1
    =


    n
    i
    i
    x
    0 1
    =


    n
    i
    i
    x
    0
    ; c
    2
    =


    n
    i
    i
    x
    0 2
    ; (4.21)
    b
    0
    =


    n
    i
    i
    i
    x
    y
    0 0
    =


    n
    i
    i
    y
    0
    ; b
    1
    =


    n
    i
    i
    i
    x
    y
    0 1
    =


    n
    i
    i
    i
    x
    y
    0
    . (4.22)
    c
    0
    c
    1
    n+1


    n
    i
    i
    x
    0
    C = = ,
    c
    1
    c
    2


    n
    i
    i
    x
    0


    n
    i
    i
    x
    0 2

    b = (b
    0
    , b
    1
    )
    T
    = (


    n
    i
    i
    y
    0
    ,


    n
    i
    i
    i
    x
    y
    0
    )
    T
    Решение системы уравнений Ca = b найдем по правилу Крамера:
    a
    0
    =
    C
    C
    1
    , a
    1
    =
    C
    C
    2
    , где C – определитель матрицы C, аC
    i
     – определитель матрицы C
    i
    , полученной из мат- рицы C заменой i-го столбца столбцом свободных членов b, i = 1, 2.
    Таким образом,
    a
    0
    =
    2 1
    2 0
    1 1
    2 0
    c
    c
    c
    c
    b
    c
    b


    , a
    1
    =
    2 1
    2 0
    1 0
    0 1
    c
    c
    c
    c
    b
    c
    b


    . (4.23)
    Алгоритм 4.1 (Алгоритм метода наименьших квадратов. Линейная аппроксимация).
    Шаг 1. Ввести исходные данные: x
    i
    , y
    i
    , i=0, 1, 2, ... , n.
    Шаг 2. Вычислить коэффициенты c
    0
    , c
    1
    , b
    0
    , b
    1 по формулам (4.21), (4.22).
    Шаг 3. Вычислить a
    0
    , a
    1
    по формулам (4.23).
    Шаг 4. Вычислить величину погрешности

    1
    =







    n
    i
    i
    i
    x
    a
    a
    y
    n
    0 2
    1 0
    )
    (
    1 1
    . (4.24)
    Шаг 5. Вывести на экран результаты: аппроксимирующую линейную функцию P
    1
    (x) = a
    0
    +
    a
    1
    x и величину погрешности

    1 2. Квадратичная аппроксимация (m = 2).
    P
    2
    (x) = a
    0
    + a
    1
    x + a
    2
    x
    2
    c
    0
    =


    n
    i
    i
    x
    0 0
    = n+1; c
    1
    =


    n
    i
    i
    x
    0 1
    =


    n
    i
    i
    x
    0
    ; c
    2
    =


    n
    i
    i
    x
    0 2
    ; c
    3
    =


    n
    i
    i
    x
    0 3
    ; c
    4
    =


    n
    i
    i
    x
    0 4
    . (4.25)
    b
    0
    =


    n
    i
    i
    i
    x
    y
    0 0
    =


    n
    i
    i
    y
    0
    ; b
    1
    =


    n
    i
    i
    i
    x
    y
    0 1
    =


    n
    i
    i
    i
    x
    y
    0
    ; b
    2
    =


    n
    i
    i
    i
    x
    y
    0 2
    . (4.26)
    c
    0
    c
    1
    c
    2
    C = c
    1
    c
    2
    c
    3
    c
    2
    c
    3
    c
    4
    b = (b
    0
    , b
    1
    , b
    2
    )
    T
    Решение системы уравнений Ca = b найдем по правилу Крамера:
    a
    i
    =
    C
    C
    i
    , i = 0, 1, где C – определитель матрицы C, аC
    i
     – определитель матрицы C
    i
    , полученной из мат- рицы C заменой i-го столбца столбцом свободных членов b.
    C = c
    0
    c
    2
    c
    4
    + 2c
    1
    c
    2
    c
    3
    c
    3 2
    – с
    2 1
    c
    4
    c
    2 3
    c
    0
    . (4.27)
    b
    0
    c
    1
    c
    2
    C
    1
     = b
    1
    c
    2
    c
    3
    = b
    0
    c
    2
    c
    4
    + b
    2
    c
    1
    c
    3
    + b
    1
    c
    2
    c
    3
    b
    2
    c
    2 2
    b
    1
    c
    1
    c
    4
    b
    0
    c
    2 3
    . (4.28)
    b
    2
    c
    3
    c
    4
    c
    0
    b
    0
    c
    2
    C
    2
     = c
    1
    b
    1
    c
    3
    = b
    1
    c
    0
    c
    4
    + b
    0
    c
    2
    c
    3
    + b
    2
    c
    1
    c
    2
    b
    1
    c
    2 2
    b
    0
    c
    1
    c
    4
    b
    2
    c
    0
    c
    3
    . (4.29)
    c
    2
    b
    2
    c
    4

    c
    0
    c
    1
    b
    0
    C
    3
     = c
    1
    c
    2
    b
    1
    = b
    2
    c
    0
    c
    2
    + b
    1
    c
    1
    c
    2
    + b
    0
    c
    1
    c
    3
    b
    0
    c
    2 2
    b
    2
    c
    2 1
    b
    1
    c
    0
    c
    3
    . (4.30)
    c
    2
    c
    3
    b
    2
    a
    0
    =
    C
    C
    1
    , a
    1
    =
    C
    C
    2
    , a
    2
    =
    C
    C
    3
    . (4.31)
    Алгоритм 4.2 (Алгоритм метода наименьших квадратов. Квадратичная аппроксимация).
    Шаг 1. Ввести исходные данные: x
    i
    , y
    i
    , i=0, 1, 2, ... , n.
    Шаг 2. Вычислить коэффициенты c
    0
    , c
    1
    , c
    2
    , c
    3
    , c
    4
    , b
    0
    , b
    1
    , b
    2
    ,
    по формулам (4.25), (4.26).
    Шаг 3. Вычислить C, C
    1
    , C
    2
    , C
    3
     по формулам (4.27) – (4.30).
    Шаг 4. Вычислить a
    0
    , a
    1
    , a
    2
    по формулам (4.31).
    Шаг 5. Вычислить величину погрешности

    2
    =








    n
    i
    i
    i
    i
    x
    a
    x
    a
    a
    y
    n
    0 2
    2 2
    1 0
    )
    (
    1 1
    . (4.32)
    Шаг 5. Вывести на экран результаты : аппроксимирующую квадратичную функцию P
    2
    (x) =
    a
    0
    + a
    1
    x + a
    2
    x
    2
    и величину погрешности

    2
    Пример 4.6.
    Построим по методу наименьших квадратов многочлены первой и второй степени и оценим степень приближения. Значения y
    i в точках x
    i
    , i =0, 1, 2, 3, 4 приведены в таблице 2.3.
    Таблица 4.1
    i
    0 1
    2 3
    4
    x
    i
    1 2
    3 4
    5
    y
    i
    –1 1
    2 4
    6
    Вычислим коэффициенты c
    0
    , c
    1
    , c
    2
    , c
    3
    , c
    4
    , b
    0
    , b
    1
    , b
    2
    ,
    по формулам (4.25), (4.26):
    c
    0
    = 5; c
    1
    = 15; c
    2
    = 55; c
    3
    = 225; c
    4
    = 979;
    b
    0
    = 12; b
    1
    = 53; b
    2
    = 235.
    1. Линейная аппроксимация (m =1).
    Система уравнений для определения коэффициентов a
    0
    и a
    1 многочлена первой степени
    P
    2
    (x) = a
    0
    + a
    1
    x + a
    2
    x
    2
    имеет вид
    5a
    0
    + 15a
    1
    = 12 15a
    0
    + 55a
    1
    = 53
    По формулам (4.23) найдем коэффициенты a
    0
    и a
    1
    :
    a
    0
    =
    2 1
    2 0
    1 1
    2 0
    c
    c
    c
    c
    b
    c
    b


     –2.7, a
    1
    =
    2 1
    2 0
    1 0
    0 1
    c
    c
    c
    c
    b
    c
    b


     1.7.
    P
    1
    (x) = a
    0
    + a
    1
    x = –2.7 + 1.7x.
    2. Квадратичная аппроксимация (m =2).
    Система уравнений для определения коэффициентов a
    0
    , a
    1
    и a
    2 многочлена второй степени
    P
    2
    (x) = a
    0
    + a
    1
    x + a
    2
    x
    2
    имеет вид
    5a
    0
    + 15a
    1
    + 55a
    2
    = 12 15a
    0
    + 55a
    1
    + 225a
    2
    = 53 55a
    0
    + 225a
    1
    + 979a
    2
    = 235
    По формулам (4.31) найдем коэффициенты a
    0
    , a
    1
    и a
    2
    :
    a
    0
     –2.20, a
    1
     1.27, a
    2
     0.07.
    P
    2
    (x) = a
    0
    + a
    1
    x + a
    2
    x
    2
    = –2.20 + 1.27x + 0.07x
    2
    Сравним значения, рассчитанные для функциональной зависимости, с исходными данны- ми. Результаты приведены в табл.2.4.

    Таблица 4.2
    i
    0 1
    2 3
    4
    x
    i
    1 2
    3 4
    5
    y
    i
    –1 1
    2 4
    6
    P
    1
    (x
    i
    )
    –1 0.7 2.4 4.1 5.8
    P
    2
    (x
    i
    )
    –1 0.62 2.24 4
    6.9
    Погрешность приближения в соответствии с формулами (4.24) и (4.32) составит

    1
    =





    4 0
    2 1
    )
    (
    5 1
    i
    i
    i
    x
    P
    y
    = 0.245.

    2
    =





    4 0
    2 2
    )
    (
    5 1
    i
    i
    i
    x
    P
    y
    = 0.084.
    Тема 5. Численное интегрирование функций одной переменной
    5.1. Постановка задачи численного интегрирования
    Далеко не все интегралы можно вычислить по известной из математического анализа фор- муле Ньютона – Лейбница:
    I =

    b
    a
    dx
    x
    f
    )
    (
    = F(b) – F(a), (5.1) где F(x) – первообразная функции f(x). Например, в элементарных функциях не выражается инте- грал
    dx
    e
    b
    x


    0 2
    . Но даже в тех случаях, когда удается выразить первообразную функцию F(x) через элементарные функции, она может оказаться очень сложной для вычислений. Кроме того, точное значение интеграла по формуле (5.1) нельзя получить, если функция f(x) задается таблицей. В этих случаях обращаются к методам численного интегрирования.
    Суть численного интегрирования заключается в том, что подынтегральную функцию f(x) заменяют другой приближенной функцией, так, чтобы, во-первых, она была близка к f(x) и, во вто- рых, интеграл от нее легко вычислялся. Например, можно заменить подынтегральную функцию интерполяционным многочленом. Широко используют квадратурные формулы:

    b
    a
    dx
    x
    f
    )
    (



    n
    i
    i
    i
    x
    f
    A
    0
    )
    (
    , (5.2) где x
    i
    – некоторые точки на отрезке [a, b],называемые узлами квадратурной формулы, A
    i
    – числовые коэффициенты, называемые весами квадратурной формулы, n  0 – целое число.
    5.2. Метод прямоугольников
    Формулу прямоугольников можно получить из геометрической интерпретации интеграла.
    Будем интерпретировать интеграл

    b
    a
    dx
    x
    f
    )
    (
    как площадь криволинейной трапеции, ограниченной графиком функции y = f(x), осью абсцисс и прямыми x = a и x = b (рис. 5.1).

    Рис. 5.1
    Разобьем отрезок [a, b] на n равных частей длиной h, так, что h =
    n
    a
    b
    . При этом получим точки a = x
    0
    < x
    1
    < x
    2
    < … < x
    n
    = b и x
    i+1
    = x
    i
    + h, i = 0, 1, … , n – 1 (рис. 5.2)
    Рис. 5.2
    Заменим приближенно площадь криволинейной трапеции площадью ступенчатой фигуры, изображенной на рис. 5.3.
    Рис. 5.3
    Эта фигура состоит из n прямоугольников. Основание i-го прямоугольника образует отре- зок [x
    i
    , x
    i+1
    ] длины h, а высота основания равна значению функции в середине отрезка [x
    i
    , x
    i+1
    ], т е.
    f








    2 1
    i
    i
    x
    x
    (рис. 5.4).

    Рис. 5.4
    Тогда получим квадратурную формулу средних прямоугольников:
    I =

    b
    a
    dx
    x
    f
    )
    (
    I
    пр
    =









     

    1 0
    1 2
    n
    i
    i
    i
    x
    x
    f
    n
    a
    b
    (5.3)
    Формулу (5.3) называют также формулой средних прямоугольников. Иногда используют формулы
    I I
    л
    пр
    =




    1 0
    )
    (
    n
    i
    i
    x
    f
    n
    a
    b
    , (5.4)
    I I
    п
    пр
    =



    n
    i
    i
    x
    f
    n
    a
    b
    1
    )
    (
    , (5.5) которые называют соответственно квадратурными формулами левых и правых прямоугольников.
    Геометрические иллюстрации этих формул приведены на рис. 5.5 и 5.6.
    Рис. 5.5

    Рис. 5. 6
    Оценка погрешности. Для оценки погрешности формулы прямоугольников воспользуемся следующей теоремой .
    Теорема 5.1. Пусть функция f дважды непрерывно дифференцируема на отрезке [a, b]. То- гда для формулы прямоугольников справедлива следующая оценка погрешности:
    | II
    пр
    | 
    24
    )
    (
    2
    a
    b
    M

    h
    2
    , (5.6) где M
    2
    =
     
    b
    a,
    max
    |f "(x)|
    Пример 5.1.
    Вычислим значение интеграла
    dx
    e
    x


    1 0
    2
    по формуле средних прямоугольников (5.3) с ша- гом h = 0.1.
    Составим таблицу значений функции e
    2
    x

    (табл. 5.1):
    Таблица 5.1
    x
    e
    2
    x

    x
    e
    2
    x

    0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 1.0000000 0.9975031 0.9900498 0.9777512 0.9607894 0.9394131 0.9139312 0.8847059 0.8521438 0.8166865 0.7788008 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 0.7389685 0.6976763 0.6554063 0.6126264 0.5697828 0.5272924 0.4855369 0.4448581 0.4055545 0.3678794
    Производя вычисления по формуле (5.3), получим:
    I
    пр
    = 0.74713088.
    Оценим погрешность полученного значения. Имеем:
    f "(x) = (e
    2
    x

    )" = (4x
    2
    – 2) e
    2
    x

    . Нетрудно убедиться, что | f "(x)|  M
    2
    = 2. Поэтому по фор- муле(5.4)
    | II
    пр
    | 
    24 1
    2 
    (0.1)
    2
     0.84 10
    -3

    5.3. Метод трапеций
    Выведем формулу трапеций так же, как и формулу прямоугольников, из геометрических соображений. Заменим график функции y = f(x) (рис.5.1) ломаной линией (рис.5.7), полученной сле- дующим образом. Из точек a = x
    0
    , x
    1
    , x
    2
    ,…, x
    n
    = b проведем ординаты до пересечения с кривой y =
    f(x). Концы ординат соединим прямолинейными отрезками.
    Рис. 5.7
    Тогда площадь криволинейной трапеции приближенно можно считать равной площади фи- гуры, составленной из трапеций. Так как площадь трапеции, построенной на отрезке [x
    i
    , x
    i+1
    ] длины
    h =
    n
    a
    b
    , равна h
    2
    )
    (
    )
    (
    1


    i
    i
    x
    f
    x
    f
    , то, пользуясь этой формулой для i = 0, 2, … , n – 1, получим
    квадратурную формулу трапеций:
    I =

    b
    a
    dx
    x
    f
    )
    (
    I
    тр
    = h













    2
    )
    (
    )
    (
    2
    )
    (
    )
    (
    2
    )
    (
    (
    1 2
    1 1
    0
    n
    n
    x
    f
    x
    f
    x
    f
    x
    f
    x
    f
    x
    f
    =












    1 1
    0
    )
    (
    2
    )
    (
    )
    (
    n
    i
    i
    n
    x
    f
    x
    f
    x
    f
    n
    a
    b
    (5.7)
    Оценка погрешности. Для оценки погрешности формулы трапеций воспользуемся следу- ющей теоремой.
    Теорема 5.2. Пусть функция f дважды непрерывно дифференцируема на отрезке [a, b]. То- гда для формулы трапеций справедлива следующая оценка погрешности:
    | II
    тр
    | 
    12
    )
    (
    2
    a
    b
    M

    h
    2
    , (5.8) где M
    2
    =
     
    b
    a,
    max
    |f "(x)|.
    Пример 5.2.
    Вычислим значение интеграла
    dx
    e
    x


    1 0
    2
    по формуле трапеций (5.7) и сравним полученный результат с результатом примера 5.1.
    Используя таблицу значений функции e
    2
    x

    из примера 5.1 и производя вычисления по фор- муле трапеций (5.7), получим:
    I
    тр
    = 0.74621079.
    Оценим погрешность полученного значения. В примере (5.1) получили оценку: | f "(x)|  M
    2
    = 2. Поэтому по формуле (5.8)
    | II
    тр
    | 
    12 1
    2 
    (0.1)
    2
     1.7 10
    -3

    Сравнивая результаты примеров 5.1 и 5.2, видим, что метод средних прямоугольников име- ет меньшую погрешность, т.е. он более точный.
    5.4. Метод Симпсона (метод парабол)
    Заменим график функции y = f(x) на отрезке [x
    i
    , x
    i+1
    ], i = 0, 2, … , n – 1, параболой, прове- денной через точки (x
    i
    , f(x
    i
    )), (x
    '
    i
    ,f(x
    '
    i
    )), (x
    i+1
    , f(x
    i+1
    )), где x
    '
    i
    - середина отрезка [x
    i
    , x
    i+1
    ]. Эта парабола есть интерполяционный многочлен второй степени L
    2
    (x) с узлами x
    i
    , x
    '
    i
    , x
    i+1
    . Нетрудно убедиться, что уравнение этой параболы имеет вид:
    y = L
    2
    (x) =
    f(x
    '
    i
    ) +
    h
    x
    f
    x
    f
    i
    i
    )
    (
    )
    (
    1


    (xx
    '
    i
    ) +
    2
    /
    )
    (
    )
    (
    2
    )
    (
    2
    '
    1
    h
    x
    f
    x
    f
    x
    f
    i
    i
    i



    (x - x
    '
    i
    )
    2
    , (5.9) где h =
    n
    a
    b
    Проинтегрировав функцию (5.9) на отрезке [x
    i
    , x
    i+1
    ], получим
    I
    i
    =

    1
    )
    (
    i
    i
    x
    x
    dx
    x
    f


    1
    )
    (
    2
    i
    i
    x
    x
    dx
    x
    L
    =
    6
    h
    ( f(x
    i
    ) + 4f(x
    '
    i
    ) + f(x
    i+1
    )). (5.10)
    Суммируя выражение (5.10) по i = 0, 1, 2, … , n – 1, получим квадратурную формулу Симп-
    сона (или формулу парабол):
    I =

    b
    a
    dx
    x
    f
    )
    (
    I
    С
    =
    6
    h
    ( f(x
    0
    ) + f(x
    n
    ) + 4



    1 0
    '
    )
    (
    n
    i
    i
    x
    f
    + 2



    1 1
    )
    (
    n
    i
    i
    x
    f
    ). (5.11)
    Оценка погрешности. Для оценки погрешности формулы Симпсона воспользуемся следу- ющей теоремой.
    Теорема 5.2. Пусть функция f имеет на отрезке [a, b] непрерывную производную четвертого порядка f
    (4)
    (x). Тогда для формулы Симпсона (5.9) справедлива следующая оценка погрешности:
    | II
    С
    | 
    2880
    )
    (
    4
    a
    b
    M

    h
    4
    , (5.12) где M
    4
    =
     
    b
    a,
    max
    | f
    (4)
    (x)|.
    Замечание. Если число элементарных отрезков, на которые делится отрезок [a, b], четно , т.е. n = 2m, то параболы можно проводить через узлы с целыми индексами, и вместо элементарного отрезка [x
    i
    , x
    i+1
    ] длины h рассматривать отрезок [x
    2i
    , x
    2i+2
    ] длины 2h. Тогда формула Симпсона при- мет вид:
    I
    3
    h
    (f(x
    0
    ) + f(x
    2m
    ) + 4



    m
    i
    i
    x
    f
    1 1
    2
    )
    (
    + 2



    1 1
    2
    )
    (
    m
    i
    i
    x
    f
    ), (5.13) а вместо оценки (5.10) будет справедлива следующая оценка погрешности:
    | II
    С
    | 
    180
    )
    (
    4
    a
    b
    M

    h
    4
    , (5.14)
    Пример 5.3.
    Вычислим значение интеграла
    dx
    e
    x


    1 0
    2
    по формуле Симпсона (5.11) и сравним получен- ный результат с результатами примеров 5.1 и 5.2.
    Используя таблицу значений функции e
    2
    x

    из примера 5.1 и производя вычисления по фор- муле Симпсона (5.11) , получим:
    I
    С
    = 0.74682418.

    Оценим погрешность полученного значения. Вычислим четвертую производную f
    (4)
    (x).
    f
    (4)
    (x) = (16x
    4
    – 48x
    2
    + 12) e
    2
    x

    , | f
    (4)
    (x)|  12.
    Поэтому
    | II
    С
    | 
    2880 1
    12 
    (0.1)
    4
     0.42  10
    -6
    Сравнивая результаты примеров 5.1, 5.2 и 5.3, видим , что метод Симпсона имеет меньшую погрешность, чем метод средних прямоугольников и метод трапеций.
    5.5. Правило Рунге практической оценки погрешности
    Оценки погрешности по формулам (5.4), (5.8) и (5.12) являются априорными. Они зависят от длины элементарного отрезка h, и при достаточно малом h справедливо приближенное равен- ство:
    II
    h
    Ch
    k
    , (5.15) где I
    h
    приближенное значение интеграла, вычисленное по одной из формул (5.3), (5.5), (5.9), C  0 и
    k > 0 – величины, не зависящие от h.
    Если уменьшить шаг h в два раза, то, в соответствии с (5.15) получим:
    II
    h/2

    k
    2 1
    Ch
    k

    k
    2 1
    ( II
    h
    ). (5.16)
    Непосредственное использование оценок погрешности (5.4), (5.8) и (5.12) неудобно, так как при этом требуется вычисление производных функции f
    (x). В вычислительной практике использу- ются другие оценки.
    Вычтем из равенства (5.15) равенство (5.16):
    I
    h/2
    I
    h

    k
    2 1
    Ch
    k
    (2
    k
    – 1). (5.17)
    Учитывая приближенное равенство (5.16), получим следующее приближенное равенство:
    II
    h/2

    1 2
    2
    /


    k
    h
    h
    I
    I
    . (5.18)
    Приближенное равенство (5.18) дает апостериорную оценку погрешности. Вычисление этой оценки называется правилом Рунге. Правило Рунге – это эмпирический способ оценки погрешно- сти, основанный на сравнении результатов вычислений , проводимых с разными шагами h.
    Для формул прямоугольников и трапеций k = 2, а для формулы Симпсона k = 4. Поэтому для этих формул приближенное равенство (5.18) принимает вид:
    II
    пр

    )
    (
    3 1
    2
    /
    h
    пр
    h
    пр
    I
    I

    , (5.19)
    II
    тр

    )
    (
    3 1
    2
    /
    h
    тр
    h
    тр
    I
    I

    , (5.20)
    II
    С

    )
    (
    15 1
    2
    /
    h
    С
    h
    С
    I
    I

    . (5.21)
    Используя правило Рунге, можно построить процедуру приближенного вычисления инте- грала с заданной точностью

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

    .
    Пример 5.4.

    Найдем значение интеграла
    dx
    e
    x


    1 0
    2
    с точностью

    = 10
    -4
    , используя формулу трапеций и применяя вышеизложенную процедуру дробления шага. В примере 5.2 было получено значение I
    1
    h
    при h
    1
    = 0.1, I
    h
    =0.74621079. Уменьшим шаг вдвое: h
    2
    = 0.05 и вычислим I
    2
    h
    = 0.74667084,

    2
    =
    3 1
    (
    I
    2
    h
    - I
    1
    h
    ) =
    3 1
    (0.74667084 – 0.74621079)  1.510
    -4
    . Так как |

    2
    | >

    , то снова дробим шаг: h
    3
    = 0.025, вычисляем I
    3
    h
    = 0.74678581,

    2
    =
    3 1
    ( I
    3
    h
    - I
    2
    h
    ) =
    3 1
    (0.74678581 – 0.74667084)  410
    -5
    . Поскольку |

    3
    |
    <

    , требуемая точность достигнута и I  0.7468  0.0001.
    1   2   3   4   5


    написать администратору сайта