Банк Рефератов - Вычисление интеграла методом Ньютона-Котеса (теория и программа на Паскале)

Главная: Рефераты

Рефераты. Вычисление интеграла методом Ньютона-Котеса (теория и программа на Паскале)

end;

readln;

setcolor(15);

end;

procedure graphik(ea:word;a,b:real;f1:string);

{процедура построения графиков}

var

f,f2:string;

d:char;

i,v,r:integer;

x1,x2,n,p,x:integer;

c,k,k1:longint;

y:array[0..1] of double;

begin

x1:=-240;

x2:=240;

c:=24;

setcolor(15);

n:=0;v:=0;r:=0;

repeat

cleardevice;

settextstyle(0,0,0);

if ea mod 2 =0 then

begin

outtextxy(10,1,'Нажмите F1 для помощи');

str(c/24:2:2,f);

f:='Масштаб '+f+':1';

end

else

begin

outtextxy(10,1,'Press F1 for help');

str(c/24:2:2,f);

f:='Scale '+f+':1';

end;

outtextxy(200,1,f);

settextstyle(3,0,1);

outtextxy(307,10,'y');

outtextxy(574,235,'x');

outtextxy(310,240,'0');

setlinestyle(1,7,100);

line(70,240,580,240);

line(320,20,320,460);

line(320,20,315,25);

line(321,20,326,25);

line(580,239,575,244);

line(580,240,575,235);

line(70,239,580,239);

line(321,20,321,460);

for i:=-9 to 10 do

begin

if ((320+i*24)71) then

line(320+i*24,240,320+i*24,242);

if ((240+i*24)19) then

line(320,240+i*24,322,240+i*24);

end;

setcolor(15);

for x:= -240+round((240+x1)/10) to 240+round((240+x1)/10) do

begin

funktia(1,x-1,x,y,c,f1);

k:=round(240-(y[0])*c);

k1:=round(240-(y[1])*c);

if ((k0)or(k10)) then

line(319-round((240+x1)/10)+x,k,320-round((240+x1)/10)+x,k1);

end;

if (v mod 2)=0 then

begin

funktia(1,a,b,y,1,f1);

k:=round(240-(y[0])*c);

k1:=round(240-(y[1])*c);

line(320-round((240+x1)/10)+round(a*c),k,320-

round((240+x1)/10)+round(a*c),240);

line(320-round((240+x1)/10)+round(b*c),k1,320-

round((240+x1)/10)+round(b*c),240);

if 320-round((240+x1)/10)+a*c560 then

begin

funktia(1,(-240-round((240+x1)/10))/c,(240-

round((240+x1)/10))/c,y,1,f1);

k1:=round(240-(y[1])*c);

line(560,k1,560,240);

end;

for x:= -240 to 240 do

begin

funktia(1,x-1,x,y,c,f1);

k1:=round(240-(y[1])*c);

if ((x/c)>a) and ((x/c)0) then

begin

if (abs(240-k1)>2) then

begin

if k17 then

setfillstyle(6,3)

else

setfillstyle(1,3);

floodfill(320-round((240+x1)/10)+x,k1,15);

end;

str(x1,f2);

outtextxy(1,450,f2);

if (n mod 2)=0 then

for i:=-9 to 10 do

begin

settextstyle(2,0,2);

setcolor(14);

if ((320+i*24)71)and(i<>0) then

begin

str((i*24+round((240+x1)/10))/c:2:2,f);

p:=247;

outtextxy(310+i*24,p,f);

str(-i*24/c:2:2,f);

outtextxy(330,240+i*24,f);

end;

for i:=-9 to 10 do

begin

setcolor(15);

if ((r mod 2)=1) and (i<>0) then

begin

if ((320+i*24)71) then

line(320+i*24,20,320+i*24,460);

if ((240+i*24)19) then

line(80,240+i*24,560,240+i*24);

end;

setcolor(15);

d:=readkey;

case d of

#75:

begin

x1:=x1-30;

x2:=x2-30;

end;

#77:

begin

x1:=x1+30;

x2:=x2+30;

end;

#80:

if c>1 then

c:=c-1;

#72:

c:=c+1;

#71:

c:=24;

#79:

n:=n+1;

#83:

r:=r+1;

#82:

v:=v+1;

#73:

begin

c:=24;

n:=0;r:=0;v:=0;x1:=-240;x2:=240;

end;

#59:

hwg(ea);

end;

until d=#13;

end;

end.

================================================

==========МОДУЛЬ UNIT==========

================================================

{$N+}

Unit k_unit;

{Модуль нахождения интеграл от многочлена q(q-1)..(q-i+1)(q-i-1)..(q-n),}

