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

Министерство образования Российской Федерации Сибирский Государственный Технологический Университет Факультет: Автоматизации и информационных технологий Кафедра: Системотехники Курсовая работа Полиномиальная интерполяция Гаусса, Ньютона, Стирлинга Руководитель: Ващенко Г. В.

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

Содержание

Аннотация Введение

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

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)

EditXn.Undo;

ZapX;

end;

procedure TForm1. EditXoExit (Sender: TObject);

begin

If StrTofloat (EditXn.Text)

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. Результаты численных экспериментов.

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

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

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

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

Заключение

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

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

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

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

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

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

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

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