Реферат: Численные методы решения типовых математических задач

--PAGE_BREAK--1. Решение систем линейных алгебраических уравнений методом простой итерации 1.1 Постановка задачи


Разработать схему алгоритма и написать программу на языке TurboPascal7.0 для решении систем линейных алгебраических уравнений, используя метод простой итерации.
1.2 Математическая формулировка задачи


Пусть А – невырожденная матрица <img width=«44» height=«17» src=«ref-1_1289180173-269.coolpic» alt=«image036001» v:shapes=«Рисунок_x0020_1»>и нужно решить систему
<img width=«231» height=«130» src=«ref-1_1289180442-1369.coolpic» v:shapes=«Рисунок_x0020_2»>
где диагональные элементы матрицы А ненулевые.


1.3 Обзор существующих численных методов решения задачи
Метод Гаусса


В методе Гаусса матрица СЛАУ с помощью равносильных преобразований преобразуется в верхнюю треугольную матрицу, получающуюся в результате прямого хода. В обратном ходе определяются неизвестные.

Пусть дана СЛАУ

Запишем расширенную матрицу системы:


<img width=«544» height=«222» src=«ref-1_1289181811-10412.coolpic» v:shapes="_x0000_i1027">
На первом шаге алгоритма Гаусса выберем диагональный элемент (если он равен 0, то первую строку переставляем с какой-либо нижележащей строкой) и объявляем его a11≠0 ведущим, а соответствующую строку и столбец, на пересечении которых он стоит — ведущими. Обнулим элементы ведущего столбца. Для этого сформируем числа (-a22/a11), (-a31/a11),…, (an1/a11).

LU-разложения матриц

Компьютерная реализация метода Гаусса часто осуществляется с использованием LU-разложения матриц.

LU – разложение матрицы A представляет собой разложение матрицы A в произведение нижней и верхней треугольных матриц, т.е.
A=LU
где L — нижняя треугольная матрица (матрица, у которой все элементы, находящиеся выше главной диагонали равны нулю, lij=0 при i>j), U- верхняя треугольная матрица (матрица, у которой все элементы, находящиеся ниже главной диагонали равны нулю, uij=0 при i>j).

LU – разложение может быть построено с использованием описанного выше метода Гаусса. Рассмотрим k — ый шаг метода Гаусса, на котором осуществляется обнуление поддиагональных элементов k — го столбца матрицы. Как было описано выше, с этой целью используется следующая операция

<img width=«563» height=«73» src=«ref-1_1289192223-4620.coolpic» v:shapes="_x0000_i1028">
В терминах матричных операций такая операция эквивалентна умножению A(k)=MkA(k-1), где элементы матрицы определяются следующим образом

В терминах матричных операций такая операция эквивалентна умножению A(k)=MkA(k-1), где элементы матрицы определяются следующим образом
<img width=«312» height=«124» src=«ref-1_1289196843-34212.coolpic» v:shapes="_x0000_i1029">
В результате прямого хода метода Гаусса получим, A(n-1)=U
<img width=«519» height=«53» src=«ref-1_1289231055-3822.coolpic» v:shapes="_x0000_i1030">

<img width=«467» height=«201» src=«ref-1_1289234877-7904.coolpic» v:shapes="_x0000_i1031">
где A(n-1)=U — верхняя треугольная матрица, а — нижняя треугольная матрица, имеющая вид.

Таким образом, искомое разложение A=LU получено.

Метод прогонки


Метод прогонки является одним из эффективных методов решения СЛАУ с трех — диагональными матрицами, возникающих при конечно-разностной аппроксимации задач для обыкновенных дифференциальных уравнений (ОДУ) и уравнений в частных производных второго порядка и является частным случаем метода Гаусса. Рассмотрим следующую СЛАУ:
<img width=«565» height=«184» src=«ref-1_1289242781-8966.coolpic» v:shapes="_x0000_i1032">
решение которой будем искать в виде
<img width=«557» height=«37» src=«ref-1_1289251747-1756.coolpic» v:shapes="_x0000_i1033">
где Qi,Pi,i=1,n — прогоночные коэффициенты, подлежащие определению. Для их определения выразим из первого уравнения СЛАУ (1.1) x1 через x2, получим:
<img width=«539» height=«62» src=«ref-1_1289253503-2563.coolpic» v:shapes="_x0000_i1034">
откуда
<img width=«214» height=«64» src=«ref-1_1289256066-13706.coolpic» v:shapes="_x0000_i1035">

Из второго уравнения СЛАУ (1.1) с помощью (1.3) выразим x2 через x3, получим:
<img width=«431» height=«69» src=«ref-1_1289269772-4109.coolpic» v:shapes="_x0000_i1036">
откуда
<img width=«391» height=«79» src=«ref-1_1289273881-3851.coolpic» v:shapes="_x0000_i1037">
Продолжая этот процесс, получим из i-го уравнения СЛАУ (1.1):
<img width=«348» height=«71» src=«ref-1_1289277732-3317.coolpic» v:shapes="_x0000_i1038">

<img width=«461» height=«88» src=«ref-1_1289281049-4084.coolpic» v:shapes="_x0000_i1039">
следовательно
<img width=«333» height=«73» src=«ref-1_1289285133-3342.coolpic» v:shapes="_x0000_i1040">
Из последнего уравнения СЛАУ имеем
<img width=«490» height=«82» src=«ref-1_1289288475-4808.coolpic» v:shapes="_x0000_i1041">
то есть

Метод простых итераций

При большом числе уравнений прямые методы решения СЛАУ (за исключением метода прогонки) становятся труднореализуемыми на ЭВМ прежде всего из-за сложности хранения и обработки матриц большой размерности. В то же время характерной особенностью ряда часто встречающихся в прикладных задачах СЛАУ является разреженность матриц. Число ненулевых элементов таких матриц мало по сравнению с их размерностью. Для решения СЛАУ с разреженными матрицами предпочтительнее использовать итерационные методы.

Методы последовательных приближений, в которых при вычислении последующего приближения решения используются предыдущие, уже известные приближенные решения, называются итерационными.
Метод Зейделя решения СЛАУ

Метод простых итераций довольно медленно сходится. Для его ускорения существует метод Зейделя, заключающийся в том, что при вычислении компонента xik+1 вектора неизвестных на (k+1)-й итерации используются x1k+1, x2k+1, .., xik+1 уже вычисленные на (k+1)-й итерации. Значения остальных компонент берутся из предыдущей итерации. Так же как и в методе простых итераций строится эквивалентная СЛАУ и за начальное приближение принимается вектор правых частей. Тогда метод Зейделя для известного вектора на k-ой итерации имеет вид:

предыдущей итерации. Так же как и в методе простых итераций строится эквивалентная СЛАУ и за начальное приближение принимается вектор правых частей. Тогда метод Зейделя для известного вектора на k-ой итерации имеет вид:
<img width=«491» height=«188» src=«ref-1_1289293283-12889.coolpic» v:shapes="_x0000_i1042">
Из этой системы видно, что <img width=«188» height=«29» src=«ref-1_1289306172-12937.coolpic» v:shapes="_x0000_i1043">, где В — нижняя треугольная матрица с диагональными элементами, равными нулю, а С — верхняя треугольная матрица с диагональными элементами, отличными от нуля, α=B+C. Следовательно

При таком способе приведения исходной СЛАУ к эквивалентному виду метод простых итераций носит название метода Якоби.
<img width=«179» height=«33» src=«ref-1_1289319109-1234.coolpic» v:shapes="_x0000_i1044">
откуда
<img width=«225» height=«33» src=«ref-1_1289320343-1802.coolpic» v:shapes="_x0000_i1045">
Таким образом, метод Зейделя является методом простых итераций с матрицей правых частей α=(E-B)-1C и вектором правых частей (E-B)-1β.

Разрешим систему относительно неизвестных при ненулевых диагональных элементах, aii≠0, i= 1,n (если какой-либо коэффициент на главной диагонали равен нулю, достаточно соответствующее уравнение поменять местами с любым другим уравнением). Получим следующие выражения для компонентов вектора β и матрицы α эквивалентной системы
<img width=«538» height=«56» src=«ref-1_1289322145-3717.coolpic» v:shapes="_x0000_i1046">
В качестве нулевого приближения x(0) вектора неизвестных примем вектор правых частей x(0)=β. Тогда метод простых итераций примет вид:
<img width=«544» height=«193» src=«ref-1_1289325862-5548.coolpic» v:shapes="_x0000_i1047">
Из (1.19) видно преимущество итерационных методов по сравнению, например, с рассмотренным выше методом Гаусса. В вычислительном процессе участвуют только произведения матрицы на вектор, что позволяет работать только с ненулевыми элементами матрицы, значительно упрощая процесс хранения и обработки матриц.

Имеет место следующее достаточное условие сходимости метода простых итераций.

Метод простых итераций (1.19) сходится к единственному решению СЛАУ при любом начальном приближении x(0), если какая-либо норма матрицы α эквивалентной системы меньше единицы <img width=«35» height=«23» src=«ref-1_1289331410-443.coolpic» v:shapes="_x0000_i1048">

Если используется метод Якоби (выражения (1.18) для эквивалентной СЛАУ), то

достаточным условием сходимости является диагональное преобладание матрицы A, т.е.

<img width=«119» height=«47» src=«ref-1_1289331853-1256.coolpic» v:shapes="_x0000_i1049">(для каждой строки матрицы A модули элементов, стоящих на главной диагонали, больше суммы модулей недиагональных элементов). Очевидно, что в этом случае ||α||cменьше единицы и, следовательно, итерационный процесс (1.19) сходится.

Приведем также необходимое и достаточное условие сходимости метода простых итераций. Для сходимости итерационного процесса (1.19) необходимо и достаточно, чтобы спектр матрицы α эквивалентной системы лежал внутри круга с радиусом, равным единице.

При выполнении достаточного условия сходимости оценка погрешности решения на k- ой итерации дается выражением:
<img width=«468» height=«54» src=«ref-1_1289333109-2834.coolpic» v:shapes="_x0000_i1050">

<img width=«295» height=«76» src=«ref-1_1289335943-16970.coolpic» v:shapes="_x0000_i1051">

Следует подчеркнуть, что это неравенство дает завышенное число итераций k, поэтому редко используется на практике
    продолжение
--PAGE_BREAK--1.4 Численный метод решения задачи


Разрешим систему относительно неизвестных при ненулевых диагональных элементах, aii≠0, i= 1,n (если какой-либо коэффициент на главной диагонали равен нулю, достаточно соответствующее уравнение поменять местами с любым другим уравнением). Получим следующие выражения для компонентов вектора β и матрицы α эквивалентной системы:


<img width=«547» height=«56» src=«ref-1_1289352913-3782.coolpic» alt=«Untitled-1 копи» v:shapes=«Рисунок_x0020_3»>
При таком способе приведения исходной СЛАУ к эквивалентному виду метод простых итераций носит название метода Якоби.

