Численные методы решения типовых математических задач (стр. 7 из 8)
,где
— приближённое значение интеграла, полученное методом Гаусса по 
точкам.4.4 Численный метод решения задачи
Пусть функция задана на стандартном интервале
. Задача состоит в том, чтобы подобрать точки 
и коэффициенты 
так, чтобы квадратурная формула 
(4.1)была точной для всех полиномов наивысшей возможной степени.
Ввиду того, что имеется
параметров 
и 
, а полином степени 
определяется коэффициентами, эта наивысшая степень в общем случае 
.Запишем полином в виде
и подставим в (4.1). Получим 
, 
.Приравнивая выражения при одинаковых коэффициентах
получим 
, 
, 
, 
. 
Итак,
и находят из системы уравнений 
, 
, 
, (4.2). . . . . . .
.Система (4.2) нелинейная, и ее решение найти довольно трудно. Рассмотрим еще один прием нахождения
и . Свойства полиномов Лежандра 
, 
таковы:
1)
, ;2)
;3) полином Лежандра
имеет 
различных и действительных корней, расположенных на интервале 
.Составим по узлам интегрирования многочлен
-й степени 
.Функция
при 
есть многочлен степени не выше 
. Значит для этой функции формула Гаусса справедлива: 
, (4.3)так как
.Разложим
в ряд по ортогональным многочленам Лежандра: 
, 
, 
,т.е. все коэффициенты
при . Значит с точностью до численного множителя совпадает с 
. Таким образом, узлами формулы Гаусса являются нули многочлена Лежандра степени .Зная
, из линейной теперь системы первых (4.2) легко найти коэффициенты являются нули многочлена Лежандра степени .Зная
, из линейной теперь системы первых (4.2) легко найти коэффициенты 
. Определитель этой системы есть определитель Вандермонда.Формулу
, в которой - нули полинома Лежандра 
, а 
определяют из (3.3), называют квадратурной формулой Гаусса.
Рис. 4.1 Основная программа метода
Рис .4.2 Интегрирование методом Гаусса
Рис 4.3 Процедура подсчета коэффициентов по методу Гаусса
4.6 Текст программы
uses crt,graph;
var aaa,bbb,kkk:real;
const
g10c1=0.9739065285/6.2012983932;
g10c2=0.8650633667/6.2012983932;
g10c3=0.6794095683/6.2012983932;
g10c4=0.4333953941/6.2012983932;
g10c5=0.1488743390/6.2012983932;
g10x1=0.0666713443/6.2012983932;
g10x2=0.1494513492/6.2012983932;
g10x3=0.2190863625/6.2012983932;
g10x4=0.2692667193/6.2012983932;
g10x5=0.2955242247/6.2012983932;
function F(x:real):real;{интегрируемая функция}
begin
F:= ln(1+x*x*x);
end;
function gauss_calc(a,b:real):real; {метод Гаусса}
var n,m,s,s1,s2,s3,s4,s5 :real;
begin
m:=(b+a)/2; n:=(b-a)/2;
s1:=g10c1*(f(m+n*g10x1)+f(m-n*g10x1));
s2:=g10c2*(f(m+n*g10x2)+f(m-n*g10x2));
s3:=g10c3*(f(m+n*g10x3)+f(m-n*g10x3));
s4:=g10c4*(f(m+n*g10x4)+f(m-n*g10x4));
s5:=g10c5*(f(m+n*g10x5)+f(m-n*g10x5));
s:=s1+s2+s3+s4+s5;
gauss_calc:=s*(b-a);
end;
{рекурсивная ф-ция подсчета с заданной точностью}
{ gc - ранее посчитаный интеграл на интервале (a,b)}
function gauss(a,b,eps,gc:real):real;
var t,ga,gb :real;
begin
t:=(a+b)/2; {разбиваем интервал на две половинки}
ga:=gauss_calc(a,t); {в каждой половинке считаем интеграл}
gb:=gauss_calc(t,b);
if abs(ga+gb-gc)>eps then {проверяем точность вычислений}
begin
ga:=gauss(a,t,eps/2,ga); {рекурсия для первой половинки}
gb:=gauss(t,b,eps/2,gb); {рекурсия для второй половинки}
end; {при этом точность повышаем, чтобы }
{при сложении ошибка не накапливалась}
gauss:=ga+gb; {интеграл = сумме интегралов половинок}
end;
procedure inputnum(prm:string;var num:real;lb,ub:real);
var q:boolean;
begin
repeat
write('Введите ',prm,' ');readln(num);
q:=(lb>=num) or (num>ub);
if q then writeln('Число должно быть от ', lb:0:0,' до ',ub:0:0);
until not q;
end;
function main_menu:integer;
var i:integer;
begin
Writeln('==========================================================');
Writeln('0 - выход');
Writeln('1 - решать тестовый пример a=0 b=1.2 k=10 eps = 0.0001');
Writeln('2 - решать пример с произвольными a, b, k, eps');
Writeln('----------------------------------------------------------');
Write('Выбор >>> ');readln(i);
Writeln('==========================================================');
main_menu:=i;
end;
{Вывод графика}
procedure outputgraph(a,b,a1,b1:real;n:integer);
var i,j,j1,k:integer;t,y1,y2,x1,x2,x,y:double;s:string;
begin
clearviewport;
x1:=a1-1;x2:=b1+1;
if x1<0.5 then x1:=-0.5;
y2:=f(ln(bbb/aaa)/(bbb-aaa))*1.2;
y1:=-y2;
{Линия y=0}
setcolor(15);
y:=0;
j:=trunc(480*(y-y2)/(y1-y2));
line(0,j,639,j);
{Линии x=a,x=b}
setcolor(5);
j:=trunc(480*(-y2)/(y1-y2));
i:=trunc(640*(a-x1)/(x2-x1));
j1:=trunc(480*(f(a)-y2)/(y1-y2));
line(i,j,i,j1);
i:=trunc(640*(b-x1)/(x2-x1));
j:=trunc(480*(-y2)/(y1-y2));
j1:=trunc(480*(f(b)-y2)/(y1-y2));