Размещено на http://www.allbest.ru/

Министерство образования Российской Федерации

Сибирский Государственный Технологический Университет

Факультет: Автоматизации и информационных технологий

Кафедра: Системотехники

Курсовая работа

Полиномиальная интерполяция Гаусса, Ньютона, Стирлинга

Руководитель: Ващенко Г.В.

Разработал: Русинов Д.В.

Содержание

Аннотация

Введение

1. Постановка задачи

2. Описание метода полиномиальной интерполяции

3. Блок-схема программного обеспечения

4. Исходные тексты основных процедур программы

5. Результаты численных экспериментов

Заключение

Список используемых источников

Аннотация

В данной курсовой работе представлены алгоритмы и программное обеспечение, реализующее решение полиномиального интерполирования методами Ньютона, Гаусса и Стирлинга.

Программное обеспечение разработано в среде программирования Delphi.

Введение

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

p(x_i)=y_i (i=0,1,2, …, n)

1. Постановка задачи

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

2. Описание метода полиномиальной интерполяции

Пусть функция y=f(x) задана на сетке равноотстоящих узлов

x_i=x_o+ih,

где i=0, 1,…, n, и для неё построена таблица конечных разностей.

Будем строить интерполяционный многочлен P_n(x) в форме:

P_n(x)=a₀+a₁(x-x₀)+a₂(x-x₀)(x-x₁)+…+a_n(x-x₀)(x-x₁)…(x-x_n-1)

3. Конструктивная схема программного обеспечения

1. Функция ввода:

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

procedure TForm1.Button1Click(Sender: TObject);

begin

x0:=strtofloat(EditXo.text);

h:=(Strtofloat(EditXn.Text)-x0)/n;

For i:=0 to n do

y[i]:=strtofloat(setka.Cells[i+1,1]);

PN1:=PolyNuton1;

PN2:=PolyNuton2;

for i:=0 to n do

y[i-(n div 2)]:=y[i];

x0:=x0+(n div 2)*h;

PG1:=PolyGauss1;

PG2:=PolyGauss2;

PS:=PolyStirling;

With Polynoms do

begin

Cells[1,0]:=PN1; Cells[1,1]:=PN2;

Cells[1,2]:=PG1; Cells[1,3]:=PG2;

Cells[1,4]:=PS;

DrawGraph;

end;

2. Функция вычисления:

Function PolyNuton1:string;

Function PolyNuton2:string;

Function PolyGauss1:string;

функции вычисления полиномов в текстовом виде

Function PolyGauss2:string;

Function PolyStirling:string;

Function PolyN1(x:extended):extended;

Function PolyN2(x:extended):extended;

Function PolyG1(x:extended):extended;

функции вычисления полиномов как функции

Function PolyG2(x:extended):extended;

Function PolyS(x:extended):extended;

Функция вывода:

procedure TForm1.Button1Click(Sender: TObject);

begin

x0:=strtofloat(EditXo.text);

h:=(Strtofloat(EditXn.Text)-x0)/n;

For i:=0 to n do

y[i]:=strtofloat(setka.Cells[i+1,1]);

PN1:=PolyNuton1;

PN2:=PolyNuton2;

for i:=0 to n do

y[i-(n div 2)]:=y[i];

x0:=x0+(n div 2)*h;

PG1:=PolyGauss1;

PG2:=PolyGauss2;

PS:=PolyStirling;

With Polynoms do

begin

Cells[1,0]:=PN1; Cells[1,1]:=PN2;

Cells[1,2]:=PG1; Cells[1,3]:=PG2;

Cells[1,4]:=PS;

end;

DrawGraph;

end;

Procedure DrawGraph - процедура рисования графиков

4. Исходные тексты основных процедур программы

unit Unit1;

interface

uses

Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,

Dialogs, Grids, Math, StdCtrls, Spin,StrUtils, ExtCtrls;

type

TArray=array[-100..200] of extended;

TForm1 = class(TForm)

Setka: TStringGrid;

SpinEditn: TSpinEdit;

Label1: TLabel;

Label2: TLabel;

EditXo: TEdit;

Label3: TLabel;

EditXn: TEdit;

Button1: TButton;

Polynoms: TStringGrid;

Image1: TImage;

CheckN1: TCheckBox;

CheckN2: TCheckBox;

CheckG1: TCheckBox;

CheckG2: TCheckBox;

CheckS: TCheckBox;

