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

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

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

 (4.1)

была точной для всех полиномов наивысшей возможной степени.

Ввиду того, что имеется  параметров и  , а полином степени  определяется  коэффициентами, эта наивысшая степень в общем случае .

Запишем полином в виде  и подставим в (4.1). Получим

 ,

 .

Приравнивая выражения при одинаковых коэффициентах получим

 ,  ,

 ,  .

Итак, и  находят из системы  уравнений

 ,

,

 , (4.2)

 . . . . . . .

  .

Система (4.2) нелинейная, и ее решение найти довольно трудно. Рассмотрим еще один прием нахожденияи . Свойства полиномов Лежандра

 ,

таковы:

1)  ,  ;

2)  ;

3) полином Лежандра  имеет  различных и действительных корней, расположенных на интервале .

Составим по узлам интегрирования многочлен -й степени

 .

Функция  при  есть многочлен степени не выше . Значит для этой функции формула Гаусса справедлива:

 , (4.3)

так как  .

Разложим  в ряд по ортогональным многочленам Лежандра:

 ,

 ,

,

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

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

Зная , из линейной теперь системы первых  (4.2) легко найти коэффициенты . Определитель этой системы есть определитель Вандермонда.

Формулу , в которой - нули полинома Лежандра , а  определяют из (3.3), называют квадратурной формулой Гаусса.

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

Рис. 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));

 line(i,j,i,j1);

 {Сам график}

 setcolor(14);

Страницы: 1, 2, 3, 4, 5, 6, 7, 8, 9