В качестве нулевого приближения x(0) вектора неизвестных примем вектор правых частей x(0)=β. Тогда метод простых итераций примет вид:
<img width=«544» height=«193» src=«ref-1_1289325862-5548.coolpic» alt=«Untitled-1 копи» v:shapes=«Рисунок_x0020_4»>
Из (1.19) видно преимущество итерационных методов по сравнению, например, с рассмотренным выше методом Гаусса. В вычислительном процессе участвуют только произведения матрицы на вектор, что позволяет работать только с ненулевыми элементами матрицы, значительно упрощая процесс хранения и обработки матриц.

Имеет место следующее достаточное условие сходимости метода простых итераций.

Метод простых итераций (1.19) сходится к единственному решению СЛАУ при любом начальном приближении x(0), если какая-либо норма матрицы α эквивалентной системы меньше единицы <img width=«35» height=«24» src=«ref-1_1289362243-435.coolpic» alt=«Untitled-1 копи» v:shapes=«Рисунок_x0020_5»>

Если используется метод Якоби (выражения (1.18) для эквивалентной СЛАУ), то

достаточным условием сходимости является диагональное преобладание матрицы A, т.е.

<img width=«119» height=«48» src=«ref-1_1289362678-1289.coolpic» alt=«Untitled-1 копи» v:shapes=«Рисунок_x0020_6»>(для каждой строки матрицы A модули элементов, стоящих на главной диагонали, больше суммы модулей недиагональных элементов). Очевидно, что в этом случае ||α||cменьше единицы и, следовательно, итерационный процесс (1.19) сходится.

Приведем также необходимое и достаточное условие сходимости метода простых итераций. Для сходимости итерационного процесса (1.19) необходимо и достаточно, чтобы спектр матрицы α эквивалентной системы лежал внутри круга с радиусом, равным единице.

При выполнении достаточного условия сходимости оценка погрешности решения на k- ой итерации дается выражением:
<img width=«468» height=«54» src=«ref-1_1289333109-2834.coolpic» alt=«Untitled-1 копи» v:shapes=«Рисунок_x0020_7»>
где x· — точное решение СЛАУ.