procedure Button1Click(Sender: TObject);

procedure FormCreate(Sender: TObject);

procedure EditXnExit(Sender: TObject);

procedure EditXoExit(Sender: TObject);

procedure SpinEditnKeyPress(Sender: TObject; var Key: Char);

procedure SpinEditnExit(Sender: TObject);

procedure PolynomsKeyPress(Sender: TObject; var Key: Char);

procedure CheckN1Click(Sender: TObject);

private

Function Fact(n:longint):longint; //факториал

Function C(n,k:integer):extended; //С из n по k

Function Delta(k:integer;y:Tarray;i:integer):extended; //Дельта к-ой степ. Yi

Procedure ZapX; //заполнить строку с x-ми

Function PolyNuton1:string;

Function PolyNuton2:string;

Function PolyGauss1:string;

Function PolyGauss2:string;

Function PolyStirling:string;

Function PolyN1(x:extended):extended;

Function PolyN2(x:extended):extended;

Function PolyG1(x:extended):extended;

Function PolyG2(x:extended):extended;

Function PolyS(x:extended):extended;

Procedure DrawGraph;

public

{ Public declarations }

end;

var

Form1: TForm1;

y:TArray;

n,i,j,k:integer;

h,x0,z,x,r:extended;

PN1,PN2,PG1,PG2,PS:string;

implementation

{$R *.dfm}

Procedure TForm1.DrawGraph;

var

scx,scy,miny,maxy:extended;

i:integer;

begin

x0:=strtofloat(EditXo.text);

h:=(Strtofloat(EditXn.Text)-x0)/n;

For i:=0 to n do

y[i]:=strtofloat(setka.Cells[i+1,1]);

scx:=(image1.Width-30)/n;

miny:=y[0];

maxy:=miny;

for i:=1 to n do

begin

If miny>y[i] then miny:=y[i];

If maxy<y[i] then maxy:=y[i];

end;

scy:=(Image1.Height-30)/(maxy-miny);

Image1.Canvas.FillRect(Rect(0,0,Width,Height));

With Image1.Canvas do

begin

Pen.Color:=clLime;

Pen.Style:=psSolid;

Pen.Width:=1;

for i:=0 to n do

begin

Moveto(20+round(scx*i),Image1.Height-15);

Lineto(20+round(scx*i),0);

TextOut(15+round(scx*i),Image1.Height-15,floatToStr(x0+i*h));

Moveto(15,Image1.Height-20-round(scy*(y[i]-miny)));

LineTo(Image1.Width,Image1.Height-20-round(scy*(y[i]-miny)));

TextOut(5,Image1.Height-27-round(scy*(y[i]-miny)),floatToStr(y[i]));

end;

end;

With Image1.Canvas do

begin

Pen.Color:=clblack;

Pen.Width:=3;

Moveto(20,Image1.Height-20-round(scy*(y[0]-miny)));

for i:=1 to n do

LineTo(20+round(scx*i),Image1.Height-20-round(scy*(y[i]-miny)));

pen.Width:=2;

For i:=21 to Image1.Width-10 do

begin

If CheckN1.Checked then

begin

Pen.Color:=clBlue;

Moveto(i-1,Image1.Height-20-round(scy*(PolyN1(x0+(i-21)*h/scx)-miny)));

Lineto(i,Image1.Height-20-round(scy*(PolyN1(x0+(i-20)*h/scx)-miny)));

end;

If CheckN2.Checked then

begin

Pen.Color:=clGreen;

Moveto(i-1,Image1.Height-20-round(scy*(PolyN2(x0+(i-21)*h/scx)-miny)));

Lineto(i,Image1.Height-20-round(scy*(PolyN2(x0+(i-20)*h/scx)-miny)));

end;

end;

for i:=0 to n do

y[i-(n div 2)]:=y[i];

x0:=x0+(n div 2)*h;

For i:=21 to Image1.Width-10 do

begin

If CheckG1.Checked then

begin

Pen.Color:=clRed;

Moveto(i-1,Image1.Height-20-round(scy*(PolyG1(x0-(n div 2)*h+(i-21)*h/scx)-miny)));

Lineto(i,Image1.Height-20-round(scy*(PolyG1(x0-(n div 2)*h+(i-20)*h/scx)-miny)));

end;

If CheckG2.Checked then

begin

Pen.Color:=rgb(50,200,150);