{где n-точность интеграла ,i-номер коофициента. }

interface

procedure rasposn(f:string;x:real;var ec:word;var t:real);

procedure hkoef(n:integer;var h:array of double);

procedure funktia(n:integer;a,b:real;var y:array of

double;c:real;f:string);

procedure koef(w:array of double;n:integer;var e:array of double);

procedure mnogochlen(n,i:integer;var c:array of double);

function facktorial(n:integer):double;

function integral(w:array of double;n:integer):double;

function mainint(n:integer;a,b:real;y:array of double):double;

implementation

procedure rasposn(f:string;x:real;var ec:word;var t:real);

{Процедура распознования функции}

var

k:word;

begin

k:=pos('x',f);

if k<>0 then

begin {Распознавание функции}

ec:=1; {Код ошибки}

t:=x;

k:=pos('abs(x)',f);

if k<>0 then t:=abs(x);

k:=pos('sin(x)',f);

if k<>0 then t:=sin(x);

k:=pos('cos(x)',f);

if k<>0 then t:=cos(x);

k:=pos('arctg(x)',f);

if k<>0 then t:=arctan(x);

k:=pos('sqr(x)',f);

if k<>0 then t:=x*x;

k:=pos('exp(x)',f);

if k<>0 then t:=exp(x);

k:=pos('cos(x)*x',f);

if k<>0 then t:=cos(x)*x;

k:=pos('ln(x)',f);

if k<>0 then

begin

if x>0 then t:=ln(x)

else

t:=0;

end;

k:=pos('sqrt(x)',f);

if k<>0 then

if x>=0 then t:=sqrt(x)

else t:=0;

k:=pos('arcctg(x)',f);

if k<>0 then t:=pi/2-arctan(x);

k:=pos('sin(x)/x',f);

if k<>0 then if x<>0 then t:=sin(x)/x;

end

else

ec:=0;

end;

procedure funktia(n:integer;a,b:real;var y:array of

double;c:real;f:string);

{Процедур подсчет Y-ков и распознавания функции}

var

t,h,x:real;

k,i:integer;

es:word;

begin

h:=(b-a)/n;

for i:=0 to n do

begin

x:=(a+h*i)/c;

rasposn(f,x,es,t);

y[i]:=t;

end;

procedure koef(w:array of double;n:integer;var e:array of double);

{Изменение коофициентов для интеграла}

var

t:integer;

begin

for t:=1 to n do

e[t]:=w[t]/(n-t+2);

end;

procedure mnogochlen(n,i:integer;var c:array of double);

{процедура нахождения коофициентов при Q^n(q в степени n )}

var

k,j:integer;

d:array[1..100] of double;

begin

d[1]:=1;

for j:=1 to n do

begin {Вычисление коэффициентов при раскрытии q*(q-1)*(q-2)*..*(q-n)}

d[j+1]:=d[j]*j*(-1);

if j>1 then

for k:=j downto 2 do

d[k]:=d[k]+d[k-1]*j*(-1);

end;

c[1]:=d[1]; {Деление многочлена на (q-i) по схеме Горнера}

for j:=1 to n+1 do

c[j]:=i*c[j-1]+d[j];

koef(c,n,c); {Изменение коэффициентов при интегрировании}

end;

function facktorial(n:integer):double;

{функция нахождения факториала }

var

t:integer;

s:double;

begin

s:=1;

if n=0 then

s:=1

else

for t:=1 to n do

s:=s*t;

facktorial:=s;

end;

function integral(w:array of double;n:integer):double;

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

var

t,p:integer;

s,c:double;

begin

s:=0;p:=n;

for t:=0 to p+1 do

s:=s+w[t]*exp((p-t+2)*ln(p)); {Подсчет интеграла}

integral:=s;

end;

procedure hkoef(n:integer;var h:array of double);

{Процедура подсчета коэф. Ньютона-Котеса}

var

p,j,d,c,i:integer;

kq:array[0..20] of double;

s:array[0..20] of double;

begin

p:=n;

if (p mod 2)=1 then {Вычисление половины от всех вычислений

коэффициентов}

d:=round((p-1)*0.5)

else

d:=round(0.5*p);

for i:=0 to n do

begin

mnogochlen(p,i,kq);

s[i]:=integral(kq,p); {Формирование массива из интегралов}

end;

for i:=0 to d do

begin

if ((p-i) mod 2) = 0 then

c:=1

else

c:=(-1);

h[i]:=(c*s[i])/(facktorial(i)*facktorial(p-i)*p);

h[p-i]:=h[i];

end;

function mainint(n:integer;a,b:real;y:array of double):double;

{функция подсчета основного интеграла}

var