Процесс итераций останавливается при выполнении условия, где εε≤)(kε — задаваемая вычислителем точность.

Принимая во внимание, что из (1.20) следует неравенство <img width=«154» height=«43» src=«ref-1_1289366801-1579.coolpic» alt=«Untitled-1 копи» v:shapes=«Рисунок_x0020_8»>, можно получить априорную оценку необходимого для достижения заданной точности числа итераций. При использовании в качестве начального приближения вектора β такая оценка определится неравенством:

<img width=«110» height=«53» src=«ref-1_1289368380-1493.coolpic» alt=«Untitled-1 копи» v:shapes=«Рисунок_x0020_9»>

откуда получаем априорную оценку числа итераций k при ||α||<1
<img width=«295» height=«76» src=«ref-1_1289335943-16970.coolpic» alt=«Untitled-1 копи» v:shapes=«Рисунок_x0020_10»>

Следует подчеркнуть, что это неравенство дает завышенное число итераций k, поэтому редко используется на практике.
1.6 Текст программы
program Yakobi;

uses crt;

const

 maxn=100;

type

 matrix=array[1..maxn,1..maxn] of real;

 vector=array[1..maxn] of real;

 vector1=array[1..maxn] of real;

var

 i,j,n,k1: integer;

 e,norma:real;

 a: matrix;

 b: vector;

 x2,x3: vector1;

 imya,dannble_i_rezultat,ekran:string;

procedure input(var kolvo:integer; var pogreshnostb:real; var matr1:matrix; var matr2:vector);

 {Ввод исходных данных}

 var i,j,code:integer;

 a_s: string;

 b_s: string;

 begin

 writeln('введите количество уравнений');

 readln(kolvo);

 writeln('введите погрешность вычисления');

 readln(pogreshnostb);

 writeln('введите матрицу коэффициентов при неизвестных');

 for i:=1 to kolvo do

 for j:=1 to kolvo do

 begin

 repeat

begin

 writeln('введите элемент [',i,',',j,'] и нажмите Enter');

 readln(a_s);

 if (i=j) and (a_s='0') then

 repeat

 begin

 writeln('Элементы на главной диагонали должны быть отличны от нуля. Повторите ввод');

 readln(a_s);

 end;

 until a_s<>'0';

 val(a_s,matr1[i,j],code);

 end;

 until code=0;

 end;

 writeln('введите вектор свободных коэффициентов');

 for j:=1 to kolvo do

 begin

 repeat

 writeln('введите элемент ',j,' и нажмите Enter');

 readln(b_s);

 val(b_s,matr2[j],code);

 untilcode=0;

 end;

 end;

<img width=«351» height=«119» src=«ref-1_1289386843-8535.coolpic» alt=«L_n(x)=\frac{(x-x_1)(x-x_2)(x-x_3) \ldots \ldots (x-x_n)}{(x_0-x_1)(x_0-x_2)(x_0-x_3) \ldots (x_0-x_n)} \cdot y_0 +\\ + \frac{(x-x_0)(x-x_2)(x-x_3) \ldots \ldots (x-x_n)}{(x_1-x_0)(x_1-x_2)(x_1-x_3) \ldots (x_1-x_n)} \cdot y_1 +\\ + \frac{(x-x_0)(x-x_1)(x-x_3) \ldots \ldots (x-x_n)}{(x_2-x_0)(x_2-x_1)(x_2-x_3) \ldots (x_2-x_n)} \cdot y_2 + \ldots\\ + \frac{(x-x_0)(x-x_1)(x-x_1) \ldots \ldots (x-x_n-1)}{(x_n-x_0)(x_n-x_1)(x_n-x_1) \ldots (x_n-x_n-1)} \cdot y_n.» v:shapes="_x0000_i1060">
    продолжение
--PAGE_BREAK--1.7 Тестовый пример


<img width=«619» height=«321» src=«ref-1_1289395378-21497.coolpic» v:shapes="_x0000_i1061">

2. Полиномиальная интерполяция функции методом Ньютона с разделенными разностями
2.1 Постановка задачи


Разработать схему алгоритма и написать программу на языке TurboPascal7.0 для интерполирования функции, заданной в узлах, используя метод Ньютона с разделенными разностями.
2.2 Математическая формулировка задачи


Дана табличная функция:





или

<img width=«157» height=«21» src=«ref-1_1289416875-1573.coolpic» alt=«y_i = f(x_i), i=\overline{0,n}.» v:shapes=«Рисунок_x0020_19»>Точки с координатами (xi, yi) называются узловыми точками или узлами.

Количество узлов в табличной функции равно N=n+1.

Необходимо найти значение этой функции в промежуточной точке, например, x=D, причем <img width=«97» height=«21» src=«ref-1_1289418448-1319.coolpic» alt=«D \in[x_0,x_n]» v:shapes=«Рисунок_x0020_20»>.

Для решения задачи строим интерполяционный многочлен.
2.3 Обзор существующих численных методов решения задачи


Интерполяция по Лагранжу

Интерполяционный многочлен может быть построен при помощи специальных интерполяционных формул Лагранжа, Ньютона, Стерлинга, Бесселя и др.

Интерполяционный многочлен по формуле Лагранжа имеет вид:

Докажем, что многочлен Лагранжа является интерполяционным многочленом, проходящим через все узловые точки, т.е. в узлах интерполирования xi выполняется условие Ln(xi) = yi. Для этого будем последовательно подставлять значения координат узловых точек таблицы в многочлен (2.1). В результате получим:
если x=x0, то Ln(x0) = y0,

если x=x1, то Ln(x1) = y1,

……………

если x=xn, то Ln(xn) = yn.
Это достигнуто за счет того, что в числителе каждой дроби при соответствующем значении уj, j=0,1,2,…,n отсутствует сомножитель (x-xi), в котором i=j, а знаменатель каждой дроби получен заменой переменной х на соответствующее значение хj.

Таким образом, интерполяционный многочлен Лагранжа приближает заданную табличную функцию, т.е. Ln(xi) = yi и мы можем использовать его в качестве вспомогательной функции для решения задач интерполирования, т.е. <img width=«101» height=«22» src=«ref-1_1289419767-456.coolpic» alt=«L_n(x_k) \approx y_k» v:shapes="_x0000_i1064">.

Чем больше узлов интерполирования на отрезке [x0,xn], тем точнее интерполяционный многочлен приближает заданную табличную функцию, т.е. тем точнее равенство:


<img width=«132» height=«22» src=«ref-1_1289420223-553.coolpic» alt=«f(x_k) \approx L_n(x_k).» v:shapes="_x0000_i1065">
Однако с увеличением числа узлов интерполирования возрастает степень интерполяционного многочлена n и в результате значительно возрастает объем вычислительной работы. Поэтому при большом числе узлов необходимо применять ЭВМ. В этом случае удобно находить значения функции в промежуточных точках, не получая многочлен в явном виде.

При решении задачи экстраполирования функции с помощью интерполяционного многочлена вычисление значения функции за пределами отрезка [x0,xn] обычно производят не далее, чем на один шаг h, равный наименьшей величине
<img width=«84» height=«20» src=«ref-1_1289420776-1160.coolpic» alt="\left|x_{i+1} – x_i\right|," v:shapes="_x0000_i1066">
так как за пределами отрезка [x0,xn] погрешности, как правило, увеличиваются.
Интерполяция по Ньютону

Интерполяционный многочлен по формуле Ньютона имеет вид:
<img width=«405» height=«87» src=«ref-1_1289421936-6432.coolpic» alt=«L_n(x) = f(x_0) + (x – x_0) \cdot f(x_0;x_1) +\\ + (x – x_0) \cdot(x – x_1) \cdot f(x_0;x_1;x_2) +\\ + (x – x_0) \cdot(x – x_1) \cdot(x – x_2) \cdot f(x_0;x_1;x_2;x_3) + \ldots +\\ + (x – x_0) \cdot(x – x_1) \cdot \ldots \cdot (x – x_{n-1}) \cdot f(x_0;x_1; \ldots;x_n),» v:shapes="_x0000_i1067">(2.2)
где n – степень многочлена,
<img width=«366» height=«22» src=«ref-1_1289428368-2670.coolpic» alt=«f(x_0), f(x_0;x_1), f(x_0;x_1;x_2), f(x_0;x_1;\ldots; x_n)» v:shapes="_x0000_i1068">- разделенные разности 0-го, 1-го, 2-го,…., n-го порядка, соответственно.


Сплайн-интерполяция

Сплайны стали широко использоваться в вычислительной математике сравнительно недавно. В машиностроительном черчении они применяются уже давно, так как сплайны – это лекала или гибкие линейки, деформация которых позволяет провести кривую через заданные точки (xi, уi).

Используя теорию изгиба бруса при малых деформациях, можно показать, что сплайн – это группа кубических многочленов, в местах сопряжения которых первая и вторая производные непрерывны. Такие функции называются кубическими сплайнами. Для их построения необходимо задать коэффициенты, которые единственным образом определяют многочлен в промежутке между данными точками.

Например, для некоторых функций (рис.) необходимо задать все кубические функции q1(x), q2(x), …qn(x).

В наиболее общем случае эти многочлены имеют вид:
<img border=«0» width=«342» height=«25» src=«ref-1_1289431038-2360.coolpic» alt=«q_i(x) =k_{1j} + k_{2i} x + k_{3i} x^2 + k_{4i} x^3, i=\overline{1,n},» v:shapes="_x0000_i1069">
где kij — коэффициенты, определяемые описанными ранее условиями, количество которых равно 4n. Для определения коэффициентов kij необходимо построить и решить систему порядка 4n.

Первые 2nусловий требуют, чтобы сплайны соприкасались в заданных точках:
<img border=«0» width=«208» height=«43» src=«ref-1_1289433398-2359.coolpic» alt=«q_i(x_i) = y_i, i=\overline{1,n};\\q_{i+1}(x_i) = y_i, i=\overline{0,n-1}.» v:shapes="_x0000_i1070">
Следующие (2n-2) условий требуют, чтобы в местах соприкосновения сплайнов были равны первые и вторые производные:
<img border=«0» width=«241» height=«44» src=«ref-1_1289435757-3063.coolpic» alt=«q_{i+1}'(x_i) = q'_i(x_i),& i=\overline{1,n-1};\\q''_{i+1}(x_i) = q_i''(x_i),& i=\overline{0,n-1}.» v:shapes="_x0000_i1071">

Система алгебраических уравнений имеет решение, если число уравнений соответствует числу неизвестных. Для этого необходимо ввести еще два уравнения. Обычно используются следующие условия:
<img border=«0» width=«130» height=«32» src=«ref-1_1289438820-1470.coolpic» alt=«q_1^''(x_0) = 0; q_n^''(x_n) = 0.» v:shapes="_x0000_i1072">
При построении алгоритма метода первые и вторые производные удобно аппроксимировать разделенными разностями соответствующих порядков.

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


Значения f(x0), f(x1), …, f(xn), т.е. значения табличной функции в узлах, называются разделенными разностями нулевого порядка (k=0).
<img border=«0» width=«235» height=«290» src=«ref-1_1289440290-11422.coolpic» alt=" " v:shapes="_x0000_i1073">

Отношение <img border=«0» width=«181» height=«29» src=«ref-1_1289451712-1803.coolpic» alt=«f(x_0;x_1) = \frac{f(x_1)-f(x_0)}{x_1 – x_0}» v:shapes="_x0000_i1074">называется разделенной разностью первого порядка (k=1) на участке [x0, x1] и равно разности разделенных разностей нулевого порядка на концах участка [x0, x1], разделенной на длину этого участка.

Для произвольного участка [xi, xi+1] разделенная разность первого порядка (k=1) равна
<img border=«0» width=«243» height=«44» src=«ref-1_1289453515-2334.coolpic» alt=«f(x_i;x_{i+1}) = \frac{f(x_{i+1})-f(x_i)}{x_{i+1} – x_i}.» v:shapes="_x0000_i1075">
Отношение <img border=«0» width=«246» height=«29» src=«ref-1_1289455849-2125.coolpic» alt=«f(x_0;x_1;x_2) = \frac{f(x_1;x_2)-f(x_0;x_1)}{x_2 – x_0}» v:shapes="_x0000_i1076">называется разделенной разностью второго порядка (k=2) на участке [x0, x2] и равно разности разделенных разностей первого порядка, разделенной на длину участка [x0, x2].

Для произвольного участка [xi, xi+2] разделенная разность второго порядка (k=2) равна
<img border=«0» width=«364» height=«44» src=«ref-1_1289457974-2916.coolpic» alt=«f(x_i;x_{i+1};x_{i+2}) = \frac{f(x_{i+1};x_{i+2}) — f(x_i;x_{i+1})}{x_{i+2} – x_i}» v:shapes="_x0000_i1077">
Таким образом, разделенная разность k-го порядка на участке [xi, xi+k] может быть определена через разделенные разности (k-1)-го порядка по рекуррентной формуле:
<img border=«0» width=«502» height=«44» src=«ref-1_1289460890-3315.coolpic» alt=«f(x_i;x_{i+1};x_{i+2}; \ldots;x_{i+k}) = \frac{ f(x_{i+1};x_{i+2}; \ldots;x_{i+k}) — f(x_i;x_{i+1}; \ldots;x_{i+k-1})}{x_{i+k} – x_i}» v:shapes="_x0000_i1078">(2.3)
Где <img border=«0» width=«70» height=«22» src=«ref-1_1289464205-297.coolpic» alt=«k=\overline{1,n},» v:shapes="_x0000_i1079"> <img border=«0» width=«103» height=«23» src=«ref-1_1289464502-386.coolpic» alt=«k=\overline{0,n-k},» v:shapes="_x0000_i1080"> n – степень многочлена.

Максимальное значение k равно n. Тогда i =0 и разделенная разность n-го порядка на участке [x0,xn] равна

<img border=«0» width=«371» height=«29» src=«ref-1_1289464888-2430.coolpic» alt=«f(x_0;x_i; \ldots;x_n) = \frac{ f(x_1;x_2; \ldots;x_n) — f(x_0;x_1; \ldots;x_n-1)}{x_n – x_0}» v:shapes="_x0000_i1081">,
т.е. равна разности разделенных разностей (n-1)-го порядка, разделенной на длину участка [x0,xn].

Разделенные разности
<img border=«0» width=«347» height=«22» src=«ref-1_1289467318-2397.coolpic» alt=«f(x_0;x_1), f(x_0;x_1;x_2), \ldots, f(x_0;x_1;\ldots; x_n)» v:shapes="_x0000_i1082">
являются вполне определенными числами, поэтому выражение (2.2) действительно является алгебраическим многочленом n-й степени. При этом в многочлене (2.2) все разделенные алгебраическим многочленом n-й степени. При этом в многочлене (2.2) все разделенные разности определены для участков [x0, x0+k], <img border=«0» width=«66» height=«22» src=«ref-1_1289469715-692.coolpic» alt=«k=\overline{1,n}» v:shapes="_x0000_i1083">.

Лемма: алгебраический многочлен (2.2), построенный по формулам Ньютона, действительно является интерполяционным многочленом, т.е. значение многочлена в узловых точках равно значению табличной функции
<img border=«0» width=«293» height=«22» src=«ref-1_1289470407-1995.coolpic» alt=«L_n(x_i) = f(x-i) = y_i; i=0,1, \ldots n.» v:shapes="_x0000_i1084">
Докажем это. Пусть х=х0, тогда многочлен (2.2) равен
<img border=«0» width=«171» height=«22» src=«ref-1_1289472402-1619.coolpic» alt=«L_n(x_0) = f(x_0) = y_0.» v:shapes="_x0000_i1085">
Пусть х=х1, тогда многочлен (2.2) равен
<img border=«0» width=«500» height=«43» src=«ref-1_1289474021-3662.coolpic» alt=«L_n(x_1) = f(x_0) + (x_1 – x_0) \cdot f(x_0;x_1) = f(x_0) + (x_1 – x_0) \cdot \frac{f(x_1) – f(x_0)}{x_1 – x_0} = f(x_1) = y_1» v:shapes="_x0000_i1086">
Пусть х=х2, тогда многочлен (2.2) равен

<img border=«0» width=«542» height=«52» src=«ref-1_1289477683-6407.coolpic» alt=«L_n(x_1) = f(x_0) + (x_2 – x_0) \cdot f(x_0;x_1) + (x_2 – x_0) \cdot (x_2 – x_1) \cdot f(x_0;x_1;x_2) =\\= f(x_0) + (x_2 – x_0) \cdot \frac{f(x_1) – f(x_0)}{x_1 – x_0}+ (x_2 – x_0) \cdot (x_2 – x_1) \cdot \frac{f(x_1;x_2) – f(x_0;x_1)}{x_2 – x_0}= f(x_1) = y_1» v:shapes="_x0000_i1087">
Заметим, что решение задачи интерполяции по Ньютону имеет некоторые преимущества по сравнению с решением задачи интерполяции по Лагранжу. Каждое слагаемое интерполяционного многочлена Лагранжа зависит от всех значений табличной функции yi, i=0,1,…n. Поэтому при изменении количества узловых точек N и степени многочлена n (n=N-1) интерполяционный многочлен Лагранжа требуется строить заново. В многочлене Ньютона при изменении количества узловых точек N и степени многочлена n требуется только добавить или отбросить соответствующее число стандартных слагаемых в формуле Ньютона (2.2). Это удобно на практике и ускоряет процесс вычислений.
2.5 Схема алгоритма
На рисунке 2.1 представлена схема алгоритма решения задачи №2.

На рисунке 2.2 представлена схема алгоритма ввода исходных данных (подпрограмма-процедура Vvod).

На рисунке 2.3 представлена схема алгоритма интерполяции функции по методу Ньютона с разделенными разностями (newt)

На рисунке 2.4 представлена схема алгоритма записи данных и результата в файл (подпрограмма-процедура zapisb_v_fail).

На рисунке 2.5 представлена схема алгоритма вывода содержимого записанного файла на экран (подпрограмма-процедура outputtoscreen).
2.6
Текст

программы

program newton;

 uses crt,graph;

 const c=10;

 type matr=array[0..c,0..c] of real;

 mas=array[0..c] of real;

 var x,y,koef_polinoma:mas;

 a:matr;

 b:mas;

 d1:real;

 n:integer;

 fail,fail1,ekran:text;

 procedure Vvod(var kolvo:integer; var uzel,fun:mas);

 {Процедура осуществляет ввод данных: пользователь вводит с клавиатуры

 узлы интерполяции и значения функции в них. Также определяется количество узлов.}

 var code,i:integer; s:string;

 begin

 writeln('введите количество узлов');

 readln(kolvo);

 kolvo:=kolvo-1;

 for i:=0 to kolvo do

 begin

 repeat

 writeln('введите',i,'-йузелинтерполирования');

 readln(s);

 val(s,uzel[i],code);

 until code=0;

 repeat

 writeln('введите значение функции, соответствующее данному узлу');

 readln(s);

 val(s,fun[i],code);

 until code=0;

end;

 end;
 procedure newt(var kolvo:integer; D:real; var koef,uzel,fun:mas)

 var L,P:real;

 begin

 L:=fun[0];

 P:=1;

 for i:=1 to kolvo do

 begin

 P:=P*(D-uzel[i-1]);

 for j:=1 to kolvo-i do

 begin

 fun[j]:=(fun[j-1]-fun[j])/(uzel[j+i]-uzel[i])

 end;

          koef[i]:=fun[0];

 L:=L+P*fun[0];

 end;

 end;

procedure newt(var kolvo:integer; D:real; var koef,uzel,fun:mas) {процедураинтерполяциифункцииметодомНьютона}

 var L,P:real;

 begin

 L:=fun[0];

 P:=1;

 for i:=1 to kolvo do

 begin

 P:=P*(D-uzel[i-1]);

 for j:=1 to kolvo-i do

 begin

 fun[j]:=(fun[j-1]-fun[j])/(uzel[j+i]-uzel[i])

 end;

          koef[i]:=fun[0];

 L:=L+P*fun[0];

 end;

 end;

 procedure zapisb(koef:mas; uzel,fun:mas; kolvo:integer; var f:text);

 {В данной процедуре осуществляется запись в файл данных и результата}

 var i:integer;

 begin

 assign(f,'interpol.txt');

 rewrite(f);

 for i:=0 to kolvo do writeln(f,'x= ',uzel[i]:8:4,' f(x)=',fun[i]:8:4);

 writeln(f,'Интерполяционныйполином');

 write(f,'p(x)=',koef[0]:8:4);

 for i:=1 to kolvo do if i>1 then write (f,'+(',koef[i]:8:4,')*x^',i)

 else write (f,'+(',koef[i]:8:4,')*x');

 close(f);

 end;

procedure vblvod(var f1,f2:text);

 {Вывод содержимого записанного файла на экран}

 var s1:string;

begin

 clrscr;

 assign(f1,'interpol.txt');

 reset(f1);

 assigncrt(f2);

 rewrite(f2);

 while not eof(f1) do
assigncrt(f2);

 rewrite(f2);

 while not eof(f1) do

 begin

 Readln(f1,s1);

 writeln(f2,s1);

 end;

 close(f2);

 close(f1);

 end;

 procedure grafik(kolvo:integer; uzlbl,funktsiya:mas; c:mas);

 {Построениеграфикаполученнойфункции}

 var driver,mode,Err,a1,b1,z,i,j:integer; s:string; xt,yt:real;

 begin

 driver:=detect;

 InitGraph(driver,mode,'d:\tp7\bp\bgi');

 err:=graphresult;

 if err<>grok then writeln('Ошибкаприинициализацииграфическогорежима')

 else

 begin

 Setcolor(9);

 line(320,0,320,480);

 line(0,240,640,240);

 settextstyle(smallfont,horizdir,3);

 setcolor(10);

 outtextxy(320,245,'0');

 a1:=0;

 b1:=480;

 z:=-10;

for i:=0 to 20 do

begin

 if z<>0 then

begin

 str(z,s);

 setcolor(10);

 outtextxy(a1,245,s);

 outtextxy(325,b1,s);

 setcolor(8);

 line(0,b1,640,b1);

 line(a1,0,a1,480);

 end;

outtextxy(325,b1,s);

 setcolor(8);

 line(0,b1,640,b1);

 line(a1,0,a1,480);

 end;

 a1:=a1+32;

 b1:=b1-24;

 z:=succ(z);

 end;

 setcolor(5);

 for i:=0 to kolvo do

 begin

 xt:=uzlbl[i];

 yt:=funktsiya[i];

 putpixel(round(320+xt*32),round(240-yt*24),5);

 end;

 moveto(round(320+uzlbl[0]*32),round(240-funktsiya[0]*24));

 setcolor(11);

 for i:=0 to kolvo do

 begin

 xt:=uzlbl[i];

 yt:=0;

 for j:=0 to kolvo do yt:=yt+c[j]*vozvedenie_v_stepenb(uzlbl[i],j);

 lineto(round(320+xt*32),round(240-yt*24));

 moveto(round(320+xt*32),round(240-yt*24));

 end;

 readln;

 closegraph;

 end;

 end;

{Основная часть программы}

 BEGIN

 CLRSCR;

Writeln('Программа осуществляет интерполирование функции, заданной в узлах');

 Vvod(N,X,Y);

 writeln(‘введите значение промежуточной точки’);

 readln(d1);
    продолжение
--PAGE_BREAK--2.7 Тестовый пример


 writeln('Нажмите Enter');

 readln;

 newt(N,d1,X,Y,koef_polinoma);

 zapisb(koef_polinoma,x,y,n,fail);

 vblvod(fail,fail1);

 writeln('Нажмите Enterдля просмотра графика функции, затем еще раз для выхода из программы');

 readln;

 grafik(N,X,Y,koef_polinoma);

 END.

readln(d1);

 writeln('НажмитеEnter');

 readln;

 newt(N,d1,X,Y,koef_polinoma);

 zapisb(koef_polinoma,x,y,n,fail);

 vblvod(fail,fail1);

 writeln('Нажмите Enter для просмотра графика функции, затем еще раз для выхода из программы');

 readln;

 grafik(N,X,Y,koef_polinoma);

 END.
Дана табличная функция:

Вычислить разделенные разности 1-го, 2-го, 3-го порядков (n=3) и занести их в диагональную таблицу.

Разделенные разности первого порядка:
<img border=«0» width=«425» height=«84» src=«ref-1_1289484090-8210.coolpic» alt=«f(x_0;x_1) = \frac{f(x_1) – f(x_0)}{x_1 – x_0}= \frac{1,098613 – 0,693147}{3 — 2} = 0,405466.\\f(x_1;x_2) = \frac{f(x_2) – f(x_1)}{x_2 – x_1}= \frac{1,386295 – 1,098613}{4 — 3} = 0,287682.\\f(x_2;x_3) = \frac{f(x_3) – f(x_2)}{x_3 – x_2}= \frac{1,609438 – 1,386295}{5 — 4} = 0,223143.» v:shapes="_x0000_i1088">
Разделенные разности второго порядка:
<img border=«0» width=«511» height=«58» src=«ref-1_1289492300-5942.coolpic» alt=«f(x_0;x_1;x_2) = \frac{f(x_1;x_2) – f(x_0;x_1)}{x_2 – x_0}= \frac {0,287682 – 0,405466}{4 — 2} = -0,058892.\\f(x_1;x_2;x_3) = \frac{f(x_2;x_3) – f(x_0;x_2)}{x_3 – x_1}= \frac {0,223143 – 0,287682}{5 — 3} = — 0,0322695.» v:shapes="_x0000_i1089">

Разделенная разность третьего порядка:
<img border=«0» width=«565» height=«43» src=«ref-1_1289498242-4625.coolpic» alt=«f(x_0;x_1;x_2;x_3) = \frac{f(x_1;x_3) – f(x_0;x_2)}{x_3 – x_0}= \frac {- 0,0322695 – (- 0,058892)}{5 — 2} = 0,00887416» v:shapes="_x0000_i1090">
Интерполяционный многочлен Ньютона для заданной табличной функции имеет вид:
<img border=«0» width=«585» height=«64» src=«ref-1_1289502867-9081.coolpic» alt=«L_3(x) = f(x_0) + (x – x_0) \cdot f(x_0;x_1) + (x-x_0)(x — x_1) \cdot f(x_0;x_1;x) +\\+ (x – x_0)(x – x_1)(x – x_2) \cdot f(x_0;x_1;x_2;x_3) = 0,693147 + (x — 2) \cdot 0,405466+\\ + (x-2)(x-3) \cdot (-0,058892) + (x-2)(x-3)(x-4) \cdot 0,0887416.» v:shapes="_x0000_i1091">
График интерполяционного многочлена будет таким:

procedure zapisb(koef:mas; uzel,fun:mas; kolvo:integer; var f:text);

{В данной процедуре осуществляется запись в файл данных и результата}

var i:integer;

begin

 assign(f,'interpol.txt');

 rewrite(f);

 for i:=0 to kolvo do writeln(f,'x= ',uzel[i]:8:4,' f(x)=',fun[i]:8:4);

 writeln(f,'Интерполяционныйполином');

 write(f,'p(x)=',koef[0]:8:4);

 for i:=1 to kolvo do if i>1 then write (f,'+(',koef[i]:8:4,')*x^',i)

else write (f,'+(',koef[i]:8:4,')*x');

 close(f);

end;

procedure vblvod(var f1,f2:text);

{Вывод содержимого записанного файла на экран}

var s1:string;

begin

 clrscr;

 assign(f1,'interpol.txt');

 reset(f1);

 assigncrt(f2);

 rewrite(f2);

 while not eof(f1) do

 begin

 Readln(f1,s1);

 writeln(f2,s1);

 end;

 close(f2);

 close(f1);

end;

procedure grafik(kolvo:integer; uzlbl,funktsiya:mas; c:mas);

{Построениеграфикаполученнойфункции}

var driver,mode,Err,a1,b1,z,i,j:integer; s:string; xt,yt:real;

begin

 driver:=detect;

 InitGraph(driver,mode,'d:\tp7\bp\bgi');

 err:=graphresult;

 if err<>grok then writeln('Ошибкаприинициализацииграфическогорежима')

 else

 begin

 Setcolor(9);

 line(320,0,320,480);

 line(0,240,640,240);

 settextstyle(smallfont,horizdir,3);

 setcolor(10);

 outtextxy(320,245,'0');

 a1:=0;

 b1:=480;

 z:=-10;

for i:=0 to 20 do

 begin

 if z<>0 then
<img border=«0» width=«515» height=«321» src=«ref-1_1289511948-12427.coolpic» v:shapes="_x0000_i1092">




3. Среднеквадратическое приближение функции
    продолжение
--PAGE_BREAK--3.1 Постановка задачи


Разработать схему алгоритма и написать программу на языке TurboPascal7.0 для выполнения среднеквадратического приближения функции, заданной в узлах.
3.2 Математическая формулировка задачи


Пусть имеется множество функций <img border=«0» width=«45» height=«25» src=«ref-1_1289524375-335.coolpic» v:shapes="_x0000_i1093">, принадлежащих линейному пространству функций. Под близостью в среднем интерполируемой и интерполирующей функций будем понимать результат оценки интеграла
<img border=«0» width=«261» height=«73» src=«ref-1_1289524710-1325.coolpic» v:shapes="_x0000_i1094"> , (3.1)
где <img border=«0» width=«45» height=«25» src=«ref-1_1289526035-324.coolpic» v:shapes="_x0000_i1095">  — весовая функция.

Такое приближение называют среднеквадратичным.
3.3 Обзор существующих численных методов решения задачи


Задача среднеквадратичного приближения возникает во многих областях прикладных исследований, например, при статистической обработке данных эксперимента с использованием регрессивного анализа, при оценивании параметров моделей, в задачах фильтрации и т.п.

Когда уровень неопределенности в задании приближаемой функции f(xi), i=1..m, достаточно велик, что характерно для обработки экспериментальных данных, бессмысленно требовать выполнения условий интерполирования; кроме того, число точек задания функции f(xi) часто весьма велико. Все это делает применение интерполирования мало перспективным по причинам плохой обусловленности задачи высокой размерности и проблем сходимости процесса интерполяции

Одной из наиболее простых и, поэтому, широко используемых приближающих функций является алгебраический полином
<img border=«0» width=«108» height=«46» src=«ref-1_1289526359-385.coolpic» v:shapes="_x0000_i1096">
Метод среднеквадратичного приближения обеспечивает построение полинома Pn(x), исходя из минимизации величины

Рассмотренный метод приближения минимизирует среднеквадратичное уклонение аппроксимирующего полинома от аппроксимируемой функции, но не гарантирует от значительных локальных ошибок. Для предотвращения подобной возможности используют полиномы наилучшего равномерного приближения.
<img border=«0» width=«180» height=«45» src=«ref-1_1289526744-540.coolpic» v:shapes="_x0000_i1097">
в пространстве параметров a0 , a1 ,...,an. Существуют различные подходы к решению задачи минимизации функции D(a). Простейший из них приводит к необходимости решения нормальной системы линейных алгебраических уравнений
<img border=«0» width=«180» height=«46» src=«ref-1_1289527284-588.coolpic» v:shapes="_x0000_i1098">
Однако, уже при n > 5 матрица такой системы оказывается настолько плохо обусловленной, что полученные из (3.4) значения ajоказываются мало пригодными для вычисления Pn(x). Поэтому, при необходимости построения полиномов наилучшего среднеквадратичного приближения более высоких степеней применяют другие алгоритмы, например, метод сингулярного разложения.
3.4 Численный метод решения задачи


Можно рассмотреть две задачи:

1 — подобрать функцию <img border=«0» width=«45» height=«25» src=«ref-1_1289524375-335.coolpic» v:shapes="_x0000_i1099"> так, чтобы выполнялось неравенство
<img border=«0» width=«48» height=«21» src=«ref-1_1289528207-255.coolpic» v:shapes="_x0000_i1100"> ; (3.5)
2 — найти наилучшее приближение, т.е. такую функцию <img border=«0» width=«45» height=«25» src=«ref-1_1289528462-336.coolpic» v:shapes="_x0000_i1101"> , чтобы было справедливым соотношение
 <img border=«0» width=«219» height=«67» src=«ref-1_1289528798-1130.coolpic» v:shapes="_x0000_i1102"><img border=«0» width=«228» height=«67» src=«ref-1_1289529928-1200.coolpic» v:shapes="_x0000_i1103"> . (3.6)
Далее займемся отысканием наилучшего приближения.

Разложим функцию <img border=«0» width=«45» height=«25» src=«ref-1_1289524375-335.coolpic» v:shapes="_x0000_i1104"> по системе линейно независимых функций <img border=«0» width=«56» height=«28» src=«ref-1_1289531463-366.coolpic» v:shapes="_x0000_i1105">:
<img border=«0» width=«165» height=«63» src=«ref-1_1289531829-887.coolpic» v:shapes="_x0000_i1106"> . (3.7)
В дальнейшем для сокращения записи будем пользоваться определением скалярного произведения в пространстве функций <img border=«0» width=«25» height=«28» src=«ref-1_1289532716-189.coolpic» v:shapes="_x0000_i1107">:

<img border=«0» width=«235» height=«67» src=«ref-1_1289532905-1260.coolpic» v:shapes="_x0000_i1108"> .
Подставляя (3.7) в условие (3.6), получим
<img border=«0» width=«517» height=«63» src=«ref-1_1289534165-2288.coolpic» v:shapes="_x0000_i1109"> .
Дифференцируя это выражение по <img border=«0» width=«25» height=«28» src=«ref-1_1289536453-201.coolpic» v:shapes="_x0000_i1110"> и приравнивая производные нулю, получим
<img border=«0» width=«336» height=«63» src=«ref-1_1289536654-1327.coolpic» v:shapes="_x0000_i1111"> . (3.8)
Определитель этой системы есть определитель Грама функций <img border=«0» width=«56» height=«28» src=«ref-1_1289531463-366.coolpic» v:shapes="_x0000_i1112">. В силу их линейной независимости этот определитель не равен нулю. Следовательно, из системы (3.8) можно найти коэффициенты <img border=«0» width=«25» height=«28» src=«ref-1_1289536453-201.coolpic» v:shapes="_x0000_i1113">, определяющие функцию <img border=«0» width=«45» height=«25» src=«ref-1_1289524375-335.coolpic» v:shapes="_x0000_i1114"> согласно (3.6) и минимизирующие интеграл от погрешности <img border=«0» width=«111» height=«32» src=«ref-1_1289538883-583.coolpic» v:shapes="_x0000_i1115"> . Таким образом, наилучшее среднеквадратичное приближение существует и оно единственно.

При использовании ортонормированной системы функций <img border=«0» width=«56» height=«28» src=«ref-1_1289531463-366.coolpic» v:shapes="_x0000_i1116"> система (3.8) упрощается:
<img border=«0» width=«297» height=«67» src=«ref-1_1289539832-1427.coolpic» v:shapes="_x0000_i1117">,

т.е. <img border=«0» width=«25» height=«28» src=«ref-1_1289536453-201.coolpic» v:shapes="_x0000_i1118"> являются коэффициентами Фурье, а наилучшее приближение есть ряд Фурье, обрываемый на каком-то члене.

Доказано, что в любом линейно нормированном пространстве при линейной аппроксимации вида (3.4) наилучшее приближение существует, хотя оно может быть не единственным.

В тех случаях, когда функции <img border=«0» width=«56» height=«28» src=«ref-1_1289531463-366.coolpic» v:shapes="_x0000_i1119"> не ортогональны, при <img border=«0» width=«63» height=«19» src=«ref-1_1289541826-287.coolpic» v:shapes="_x0000_i1120"> определитель Грама уменьшается, приближаясь к нулю. Тогда система становится плохо обусловленной и ее решение дает большую погрешность. В этой ситуации обычно берут не более пяти-шести членов в сумме (3.7).

В качестве <img border=«0» width=«56» height=«28» src=«ref-1_1289531463-366.coolpic» v:shapes="_x0000_i1121"> чаще всего используют полиномы Лежандра, Чебышева, Лагерра, Эрмита, ортогональные с заданным весом.

Рассмотрим частный случай, когда необходимо найти наилучшее приближение функции, заданной таблично. Для вещественных функций, заданных на конечном множестве точек, скалярное произведение определяется формулой
<img border=«0» width=«300» height=«63» src=«ref-1_1289542479-1339.coolpic» v:shapes="_x0000_i1122"> , (3.9)
где <img border=«0» width=«23» height=«20» src=«ref-1_1289543818-204.coolpic» v:shapes="_x0000_i1123">- число заданных узлов.

Условие наилучшего среднеквадратичного приближения записывается следующим образом:
 <img border=«0» width=«245» height=«63» src=«ref-1_1289544022-1105.coolpic» v:shapes="_x0000_i1124"> . (3.10)

Полагая <img border=«0» width=«165» height=«63» src=«ref-1_1289531829-887.coolpic» v:shapes="_x0000_i1125">, где <img border=«0» width=«56» height=«20» src=«ref-1_1289546014-307.coolpic» v:shapes="_x0000_i1126">, и подставляя этот многочлен в (3.10), придем к системе (3.8), в которой скалярные произведения вычисляют согласно (3.9). Описанная процедура аппроксимации носит название метода наименьших квадратов.

Наиболее употребительный вариант метода наименьших квадратов соответствует случаю степенного вида функций <img border=«0» width=«56» height=«28» src=«ref-1_1289531463-366.coolpic» v:shapes="_x0000_i1127">, т.е. <img border=«0» width=«101» height=«33» src=«ref-1_1289546687-469.coolpic» v:shapes="_x0000_i1128">, причем <img border=«0» width=«83» height=«20» src=«ref-1_1289547156-386.coolpic» v:shapes="_x0000_i1129">.

Система уравнений (3.8) при этом принимает вид
<img border=«0» width=«207» height=«63» src=«ref-1_1289547542-949.coolpic» v:shapes="_x0000_i1130"> , <img border=«0» width=«83» height=«20» src=«ref-1_1289547156-386.coolpic» v:shapes="_x0000_i1131"> , (3.11)

где <img border=«0» width=«393» height=«63» src=«ref-1_1289548877-1473.coolpic» v:shapes="_x0000_i1132">.


3.5 Схема алгоритма
<img border=«0» width=«413» height=«813» src=«ref-1_1289550350-12147.coolpic» v:shapes="_x0000_i1133">

Рис. 3.1 Основная программа


<img border=«0» width=«372» height=«888» src=«ref-1_1289562497-32891.coolpic» v:shapes="_x0000_i1134">

Рис. 3.2 Процедура ввода данных


<img border=«0» width=«200» height=«875» src=«ref-1_1289595388-11219.coolpic» v:shapes="_x0000_i1135">

Рис 3.3 Процедура среднеквадратичного приближения


program srpribl; {$R+}

uses graph;

label 1,2,3,4;

const m=2;

type mas= array [1..21] of real;

         mas1= array [1..m] of real;

         mas2= array [1..m,1..m+1] of real;

var i,j:byte;

 y1,x1:mas;

 xx1:mas1;

 a1:mas2;

procedure vvod (x,y: mas);

begin

 x[1]:=0;y[1]:=0.290234387293458; x[2]:=0.25;y[2]:=0.201247759722173;

 x[3]:=0.5;y[3]:=0.0712686786428094;x[4]:=0.75; y[4]:=0.044294935464859;

 x[5]:=1;y[5]:=-0.0745576142333448; x[6]:=1.25;y[6]:=-0.0807833694852889;

 x[7]:=1.5;y[7]:=-0.233371530473232;x[8]:=1.75;y[8]:=-0.283957027737051;

 x[9]:=2;y[9]:=-0.308697660081089;x[10]:=2.25;y[10]:=-0.42091366359964;

 x[11]:=2.5;y[11]:=-0.516991161741316;x[12]:=2.75;y[12]:=-0.427710095947851;

 x[13]:=3;y[13]:=-0.374748698528856;x[14]:=3.25;y[14]:=-0.229905794281513;

 x[15]:=3.5;y[15]:=-0.205365492492496639;x[16]:=3.75;y[16]:=-0.129155068378896;

 x[17]:=4;y[17]:=-0.0438349825330079;x[18]:=4.25;y[18]:=0.0294586319476366;

 x[19]:=4.5;y[19]:=0.132592592108995;x[20]:=4.75;y[20]:=0.206369574274868;

 x[21]:=5;y[21]:=0.26959548862651;

end;

procedure srpribl (xx:mas1;a:mas2);

var l:real;

 s,z,k1:integer;

begin

 for s:=1 to m-1 do

 for z:=s+1 to m do

 begin

 a[z,s]:=-a[z,s]/a[s,s];

 for k1:=s+1 to m+1 do a[z,k1]:=a[z,k1]+a[z,s]*a[s,k1]

 end;

 xx[m]:=a[m,m+1]/a[m,m]; writeln(' xx[',m,']=',xx[m]:2:3);

 for i:=m-1 downto 1 do

 begin

 l:=a[i,m+1];

 for j:=i+1 to m do l:=l-xx[j]*a[i,j];

 xx[i]:=l/a[i,i]; writeln(' xx[',i,']=',xx[i]:2:3)

 end

end;
procedure Grafik (var x,y:mas;xx:mas1)

var ec,gd,gm:integer;

begin

gd:=detect;

initgraph (gd,gm,' ');

ec:=graphresult;

if ec<>grok then begin

writeln ('Oshibka v inicializacii grafika');

halt (1);

end;
setcolor(white);

line (0,420,620,420);

line (0,0,0,420);

setcolor (white);

for i:=1 to 21 do begin

circle (round (x[i]*150),round (y[i]*100),1);

end;

setcolor (yellow);

for i:=1 to m do begin

circle (round (x[i]*150),round (xx[i]*100),1);

end;

setcolor (green);

for i:=2 to m do begin

line (round (x[i-1]*150),round(xx[i-1]*100),round (x[i]*150),

round (xx[i]*100));

end;

end;

begin

 vvod(x1,y1);

 for i:=1 to 2 do

 for j:=1 to 3 do a[i,j]:=0;

 a[1,1]:=21;

 for i:=1 to 21 do

 begin

 a[1,2]:=a[1,2]+x[i];

 a[2,1]:=a[2,1]+x[i];

 a[2,2]:=a[2,2]+x[i]*x[i];

 a[1,3]:=a[1,3]+y[i];

 a[2,3]:=a[2,3]+y[i]*x[i]

 end;

 srpribl(xx1,a1);

 for i:=1 to 21 do

 writeln(y[i]:2:3,' ',xx[1]+xx[2]*x[i]:2:3);

Grafik(x1,y1,xx1);

end.
3.7 Тестовый пример

<img border=«0» width=«535» height=«321» src=«ref-1_1289606607-17500.coolpic» v:shapes="_x0000_i1136">
Найти тригонометрический многочлен наилучшего среднеквадратичного приближения наименьшей степени со среднеквадратичным отклонением меньшим <img border=«0» width=«32» height=«17» src=«ref-1_1289624107-128.coolpic» alt="[Graphics:1.gif]" v:shapes="_x0000_i1137">для функции <img border=«0» width=«34» height=«17» src=«ref-1_1289624235-119.coolpic» alt="[Graphics:2.gif]" v:shapes="_x0000_i1138">

Введем функцию <img border=«0» width=«11» height=«17» src=«ref-1_1289624354-92.coolpic» alt="[Graphics:3.gif]" v:shapes="_x0000_i1139">
<img border=«0» width=«606» height=«34» src=«ref-1_1289624446-624.coolpic» hspace=«52» alt="[Graphics:4.gif]" v:shapes="_x0000_i1140">

Вычислим коэффициенты Фурье
<img border=«0» width=«564» height=«116» src=«ref-1_1289625070-4161.coolpic» hspace=«52» alt="[Graphics:5.gif]" v:shapes="_x0000_i1141">
Вычислим частичные суммы ряда Фурье
<img border=«0» width=«497» height=«53» src=«ref-1_1289629231-3197.coolpic» hspace=«52» alt="[Graphics:6.gif]" v:shapes="_x0000_i1142">
Вычислим среднеквадратичноеотклонение
<img border=«0» width=«503» height=«37» src=«ref-1_1289632428-3048.coolpic» hspace=«52» alt="[Graphics:7.gif]" v:shapes="_x0000_i1143">
Найдем минимальное <img border=«0» width=«11» height=«17» src=«ref-1_1289635476-96.coolpic» alt="[Graphics:8.gif]" v:shapes="_x0000_i1144">, при котором <img border=«0» width=«32» height=«17» src=«ref-1_1289635572-129.coolpic» alt="[Graphics:9.gif]" v:shapes="_x0000_i1145">будет меньше <img border=«0» width=«32» height=«17» src=«ref-1_1289624107-128.coolpic» alt="[Graphics:10.gif]" v:shapes="_x0000_i1146">
<img border=«0» width=«533» height=«35» src=«ref-1_1289635829-2927.coolpic» hspace=«52» alt="[Graphics:11.gif]" v:shapes="_x0000_i1147">

<img border=«0» width=«551» height=«33» src=«ref-1_1289638756-1127.coolpic» hspace=«52» alt="[Graphics:12.gif]" v:shapes="_x0000_i1148">
Следовательно многочлен степени <img border=«0» width=«18» height=«17» src=«ref-1_1289639883-104.coolpic» alt="[Graphics:13.gif]" v:shapes="_x0000_i1149">является наименьшим многочленом, удовлетворяющим нашим условиям. Построим график этого многочлениа и исходной функции
<img border=«0» width=«551» height=«35» src=«ref-1_1289639987-4679.coolpic» hspace=«52» alt="[Graphics:14.gif]" v:shapes="_x0000_i1150">

<img border=«0» width=«395» height=«244» src=«ref-1_1289644666-1977.coolpic» hspace=«62» alt="[Graphics:15.gif]" v:shapes="_x0000_i1151">
<img border=«0» width=«569» height=«33» src=«ref-1_1289646643-1613.coolpic» hspace=«52» alt="[Graphics:16.gif]" v:shapes="_x0000_i1152">
Построим график среднеквадратичного отклонения
<img border=«0» width=«531» height=«55» src=«ref-1_1289648256-5070.coolpic» hspace=«52» alt="[Graphics:17.gif]" v:shapes="_x0000_i1153">
<img border=«0» width=«389» height=«240» src=«ref-1_1289653326-1469.coolpic» hspace=«62» alt="[Graphics:18.gif]" v:shapes="_x0000_i1154">
<img border=«0» width=«569» height=«33» src=«ref-1_1289654795-1622.coolpic» hspace=«52» alt="[Graphics:19.gif]" v:shapes="_x0000_i1155">

<img border=«0» width=«524» height=«32» src=«ref-1_1289656417-1844.coolpic» hspace=«52» alt="[Graphics:20.gif]" v:shapes="_x0000_i1156">
Найдем минимальное <img border=«0» width=«9» height=«14» src=«ref-1_1289658261-265.coolpic» alt="[Graphics:8.gif]" v:shapes="_x0000_i1157">, при котором <img border=«0» width=«32» height=«14» src=«ref-1_1289658526-459.coolpic» alt="[Graphics:9.gif]" v:shapes="_x0000_i1158">будет меньше <img border=«0» width=«32» height=«14» src=«ref-1_1289658985-444.coolpic» alt="[Graphics:10.gif]" v:shapes="_x0000_i1159">

<img border=«0» width=«541» height=«35» src=«ref-1_1289659429-2413.coolpic» hspace=«52» alt="[Graphics:11.gif]" v:shapes="_x0000_i1160">

<img border=«0» width=«550» height=«32» src=«ref-1_1289661842-826.coolpic» hspace=«52» alt="[Graphics:12.gif]" v:shapes="_x0000_i1161">
Следовательно, многочлен степени <img border=«0» width=«17» height=«14» src=«ref-1_1289662668-332.coolpic» alt="[Graphics:13.gif]" v:shapes="_x0000_i1162">является наименьшим многочленом, удовлетворяющим нашим условиям. Построим график этого многочлениа и исходной функции
<img border=«0» width=«550» height=«35» src=«ref-1_1289663000-3865.coolpic» hspace=«52» alt="[Graphics:14.gif]" v:shapes="_x0000_i1163">
<img border=«0» width=«395» height=«245» src=«ref-1_1289644666-1977.coolpic» hspace=«62» alt="[Graphics:15.gif]" v:shapes="_x0000_i1164">
<img border=«0» width=«611» height=«32» src=«ref-1_1289668842-585.coolpic» hspace=«52» alt="[Graphics:16.gif]" v:shapes="_x0000_i1165">
Построим график среднеквадратичного отклонения
<img border=«0» width=«530» height=«55» src=«ref-1_1289669427-4270.coolpic» hspace=«52» alt="[Graphics:17.gif]" v:shapes="_x0000_i1166">


<img border=«0» width=«389» height=«239» src=«ref-1_1289653326-1469.coolpic» hspace=«62» alt="[Graphics:18.gif]" v:shapes="_x0000_i1167">
<img border=«0» width=«547» height=«32» src=«ref-1_1289675166-1418.coolpic» hspace=«52» alt="[Graphics:19.gif]" v:shapes="_x0000_i1168">
<img border=«0» width=«547» height=«32» src=«ref-1_1289676584-1854.coolpic» hspace=«52» alt="[Graphics:20.gif]" v:shapes="_x0000_i1169">


4. Численное интегрирование функций методом Гаусса

    продолжение
--PAGE_BREAK--4.1 Постановка задачи


Разработать схему алгоритма и написать программу на языке TurboPascal7.0 для вычисления интегралов функции, используя метод Гаусса.
4.2 Математическая формулировка задачи


Пусть функция задана на стандартном интервале <img border=«0» width=«44» height=«21» src=«ref-1_1289678438-220.coolpic» v:shapes="_x0000_i1170">. Задача состоит в том, чтобы подобрать точки <img border=«0» width=«99» height=«28» src=«ref-1_1289678658-346.coolpic» v:shapes="_x0000_i1171"> и коэффициенты <img border=«0» width=«116» height=«28» src=«ref-1_1289679004-382.coolpic» v:shapes="_x0000_i1172"> так, чтобы квадратурная формула
<img border=«0» width=«188» height=«67» src=«ref-1_1289679386-1017.coolpic» v:shapes="_x0000_i1173"> 
была точной для всех полиномов наивысшей возможной степени.
4.3 Обзор существующих численных методов решения задачи


Одномерный случай

Основная идея большинства методов численного интегрирования состоит в замене подынтегральной функции на более простую, интеграл от которой легко вычисляется аналитически. При этом для оценки значения интеграла получаются формулы вида
<img border=«0» width=«135» height=«48» src=«ref-1_1289680403-972.coolpic» alt=«I \approx \sum_{i=1}^{n} w_i\, f(x_i),» v:shapes="_x0000_i1174">

где <img border=«0» width=«17» height=«12» src=«ref-1_1289681375-1014.coolpic» alt=«n\,\!» v:shapes="_x0000_i1175">— число точек, в которых вычисляется значение подынтегральной функции. Точки <img border=«0» width=«16» height=«12» src=«ref-1_1289682389-238.coolpic» alt=«x_i\,\!» v:shapes="_x0000_i1176"> называются узлами метода, числа <img border=«0» width=«19» height=«13» src=«ref-1_1289682627-274.coolpic» alt=«w_i\,\!» v:shapes="_x0000_i1177">— весами узлов. При замене подынтегральной функции на полином нулевой, первой и второй степени получаются соответственно методы прямоугольников, трапеций и парабол (Симпсона). Часто формулы для оценки значения интеграла называют квадратурными формулами.
Метод прямоугольников

Метод прямоугольниковполучается при замене подынтегральной функции на константу. В качестве константы можно взять значение функции в любой точке отрезка <img border=«0» width=«34» height=«20» src=«ref-1_1289682901-356.coolpic» alt="\left[a, b\right]\,\!" v:shapes="_x0000_i1178">. Наиболее часто используются значения функции в середине отрезка и на его концах. Соответствующие модификации носят названия методов средних прямоугольников, левых прямоугольников и правых прямоугольников. Формула для приближенного вычисления значения определённого интеграла методом прямоугольников имеет вид
<img border=«0» width=«131» height=«20» src=«ref-1_1289683257-755.coolpic» alt=«I \approx f(x) (b-a)» v:shapes="_x0000_i1179">,
где <img border=«0» width=«88» height=«43» src=«ref-1_1289684012-1317.coolpic» alt=«x=\frac{a+b}{2}» v:shapes="_x0000_i1180">, <img border=«0» width=«10» height=«9» src=«ref-1_1289685329-210.coolpic» alt=«a\ » v:shapes="_x0000_i1181"> или <img border=«0» width=«8» height=«15» src=«ref-1_1289685539-229.coolpic» alt=«b\ » v:shapes="_x0000_i1182">, соответственно.
Метод трапеций

Если функцию на каждом из частичных отрезков аппроксимировать прямой, проходящей через конечные значения, то получим метод трапеций.

Площадь трапеции на каждом отрезке:
<img border=«0» width=«264» height=«43» src=«ref-1_1289685768-1228.coolpic» alt="~I_i \approx \frac{f(x_{i-1})+f(x_{i})}{2} (x_{i}-x_{i-1})" v:shapes="_x0000_i1183">


Погрешность аппроксимации на каждом отрезке:
<img border=«0» width=«161» height=«46» src=«ref-1_1289686996-1120.coolpic» alt="~\left| R_{i} \right| \leqslant \frac{\left( b-a \right)^3}{12n^2} M_{2,i}" v:shapes="_x0000_i1184">, где <img border=«0» width=«193» height=«34» src=«ref-1_1289688116-1111.coolpic» alt=«M_{2,i}=\max_{x\mathcal{2}[x_{i-1},x_i]} \left| f''(x) \right| » v:shapes="_x0000_i1185">
Полная формула трапеций в случае деления всего промежутка интегрирования на отрезки одинаковой длины h:
<img border=«0» width=«293» height=«51» src=«ref-1_1289689227-1794.coolpic» alt="~I \approx h\left( \frac{f(x_{0})+f(x_{n})}{2} + \sum_{i=1}^{n-1}f(x_{i})\right)" v:shapes="_x0000_i1186">,
где <img border=«0» width=«82» height=«41» src=«ref-1_1289691021-495.coolpic» alt=«h=\frac{b-a}{n}» v:shapes="_x0000_i1187">

Погрешность формулы трапеций:
<img border=«0» width=«146» height=«46» src=«ref-1_1289691516-1050.coolpic» alt="~\left| R \right| \leqslant \frac{\left( b-a \right)^3}{12n^2} M_{2}" v:shapes="_x0000_i1188">, где <img border=«0» width=«156» height=«34» src=«ref-1_1289692566-1036.coolpic» alt=«M_{2}=\max_{x\mathcal{2}[a,b]} \left| f''(x) \right| » v:shapes="_x0000_i1189">
Метод парабол (метод Симпсона)

Использовав три точки отрезка интегрирования можно заменить подынтегральную функцию параболой. Обычно в качестве таких точек используют концы отрезка и его среднюю точку. В этом случае формула имеет очень простой вид
<img border=«0» width=«329» height=«51» src=«ref-1_1289693602-1930.coolpic» alt=«I \approx \frac{b-a}{6}\left(f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right)» v:shapes="_x0000_i1190">.

Если разбить интервал интегрирования на 2N равных частей, то имеем
<img border=«0» width=«521» height=«42» src=«ref-1_1289695532-4503.coolpic» alt=«I \approx \frac{b-a}{6N}\left(f_0 + 4 \left(f_1 + f_3 +… +f_{2N-1}\right) + 2 \left(f_2 + f_4 +… +f_{2N-2}\right) + f_{2N}\right)» v:shapes="_x0000_i1191">,

где <img border=«0» width=«181» height=«51» src=«ref-1_1289700035-1206.coolpic» alt=«f_i = f\left(a+\frac{(b-a)i}{2N}\right)» v:shapes="_x0000_i1192">.
Если разбить интервал интегрирования на 2N равных частей, то имеем
<img border=«0» width=«546» height=«42» src=«ref-1_1289701241-4584.coolpic» alt=«I \approx \frac{b-a}{6N}\left(f_0 + 4 \left(f_1 + f_3 +… +f_{2N-1}\right) + 2 \left(f_2 + f_4 +… +f_{2N-2}\right) + f_{2N}\right)» v:shapes="_x0000_i1193">,

где <img border=«0» width=«181» height=«51» src=«ref-1_1289700035-1206.coolpic» alt=«f_i = f\left(a+\frac{(b-a)i}{2N}\right)» v:shapes="_x0000_i1194">.
Увеличение точности

Приближение функции одним полиномом на всем отрезке интегрирования, как правило, приводит к большой ошибке в оценке значения интеграла.

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

При стремлении количества разбиений к бесконечности, оценка интеграла стремится к его истинному значению для любого численного метода.

Приведённые выше методы допускают простую процедуру уменьшения шага в два раза, при этом на каждом шаге требуется вычислять значения функции только во вновь добавленных узлах. Для оценки погрешности вычислений используется правило Рунге.
Метод Гаусса

Описанные выше методы используют фиксированные точки отрезка (концы и середину) и имеют низкий порядок точности (0 — методы правых и левых прямоугольников, 1 — методы средних прямоугольников и трапеций, 3 — метод парабол (Симпсона)). Если мы можем выбирать точки, в которых мы вычисляем значения функции <img border=«0» width=«35» height=«21» src=«ref-1_1289707031-369.coolpic» alt=«f(x)\,\!» v:shapes="_x0000_i1195">, то можно при том же количестве вычислений подынтегральной функции получить методы более высокого порядка точности. Так для двух (как в методе трапеций) вычислений значений подынтегральной функции, можно получить метод уже не 1-го, а 3-го порядка точности:
<img border=«0» width=«452» height=«51» src=«ref-1_1289707400-2564.coolpic» alt=«I \approx \frac{b-a}{2}\left(f\left(\frac{a+b}{2} — \frac{b-a}{2\sqrt{3}} \right)+f\left(\frac{a+b}{2} + \frac{b-a}{2\sqrt{3}} \right) \right)» v:shapes="_x0000_i1196">.
В общем случае, используя <img border=«0» width=«17» height=«12» src=«ref-1_1289681375-1014.coolpic» alt=«n\,\!» v:shapes="_x0000_i1197">точек, можно получить метод с порядком точности <img border=«0» width=«54» height=«15» src=«ref-1_1289710978-352.coolpic» alt=«2n-1\,\!» v:shapes="_x0000_i1198">. Значения узлов метода Гаусса по <img border=«0» width=«17» height=«12» src=«ref-1_1289681375-1014.coolpic» alt=«n\,\!» v:shapes="_x0000_i1199">точкам являются корнями полинома Лежандра степени <img border=«0» width=«17» height=«12» src=«ref-1_1289681375-1014.coolpic» alt=«n\,\!» v:shapes="_x0000_i1200">.

Значения узлов метода Гаусса и их весов приводятся в справочниках специальных функций. Наиболее известен метод Гаусса по пяти точкам.

В общем случае, используя <img border=«0» width=«17» height=«12» src=«ref-1_1289681375-1014.coolpic» alt=«n\,\!» v:shapes="_x0000_i1201">точек, можно получить метод с порядком точности <img border=«0» width=«54» height=«15» src=«ref-1_1289710978-352.coolpic» alt=«2n-1\,\!» v:shapes="_x0000_i1202">. Значения узлов метода Гаусса по <img border=«0» width=«17» height=«12» src=«ref-1_1289681375-1014.coolpic» alt=«n\,\!» v:shapes="_x0000_i1203">точкам являются корнями полинома Лежандра степени <img border=«0» width=«17» height=«12» src=«ref-1_1289681375-1014.coolpic» alt=«n\,\!» v:shapes="_x0000_i1204">.

Значения узлов метода Гаусса и их весов приводятся в справочниках специальных функций. Наиболее известен метод Гаусса по пяти точкам.

Метод Гаусса-Кронрода

Недостаток метода Гаусса состоит в том, что он не имеет лёгкого (с вычислительной точки зрения) пути оценки погрешности полученного значения интеграла. Использование правила Рунге требует вычисления подынтегральной функции примерно в таком же числе точек, не давая при этом практически никакого выигрыша точности, в отличие от простых методов, где точность увеличивается в разы при каждом новом разбиении. Кронродом был предложен следующий метод оценки значения интеграла
<img border=«0» width=«239» height=«52» src=«ref-1_1289716752-1585.coolpic» alt=«I \approx \sum_{i=1}^{n} a_i\, f(x_i) + \sum_{i=1}^{n+1} b_i\, f(y_i)» v:shapes="_x0000_i1205">,

где <img border=«0» width=«16» height=«12» src=«ref-1_1289682389-238.coolpic» alt=«x_i\,\!» v:shapes="_x0000_i1206">— узлы метода Гаусса по <img border=«0» width=«17» height=«12» src=«ref-1_1289681375-1014.coolpic» alt=«n\,\!» v:shapes="_x0000_i1207">точкам, а <img border=«0» width=«54» height=«16» src=«ref-1_1289719589-395.coolpic» alt=«3n+2\,\!» v:shapes="_x0000_i1208"> параметров <img border=«0» width=«15» height=«12» src=«ref-1_1289719984-239.coolpic» alt=«a_i\,\!» v:shapes="_x0000_i1209">, <img border=«0» width=«13» height=«17» src=«ref-1_1289720223-249.coolpic» alt=«b_i\,\!» v:shapes="_x0000_i1210">, <img border=«0» width=«15» height=«13» src=«ref-1_1289720472-251.coolpic» alt=«y_i\,\!» v:shapes="_x0000_i1211">подобраны таким образом, чтобы порядок точности метода был равен <img border=«0» width=«54» height=«16» src=«ref-1_1289720723-371.coolpic» alt=«3n+1\,\!» v:shapes="_x0000_i1212">.

Тогда для оценки погрешности можно использовать эмпирическую формулу:
<img border=«0» width=«170» height=«24» src=«ref-1_1289721094-838.coolpic» alt="\Delta = \left(200 |I — I_G|\right)^{1.5}" v:shapes="_x0000_i1213">,
где <img border=«0» width=«20» height=«17» src=«ref-1_1289721932-259.coolpic» alt=«I_G\,\!» v:shapes="_x0000_i1214">— приближённое значение интеграла, полученное методом Гаусса по <img border=«0» width=«14» height=«12» src=«ref-1_1289722191-310.coolpic» alt=«n\,\!» v:shapes="_x0000_i1215">точкам.


    продолжение
--PAGE_BREAK--4.4 Численный метод решения задачи


Пусть функция задана на стандартном интервале <img border=«0» width=«51» height=«24» src=«ref-1_1289722501-232.coolpic» v:shapes="_x0000_i1216">. Задача состоит в том, чтобы подобрать точки <img border=«0» width=«99» height=«28» src=«ref-1_1289678658-346.coolpic» v:shapes="_x0000_i1217"> и коэффициенты <img border=«0» width=«116» height=«28» src=«ref-1_1289679004-382.coolpic» v:shapes="_x0000_i1218"> так, чтобы квадратурная формула
<img border=«0» width=«188» height=«67» src=«ref-1_1289679386-1017.coolpic» v:shapes="_x0000_i1219"> (4.1)
была точной для всех полиномов наивысшей возможной степени.

Ввиду того, что имеется <img border=«0» width=«27» height=«20» src=«ref-1_1289724478-228.coolpic» v:shapes="_x0000_i1220"> параметров <img border=«0» width=«24» height=«28» src=«ref-1_1289724706-181.coolpic» v:shapes="_x0000_i1221">и <img border=«0» width=«16» height=«28» src=«ref-1_1289724887-171.coolpic» v:shapes="_x0000_i1222"> <img border=«0» width=«120» height=«25» src=«ref-1_1289725058-470.coolpic» v:shapes="_x0000_i1223">, а полином степени <img border=«0» width=«53» height=«20» src=«ref-1_1289725528-270.coolpic» v:shapes="_x0000_i1224"> определяется <img border=«0» width=«27» height=«20» src=«ref-1_1289724478-228.coolpic» v:shapes="_x0000_i1225"> коэффициентами, эта наивысшая степень в общем случае <img border=«0» width=«96» height=«20» src=«ref-1_1289726026-379.coolpic» v:shapes="_x0000_i1226">.

Запишем полином в виде <img border=«0» width=«144» height=«63» src=«ref-1_1289726405-748.coolpic» v:shapes="_x0000_i1227"> и подставим в (4.1). Получим

<img border=«0» width=«256» height=«67» src=«ref-1_1289727153-1251.coolpic» v:shapes="_x0000_i1228"> ,

<img border=«0» width=«240» height=«67» src=«ref-1_1289728404-1242.coolpic» v:shapes="_x0000_i1229"> .
Приравнивая выражения при одинаковых коэффициентах <img border=«0» width=«25» height=«28» src=«ref-1_1289536453-201.coolpic» v:shapes="_x0000_i1230">получим
<img border=«0» width=«139» height=«67» src=«ref-1_1289729847-735.coolpic» v:shapes="_x0000_i1231"> , <img border=«0» width=«164» height=«25» src=«ref-1_1289730582-560.coolpic» v:shapes="_x0000_i1232"> ,

<img border=«0» width=«413» height=«80» src=«ref-1_1289731142-1727.coolpic» v:shapes="_x0000_i1233"> , <img border=«0» width=«116» height=«25» src=«ref-1_1289732869-454.coolpic» v:shapes="_x0000_i1234"> . <img border=«0» width=«12» height=«23» src=«ref-1_1289733323-73.coolpic» v:shapes="_x0000_i1235">
Итак, <img border=«0» width=«24» height=«28» src=«ref-1_1289724706-181.coolpic» v:shapes="_x0000_i1236">и <img border=«0» width=«16» height=«28» src=«ref-1_1289724887-171.coolpic» v:shapes="_x0000_i1237"> находят из системы <img border=«0» width=«27» height=«20» src=«ref-1_1289724478-228.coolpic» v:shapes="_x0000_i1238"> уравнений
<img border=«0» width=«85» height=«63» src=«ref-1_1289733976-433.coolpic» v:shapes="_x0000_i1239"> ,

<img border=«0» width=«99» height=«63» src=«ref-1_1289734409-507.coolpic» v:shapes="_x0000_i1240">,

<img border=«0» width=«101» height=«63» src=«ref-1_1289734916-556.coolpic» v:shapes="_x0000_i1241"> , (4.2)

 … .

 <img border=«0» width=«123» height=«63» src=«ref-1_1289735472-566.coolpic» v:shapes="_x0000_i1242"> .

Система (4.2) нелинейная, и ее решение найти довольно трудно. Рассмотрим еще один прием нахождения<img border=«0» width=«24» height=«28» src=«ref-1_1289724706-181.coolpic» v:shapes="_x0000_i1243">и <img border=«0» width=«16» height=«28» src=«ref-1_1289724887-171.coolpic» v:shapes="_x0000_i1244">. Свойства полиномов Лежандра
<img border=«0» width=«243» height=«61» src=«ref-1_1289736390-1121.coolpic» v:shapes="_x0000_i1245"> , <img border=«0» width=«103» height=«25» src=«ref-1_1289737511-404.coolpic» v:shapes="_x0000_i1246">
таковы:
1) <img border=«0» width=«239» height=«33» src=«ref-1_1289737915-716.coolpic» v:shapes="_x0000_i1247"> , <img border=«0» width=«103» height=«25» src=«ref-1_1289737511-404.coolpic» v:shapes="_x0000_i1248"> ;