Moveto(i-1,Image1.Height-20-round(scy*(PolyG2(x0-(n div 2)*h+(i-21)*h/scx)-miny)));

Lineto(i,Image1.Height-20-round(scy*(PolyG2(x0-(n div 2)*h+(i-20)*h/scx)-miny)));

end;

If CheckS.Checked then

begin

Pen.Color:=clgray;

Moveto(i-1,Image1.Height-20-round(scy*(PolyS(x0-(n div 2)*h+(i-21)*h/scx)-miny)));

Lineto(i,Image1.Height-20-round(scy*(PolyS(x0-(n div 2)*h+(i-20)*h/scx)-miny)));

end;

end;

end;

end;

Procedure TForm1.ZapX;

begin

x0:=strtofloat(EditXo.text);

h:=(Strtofloat(EditXn.Text)-x0)/n;

For i:=1 to n+1 do

Setka.Cells[i,0]:=floattostr(x0+pred(i)*h);

end;

Function TForm1.Fact(n:integer):longint;

var i:integer;

begin

Result:=1;

For i:=1 to n do

Result:=Result*i;

end;

function Tform1.C(n,k:integer):extended;

begin

Result:=fact(n)/(fact(k)*fact(n-k));

end;

Function Tform1.Delta(k:integer;y:TArray;i:integer):extended;

var j:integer;

begin

result:=0;

For j:=0 to k do

Result:=result+intpower(-1,j)*C(k,j)*y[k+i-j];

end;

Function TForm1.PolyNuton1:string;

begin

Result:=floattostr(y[0]);

For i:=1 to n do

begin

z:=Delta(i,y,0)/(fact(i)*IntPower(h,i));

if z=0 then continue;

If (z>0) then Result:=Result+'+';

If z<>1 then Result:=Result+floattostr(z);

For j:=1 to i do

begin

z:=x0+pred(j)*h;

if z=0

then Result:=Result+'x'

else if z>0 then Result:=Result+'(x-'+floattostr(z)+')'

else Result:=Result+'(x'+floattostr(z)+')'

end;

end;

end;

Function Tform1.PolyN1(x:extended):extended;

begin

Result:=(y[0]);

For i:=1 to n do

begin

z:=Delta(i,y,0)/(fact(i)*IntPower(h,i));

if z=0 then continue;

r:=z;

For j:=1 to i do

begin

z:=x0+pred(j)*h;

r:=r*(x-z)

end;

Result:=Result+r;

end;

end;

function TForm1.PolyNuton2:string;

begin

Result:=floattostr(y[n]);

For i:=1 to n do

begin

z:=Delta(i,y,n-i)/(fact(i)*IntPower(h,i));

if z=0 then continue;

If (z>0) then Result:=Result+'+';

If z<>1 then Result:=Result+floattostr(z);

For j:=n downto n-pred(i) do

begin

z:=x0+j*h;

if z=0

then Result:=Result+'x'

else if z>0 then Result:=Result+'(x-'+floattostr(z)+')'

else Result:=Result+'(x'+floattostr(z)+')'

end;

end;

end;

Function Tform1.PolyN2(x:extended):extended;

begin

Result:=(y[n]);

For i:=1 to n do

begin

z:=Delta(i,y,n-i)/(fact(i)*IntPower(h,i));

if z=0 then continue;

r:=z;

For j:=n downto n-pred(i) do

begin

z:=x0+j*h;

r:=r*(x-z)

end;

Result:=Result+r;

end;

end;

function TForm1.PolyGauss1:string;

begin

Result:=floattostr(y[0]);

For i:=1 to n do

begin

z:=Delta(i,y,-(i div 2))/(fact(i)*IntPower(h,i));

if z=0 then continue;

If (z>0) then Result:=Result+'+';

If z<>1 then Result:=Result+floattostr(z);

For j:=-(i div 2)+1-i mod 2 to (i div 2) do

begin

z:=x0+j*h;

if z=0

then Result:=Result+'x'

else if z>0 then Result:=Result+'(x-'+floattostr(z)+')'

else Result:=Result+'(x'+floattostr(z)+')'

end;

end;

end;

Function Tform1.PolyG1(x:extended):extended;

begin

Result:=(y[0]);

For i:=1 to n do

begin

z:=Delta(i,y,-(i div 2))/(fact(i)*IntPower(h,i));

if z=0 then continue;

r:=z;

For j:=-(i div 2)+1-i mod 2 to (i div 2) do

begin

z:=x0+j*h;