sum:double;

p,i:integer;

kq,h:array[0..20] of double;

begin

p:=n;

hkoef(n,h);

sum:=0;

for i:=0 to p do

sum:=sum+h[i]*y[i]; {Сумма произведений y-ков на коэффициенты}

mainint:=sum*(b-a);

end;

end.

================================================

=======ОСНОВНАЯ ПРОГРАММА=======

================================================

{$N+}

program Newton_Cotes_metod;{Программа нахождения определенного интеграла}

uses {методом Ньютона-Котеса }

k_unit,k_graph,graph,crt;

const

t=15;

var

c:char;

a1,b1,a,b:real;

n1,v,r,n:integer;

h,y:array[0..t] of double;

ea,k:word;

int:double;

f:string;

begin

ea:=10;

v:=detect;

initgraph(v,r,'');

cleardevice;

newsc(ea);

winwin1;

setcolor(15);

outtextxy(380,430,'Нажмите F2 для смены языка.');

repeat

win1(ea);

settextstyle(3,0,1);

outtextxy(178,340,'Press Enter...');

delay(13000);

bar(178,340,350,365);

delay(13000);

if keypressed then {Смена языка}

begin

c:=readkey;

if c=#60 then

begin

ea:=ea+1;

newsc(ea);

winwin1;

setcolor(15);

if ea mod 2 =0 then

outtextxy(380,430,'Нажмите F2 для смены языка.')

else

outtextxy(380,430,'Press F2 key to change language.');

end;

until c=#13;

repeat

newsc(ea);

win2(ea,k); {Ввод способа задания функции}

case k of

wwod1(ea,y,n,a,b);

begin

wwod2(ea,ea,n1,a1,b1,f);

n:=n1;a:=a1;b:=b1;

k:=4;

end;

if k=4 then

funktia(n,a,b,y,1,f);

int:=mainint(n,a,b,y); {Вычисление интеграла}

hkoef(n,h);

proline(ea);

win3(ea,n,a,b,int,f,h,k); {Последнее меню вывода результатов}

until k<>4;

closegraph;

end.

[pic]

Рассмотрим результаты тестовых испытаний для функций sin(x) на

интервале [-5;3] и exp(x) на интервале [2;8]

| |n=1 |n=2 |n=3 |n=4 |n=5 |n=7 |

|Sin(x) |4,040017 |3,02112 |0,087629 |1,779012 |1,537481 |1,246 |

|Exp(x) |8965,041 |3581,999 |3271,82 |3002,908 |2990,644 |2974,322 |

|N=9 |n=12 |

|1,273561 |1,27366 |

|2973,593 |2973,569 |

Видно, что при увеличении числа узлов интерполяции точность растет,

однако при больших n (n>15) наблюдался обратный эффект.

Рекомендуемый диапозон n: от 7 до 13.

1) Интерфейс программы составлен на 2 языках: русском и английском. Переход

с одного языка на другой осуществляется в вводном окне путем нажатия

клавишы F2. Сменить язык можно только в этой части программы.

2) При вводе значений функции вручную необходимо вводить только цифры и

после каждого ввода нажимать клавишу ENTER.

3) При испытании программы под разные операционные системы(Dos, Windows 98-

2k,NT, из под паскаля) происходил непонятный баг с неверным выводом на

экран значений коэффициентов Ньютона-Котеса, хотя интеграл считался

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

Dos.

4) При вводе параметров в “Меню задания параметров нахождения интеграла”

желательно их вводить постепенно сверху вниз, т.е. сначала ввести

количество узлов интерполяции, затем пределы интегрирования, а уж потом

вводить саму функцию.

5) Данная версия программы не способна распознавать все функции. Она может

распознать только стандартные функции Турбо Паскаля и еще несколько

дотполнительных: sin(x)/x, cos(x)*x ,arcctg(x). Для работы со

специфическими функциями необходимо в модуле K-unit в процедуре RASPOSN в

конце, перед end else, добавить :

k:=pos(‘Формула f(x)’,f);

if k<>0 then t:= ‘Формула f(x)’;

где ‘Формула f(x)’ – желаемая формула для распознования.

6) Вся помощь по вводу и работе с пограммой выводится в окне помощи.

Для нахождения интеграла существует много методов, однако, метод

Ньютона-Котеса один из самых быстрых: достаточно знать значения

коэффициентов для n=4, чтобы с точностью до сотых мгновенно посчитать

интеграл. Быстрота и простота –главные части этого метода.

В.И. Грызлов «Турбо Паскаль 7.0» Москва: ДМК 2000г.

Данилина «Численные методы» Москва: Высшая школа 1978г.

-----------------------

[pic]

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