Математическая модель для решения задачи

В формуле Гаусса на каждом интервале интегрирования значение функции f(x) вычисляется не в равномерно распределенных по интервалу узлах, а в абсциссах, выбранных из условия обеспечения минимума погрешности интерполяции:

где n- число интервалов интегрирования, m - число вычисляемых на каждом интервале значений функции. , - границы интервалов интегрирования; и - коэффициенты значения которых определяются величиной m. Для m=3 A1=5/9, A2=8/9, A3=5/9, , t2=0, t3=-t1

Блок-схема программы

Текст программы

Uses

Crt, Graph;

const

A1=5/9;

A2=8/9;

A3=5/8;

t=-0.77459;

t2=0;

t3=-t; {константы для взятия интеграла методом Гаусса}

var

GraphDriver, GraphMode : Integer;

x, x_nach, x_kon, eps: real; {пределы интегрирования и точность вычисления}

MaxX, MaxY: integer;

y_min_FUN, y_max_FUN: real;

yFUN: array [1..1100] of real; {массив для хранения точек графика}

c: char;

i, j, n: integer; {указатели на эл.массива и т.д.}

{ * * * * * * * функция y(x) * * * * * * * * * * * *--}

Function FUN(x: real): real;

begin

FUN := x*sin(x)/(1+x*x);

end;

{ * * * * * * *--вычисление интеграла методом Гаусса * * * * *}

function Gauss (a, b, eps: real; var k:integer): real;

var

S, z, h, c, d, l, x, x1, x2, x3: real;

i, n: integer;

begin

n:=1;

S:=0;

k:=0;

repeat

k:=k+1; {увеличиваем число итераций}

z:=S; {предыдущее значение интеграла равно текущему}

n:=n*2; {в два раза увеличиваем число интервалов интегрирования}

h:=(b-a)/n;

x:=a;

S:=0;

c:=h/2;

l:=c*t; {определение шага интегрирования,

начального значения x, сам интеграл сначала равен 0,

вспомогательные переменные считаем }

for i:=0 to n-1 do {перебираем все интервалы интегрирования}

begin

d:=x+c;

x1:=d-l;

x2:=d;

x3:=d+l; {вычисляем значения абсцисс узлов,

выбранных из условия обеспечения минимума погрешности интерполяции}

S:=S+A1*(FUN(x1)+FUN(x3))+A2*FUN(x2); {добавляем к сумме}

x:=x+h; {переходим на новый интервал интегрирования}

end;

S:=S*c; {умножаем полученную сумму на h/2}

until (abs(z-S)<eps*abs(S)) or (k>=14);

Gauss := S;

end;

{ * * * * * *--построение Декарт. системы координат * * *-}

procedure Dekart (x_min, x_max, y_min, y_max: real);

var

x_i, y_i: real;

s: string;

i : real;

ii: integer;

nulx: integer;

koord_x, koord_y: integer;

begin

GraphDriver := Detect; {инициализация графики}

InitGraph(GraphDriver, GraphMode, '');

if GraphResult <> grOk then Halt(1);

nulx := round (MaxY-MaxY*(0-y_min)/(y_max - y_min));

setcolor(DarkGray);

Setlinestyle(1,1,1);

i := 0;

while i <= MaxY

do begin

line (0, round(i), MaxX, round(i)); {деление экрана на 10 частей}

i := i + (MaxY/20);

end;

i := 0;

while i <= MaxX

do begin

line (round(i), 0, round(i), MaxY);

i := i + (MaxX/20);

end;

Setlinestyle(0,1,1);

setcolor (white);

x_i := x_min;

for ii := 1 to 10 {разметка по осям}

do begin

koord_x := round (MaxX*(x_i-x_min)/(x_max-x_min));

koord_y := nulx+10;

Str (x_i:1:2, s);

OutTextXY (koord_x, koord_y, s);

x_i := x_i + (x_max-x_min)/10;

end;

y_i := y_min;

for ii := 1 to 10

do begin

koord_x := 5;

koord_y := round (MaxY-MaxY*(y_i-y_min)/(y_max - y_min));

Str (y_i:1:2, s);

OutTextXY (koord_x, koord_y-10, s);

y_i := y_i + (y_max - y_min)/10;

end;

setcolor (white);

OutTextXY (MaxX-15, nulx-12, 'X');

OutTextXY (MaxX-20, nulx-4, ' >');

OutTextXY (0, 4, '^ Y');

line (0, nulx, MaxX, nulx);

line (2, 0, 2, MaxY);

end;

{ * * *--расчет точек, мин и макс для графика функции *--}

Procedure CalcF (x_nach, x_kon: real);

var

x: real;

dx: real;

i: integer;

begin

x := x_nach;

i := 1;

dx := (x_kon-x_nach)/(MaxX-1);

y_min_FUN := FUN(x); y_max_FUN := FUN(x);

while x <= x_kon do

begin

yFUN[i] := FUN(x);

if yFUN[i] < y_min_FUN

then y_min_FUN := yFUN[i];

if yFUN[i] > y_max_FUN

then y_max_FUN := yFUN[i];

x := x + dx;

i := i + 1

end

end;

{ * * * *--построение функции y(x) * * * * * * * *}

Procedure GFUN;

var

i: integer;

koord_x1, koord_y1, koord_x2, koord_y2: integer;

begin

setcolor (White);

OutTextXY (40, 40, ' sin(x) ');

OutTextXY (40, 48, ' Y(x) = x* * * . Для продолжения нажмите ввод...');

OutTextXY (40, 56, ' (1+x^2) ');

for i := 1 to MaxX-2 do

begin

koord_x1 := i;

koord_y1 := round (MaxY-MaxY*(yFUN[i]-y_min_FUN)/(y_max_FUN - y_min_FUN));

koord_x2 := i+1;

koord_y2 := round (MaxY-MaxY*(yFUN[i+1]-y_min_FUN)/(y_max_FUN - y_min_FUN));

line (koord_x1, koord_y1, koord_x2, koord_y2)

end

end;

{======================================================================}

{ * * * * * * * * *--главная программа * * * * * * * *}

begin

GraphDriver := Detect;

InitGraph(GraphDriver, GraphMode, '');

if GraphResult <> grOk then Halt(1);

MaxX := GetMaxX;

MaxY := GetMaxY;

CloseGraph;

repeat

repeat

clrscr;

writeln ('Введите значения для вычисления графика:');

write (' X нач. = ');

readln (x_nach);

write (' Х кон. = ');

readln (x_kon);

write (' Точность = ');

readln (eps);

until (x_nach>=0) and (x_kon <= 22) and (x_nach<x_kon);

writeln ('Интеграл = ', Gauss ( x_nach, x_kon, eps, n):7:5,' ',n,' итераций');

writeln ('Для продолжения нажмите ввод...');

readln;

CalcF (x_nach, x_kon);

Dekart (x_nach, x_kon, y_min_FUN, y_max_FUN);

GFUN;

readln;

CloseGraph;

GotoXY(23,10);

writeln ('Для продолжения нажмите пробел,');

GotoXY(23,12);

writeln ('для выхода из программы - ESC.');

c := ReadKey;

until ord(c) = 27;

End.