r:=r*(x-z)

end;

Result:=Result+r;

end;

end;

function TForm1.PolyGauss2:string;

begin

Result:=floattostr(y[0]);

For i:=1 to n do

begin

z:=Delta(i,y,-(succ(i) div 2))/(fact(i)*IntPower(h,i));

if z=0 then continue;

If (z>0) then Result:=Result+'+';

If z<>1 then Result:=Result+floattostr(z);

For j:=-(i div 2) to (i div 2)-1+i mod 2 do

begin

z:=x0+j*h;

if z=0

then Result:=Result+'x'

else if z>0 then Result:=Result+'(x-'+floattostr(z)+')'

else Result:=Result+'(x'+floattostr(z)+')'

end;

end;

end;

Function Tform1.PolyG2(x:extended):extended;

begin

Result:=(y[0]);

For i:=1 to n do

begin

z:=Delta(i,y,-(succ(i) div 2))/(fact(i)*IntPower(h,i));

if z=0 then continue;

r:=z;

For j:=-(i div 2) to (i div 2)-1+i mod 2 do

begin

z:=x0+j*h;

r:=r*(x-z)

end;

Result:=Result+r;

end;

end;

function TForm1.PolyStirling:string;

begin

Result:=floattostr(y[0]);

For i:=1 to n do

begin

z:=(Delta(i,y,-(i div 2))+Delta(i,y,-(succ(i) div 2)))/(2*fact(i)*IntPower(h,i));

If not(odd(i)) then z:=z/2;

if z=0 then continue;

If (z>0) then Result:=Result+'+';

If z<>1 then Result:=Result+floattostr(z);

For j:=-(i div 2)+1-i mod 2 to (i div 2)-1+i mod 2 do

begin

z:=x0+j*h;

if z=0

then Result:=Result+'x'

else if z>0 then Result:=Result+'(x-'+floattostr(z)+')'

else Result:=Result+'(x'+floattostr(z)+')'

end;

If not(odd(i)) then

begin

z:=-2*x0;

If z=0 then Result:=result+'2x'

else

if z<0

then Result:=Result+'(2x'+floattostr(z)+')'

else Result:=Result+'(2x+'+floattostr(z)+')';

end;

end;

end;

Function Tform1.PolyS(x:extended):extended;

begin

Result:=(y[0]);

For i:=1 to n do

begin

z:=(Delta(i,y,-(i div 2))+Delta(i,y,-(succ(i) div 2)))/(2*fact(i)*IntPower(h,i));

If not(odd(i)) then z:=z/2;

if z=0 then continue;

r:=z;

For j:=-(i div 2)+1-i mod 2 to (i div 2)-1+i mod 2 do

begin

z:=x0+j*h;

r:=r*(x-z)

end;

If not(odd(i)) then r:=r*(2*(x-x0));

Result:=Result+r;

end;

end;

procedure TForm1.Button1Click(Sender: TObject);

begin

x0:=strtofloat(EditXo.text);

h:=(Strtofloat(EditXn.Text)-x0)/n;

For i:=0 to n do

y[i]:=strtofloat(setka.Cells[i+1,1]);

PN1:=PolyNuton1;

PN2:=PolyNuton2;

for i:=0 to n do

y[i-(n div 2)]:=y[i];

x0:=x0+(n div 2)*h;

PG1:=PolyGauss1;

PG2:=PolyGauss2;

PS:=PolyStirling;

With Polynoms do

begin

Cells[1,0]:=PN1; Cells[1,1]:=PN2;

Cells[1,2]:=PG1; Cells[1,3]:=PG2;

Cells[1,4]:=PS;

end;

DrawGraph;

end;

procedure TForm1.FormCreate(Sender: TObject);

begin

n:=4;

With Polynoms do

begin

ColWidths[1]:=1593;

Cells[0,0]:='I полином Ньютона';

Cells[0,1]:='II полином Ньютона';

Cells[0,2]:='I полином Гаусса';

Cells[0,3]:='II полином Гаусса';

Cells[0,4]:='полином Стирлинга';

end;

With setka do

begin

Cells[0,0]:=' x';

Cells[0,1]:=' y';

zapX

end;

end;

procedure TForm1.EditXnExit(Sender: TObject);

begin

If StrTofloat(EditXn.Text)<StrTofloat(EditXo.Text) then

EditXn.Undo;

ZapX;

end;

procedure TForm1.EditXoExit(Sender: TObject);