2) <img border=«0» width=«365» height=«67» src=«ref-1_1289739035-1522.coolpic» v:shapes="_x0000_i1249"> ;
3) полином Лежандра <img border=«0» width=«53» height=«28» src=«ref-1_1289740557-338.coolpic» v:shapes="_x0000_i1250"> имеет <img border=«0» width=«15» height=«15» src=«ref-1_1289740895-171.coolpic» v:shapes="_x0000_i1251"> различных и действительных корней, расположенных на интервале <img border=«0» width=«51» height=«24» src=«ref-1_1289722501-232.coolpic» v:shapes="_x0000_i1252">.

Составим по узлам интегрирования многочлен <img border=«0» width=«15» height=«15» src=«ref-1_1289740895-171.coolpic» v:shapes="_x0000_i1253">-й степени
<img border=«0» width=«173» height=«63» src=«ref-1_1289741469-756.coolpic» v:shapes="_x0000_i1254"> .
Функция <img border=«0» width=«173» height=«28» src=«ref-1_1289742225-758.coolpic» v:shapes="_x0000_i1255"> при <img border=«0» width=«80» height=«20» src=«ref-1_1289742983-339.coolpic» v:shapes="_x0000_i1256"> есть многочлен степени не выше <img border=«0» width=«53» height=«20» src=«ref-1_1289725528-270.coolpic» v:shapes="_x0000_i1257">. Значит для этой функции формула Гаусса справедлива:
<img border=«0» width=«356» height=«67» src=«ref-1_1289743592-1632.coolpic» v:shapes="_x0000_i1258"> , (4.3)