begin

If StrTofloat(EditXn.Text)<StrTofloat(EditXo.Text) then

EditXo.Undo;

ZapX;

end;

procedure TForm1.SpinEditnKeyPress(Sender: TObject; var Key: Char);

begin

If not(key in ['0'..'9','-',#8,DecimalSeparator]) then key:=#0;

end;

procedure TForm1.SpinEditnExit(Sender: TObject);

begin

setka.ColCount:=SpinEditn.Value+1;

n:=SpinEditn.Value-1;

ZapX;

end;

procedure TForm1.PolynomsKeyPress(Sender: TObject; var Key: Char);

begin

key:=#0;

procedure TForm1.CheckN1Click(Sender: TObject);

begin

DrawGraph;

end;

end.5. Результаты численных экспериментов.

Хорошо обусловленные матрицы(Ввод пользователем)

№	Размерность матрицы	Вектор невязки(r)	Норма матрицы ||r||-равномерное
1	3	r1=0,00001 r2=0 r3=0	||r||=101
2	4	r1=0 r2=0,00002 r3=0,00003 r4=0,00003	||r||=153
3	10	r1=0,00002 r2=0,00004 r3=0,00003 r4=0,00006 r5=0,00001 r6=0,00003 r7=0,00012 r8=0,00001 r9=0,00004 r10=0,00021	||r||=362

Хорошо обусловленные матрицы(случайный выбор)

№	Размерность матрицы	Вектор невязки(r)	Норма матрицы ||r||-равномерное	Диапазон случайной величины
1	3	r1=0 r2=0 r3=0	||r||=13	От 1 до 5
2	4	r1=0 r2=0 r3=0 r4=0	||r||=36	От 1 до10
3	10	r1=0,00014 r2=0,00004 r3=0,00013 r4=0,00009 r5=0,00032 r6=0,00001 r7=0,00012 r8=0,00021 r9=0,00002 r10=0,00022	||r||=472	От 1 до 100

Плохо обусловленные матрицы(Ввод пользователем)

№	Размерность матрицы	Вектор невязки(r)	Норма матрицы ||r||-равномерное
1	3	r1=0,00001 r2=0,00021 r3=0,00003	||r||=1,45
2	4	r1=0,00012 r2=0,00002 r3=0,00005 r4=0,00021	||r||=2,67
3	10	r1=0,00003 r2=0,00032 r3=0,00023 r4=0,00078 r5=0,00005 r6=0,00036 r7=0,00011 r8=0,00054 r9=0,00262 r10=0,00050	||r||=3,49

Плохо обусловленные матрицы(заполнение по правилу 1/I+j)

№	Размерность матрицы	Вектор невязки(r)	Норма матрицы ||r||-равномерное
1	3	r1=0 r2=0 r3=0	||r||=1,83
2	4	r1=0 r2=0 r3=0 r4=0	||r||=2,08
3	10	r1=0 r2=0 r3=0 r4=0 r5=0 r6=0 r7=0 r8=0 r9=0 r10=0	||r||=2,92

Заключение

программный полиномиальный интерполяция линейный

В данной работе был разработан алгоритм и программа для решения уравнения вида A*х=B методом оптимального исключения. Метод оптимального исключения дает возможность решения систем линейных уравнений размерностью на два порядка выше, чем в обычных методах (как например методы Гаусса, Жордана-Гаусса и т. д), за счет возможности подзагрузки в оперативную память порции уравнений во время решения системы.

Численные эксперименты, проведенные над линейными системами с хорошо и плохо обусловленными матрицами, различных размерностей, показали, что, не смотря на сложность вычислительной схемы (алгоритма) метод дает ряд преимуществ перед обычными методами (методы Гаусса, Жордана-Гаусса и т. д), а именно уменьшается погрешность (вектор невязки) и возрастает скорость вычисления.

Численные эксперименты проведены на процессоре Pentium 4 с одноразрядной точностью(real).

Список используемых источников

1. Фаддеев Д.К., Фаддеева В.Н. Вычислительные методы линейной алгебры.-М.: Физматгиз, 1963.- 734 с.

2. Культин Н.Б. Основы программирования 0в Delphi 7. - СПб.: БХВ-Петербург, 2003. -608с.: ил.

3. Бобровский И.С. Delphi 7. Учебный курс. - СПб.: Питер, 2003. - 736 с.: ил.

Размещено на Allbest.ru