так как <img border=«0» width=«89» height=«28» src=«ref-1_1289745224-430.coolpic» v:shapes="_x0000_i1259"> .

Разложим <img border=«0» width=«56» height=«28» src=«ref-1_1289745654-343.coolpic» v:shapes="_x0000_i1260"> в ряд по ортогональным многочленам Лежандра:
<img border=«0» width=«172» height=«63» src=«ref-1_1289745997-902.coolpic» v:shapes="_x0000_i1261"> ,

<img border=«0» width=«459» height=«67» src=«ref-1_1289746899-2050.coolpic» v:shapes="_x0000_i1262"> ,

<img border=«0» width=«80» height=«20» src=«ref-1_1289742983-339.coolpic» v:shapes="_x0000_i1263">,
т.е. все коэффициенты <img border=«0» width=«60» height=«28» src=«ref-1_1289749288-301.coolpic» v:shapes="_x0000_i1264"> при <img border=«0» width=«80» height=«20» src=«ref-1_1289742983-339.coolpic» v:shapes="_x0000_i1265"> . Значит <img border=«0» width=«56» height=«28» src=«ref-1_1289745654-343.coolpic» v:shapes="_x0000_i1266"> с точностью до численного множителя совпадает с <img border=«0» width=«53» height=«28» src=«ref-1_1289740557-338.coolpic» v:shapes="_x0000_i1267">. Таким образом, узлами формулы Гаусса являются нулимногочлена Лежандра степени <img border=«0» width=«15» height=«15» src=«ref-1_1289740895-171.coolpic» v:shapes="_x0000_i1268">.

Зная <img border=«0» width=«16» height=«28» src=«ref-1_1289724887-171.coolpic» v:shapes="_x0000_i1269">, из линейной теперь системы первых <img border=«0» width=«15» height=«15» src=«ref-1_1289740895-171.coolpic» v:shapes="_x0000_i1270"> (4.2) легко найти коэффициенты являются нули многочлена Лежандра степени <img border=«0» width=«15» height=«15» src=«ref-1_1289740895-171.coolpic» v:shapes="_x0000_i1271">.

Зная <img border=«0» width=«16» height=«28» src=«ref-1_1289724887-171.coolpic» v:shapes="_x0000_i1272">, из линейной теперь системы первых <img border=«0» width=«15» height=«15» src=«ref-1_1289740895-171.coolpic» v:shapes="_x0000_i1273"> (4.2) легко найти коэффициенты <img border=«0» width=«24» height=«28» src=«ref-1_1289724706-181.coolpic» v:shapes="_x0000_i1274"><img border=«0» width=«120» height=«25» src=«ref-1_1289725058-470.coolpic» v:shapes="_x0000_i1275">. Определитель этой системы есть определитель Вандермонда.

Формулу <img border=«0» width=«185» height=«67» src=«ref-1_1289752286-1017.coolpic» v:shapes="_x0000_i1276">, в которой <img border=«0» width=«16» height=«28» src=«ref-1_1289724887-171.coolpic» v:shapes="_x0000_i1277">- нули полинома Лежандра <img border=«0» width=«49» height=«28» src=«ref-1_1289753474-332.coolpic» v:shapes="_x0000_i1278">, а <img border=«0» width=«24» height=«28» src=«ref-1_1289753806-182.coolpic» v:shapes="_x0000_i1279"> определяют из (3.3), называют квадратурной формулой Гаусса.


4.5 Схема алгоритма

 <img border=«0» width=«98» height=«778» src=«ref-1_1289753988-7292.coolpic» v:shapes="_x0000_i1280">

Рис. 4.1 Основная программа метода


<img border=«0» width=«138» height=«778» src=«ref-1_1289761280-8934.coolpic» v:shapes="_x0000_i1281">

Рис .4.2 Интегрирование методом Гаусса


<img border=«0» width=«423» height=«401» src=«ref-1_1289770214-9349.coolpic» v:shapes="_x0000_i1282">

Рис 4.3 Процедура подсчета коэффициентов по методу Гаусса


    продолжение
--PAGE_BREAK--