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

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

2.3 Обзор существующих численных методов решения задачи

Интерполяция по Лагранжу

Интерполяционный многочлен может быть построен при помощи специальных интерполяционных формул Лагранжа, Ньютона, Стерлинга, Бесселя и др.

Интерполяционный многочлен по формуле Лагранжа имеет вид:

Докажем, что многочлен Лагранжа является интерполяционным многочленом, проходящим через все узловые точки, т.е. в узлах интерполирования xi выполняется условие Ln(xi) = yi. Для этого будем последовательно подставлять значения координат узловых точек таблицы в многочлен (2.1). В результате получим:

если x=x0, то Ln(x0) = y0,

если x=x1, то Ln(x1) = y1,

……………

если x=xn, то Ln(xn) = yn.

Это достигнуто за счет того, что в числителе каждой дроби при соответствующем значении уj, j=0,1,2,…,n отсутствует сомножитель (x-xi), в котором i=j, а знаменатель каждой дроби получен заменой переменной х на соответствующее значение хj.

Таким образом, интерполяционный многочлен Лагранжа приближает заданную табличную функцию, т.е. Ln(xi) = yi и мы можем использовать его в качестве вспомогательной функции для решения задач интерполирования, т.е. .

Чем больше узлов интерполирования на отрезке [x0,xn] , тем точнее интерполяционный многочлен приближает заданную табличную функцию, т.е. тем точнее равенство:

Однако с увеличением числа узлов интерполирования возрастает степень интерполяционного многочлена n и в результате значительно возрастает объем вычислительной работы. Поэтому при большом числе узлов необходимо применять ЭВМ. В этом случае удобно находить значения функции в промежуточных точках, не получая многочлен в явном виде.

При решении задачи экстраполирования функции с помощью интерполяционного многочлена вычисление значения функции за пределами отрезка [x0,xn] обычно производят не далее, чем на один шаг h, равный наименьшей величине

так как за пределами отрезка [x0,xn] погрешности, как правило, увеличиваются.

Интерполяция по Ньютону

Интерполяционный многочлен по формуле Ньютона имеет вид:

(2.2)

где n - степень многочлена,

- разделенные разности 0-го, 1-го, 2-го,…., n-го порядка, соответственно.

Сплайн-интерполяция

Сплайны стали широко использоваться в вычислительной математике сравнительно недавно. В машиностроительном черчении они применяются уже давно, так как сплайны - это лекала или гибкие линейки, деформация которых позволяет провести кривую через заданные точки (xi, уi).

Используя теорию изгиба бруса при малых деформациях, можно показать, что сплайн - это группа кубических многочленов, в местах сопряжения которых первая и вторая производные непрерывны. Такие функции называются кубическими сплайнами. Для их построения необходимо задать коэффициенты, которые единственным образом определяют многочлен в промежутке между данными точками.

Например, для некоторых функций (рис.) необходимо задать все кубические функции q1(x), q2(x), …qn(x).

В наиболее общем случае эти многочлены имеют вид:

где kij - коэффициенты, определяемые описанными ранее условиями, количество которых равно 4n. Для определения коэффициентов kij необходимо построить и решить систему порядка 4n.

Первые 2n условий требуют, чтобы сплайны соприкасались в заданных точках:

Следующие (2n-2) условий требуют, чтобы в местах соприкосновения сплайнов были равны первые и вторые производные:

Система алгебраических уравнений имеет решение, если число уравнений соответствует числу неизвестных. Для этого необходимо ввести еще два уравнения. Обычно используются следующие условия:

При построении алгоритма метода первые и вторые производные удобно аппроксимировать разделенными разностями соответствующих порядков.

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

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

Значения f(x0), f(x1), … , f(xn) , т.е. значения табличной функции в узлах, называются разделенными разностями нулевого порядка (k=0).

Отношение называется разделенной разностью первого порядка (k=1) на участке [x0, x1] и равно разности разделенных разностей нулевого порядка на концах участка [x0, x1], разделенной на длину этого участка.

Для произвольного участка [xi, xi+1] разделенная разность первого порядка (k=1) равна

Отношение называется разделенной разностью второго порядка (k=2) на участке [x0, x2] и равно разности разделенных разностей первого порядка, разделенной на длину участка [x0, x2].

Для произвольного участка [xi, xi+2] разделенная разность второго порядка (k=2) равна

Таким образом, разделенная разность k-го порядка на участке [xi, xi+k] может быть определена через разделенные разности (k-1)-го порядка по рекуррентной формуле:

(2.3)

Где n - степень многочлена.

Максимальное значение k равно n. Тогда i =0 и разделенная разность n-го порядка на участке [x0,xn] равна

,

т.е. равна разности разделенных разностей (n-1)-го порядка, разделенной на длину участка [x0,xn].

Разделенные разности

являются вполне определенными числами, поэтому выражение (2.2) действительно является алгебраическим многочленом n-й степени. При этом в многочлене (2.2) все разделенные алгебраическим многочленом n-й степени. При этом в многочлене (2.2) все разделенные разности определены для участков [x0, x0+k], .

Лемма: алгебраический многочлен (2.2), построенный по формулам Ньютона, действительно является интерполяционным многочленом, т.е. значение многочлена в узловых точках равно значению табличной функции

Докажем это. Пусть х=х0 , тогда многочлен (2.2) равен

Пусть х=х1, тогда многочлен (2.2) равен

Пусть х=х2, тогда многочлен (2.2) равен

Заметим, что решение задачи интерполяции по Ньютону имеет некоторые преимущества по сравнению с решением задачи интерполяции по Лагранжу. Каждое слагаемое интерполяционного многочлена Лагранжа зависит от всех значений табличной функции yi, i=0,1,…n. Поэтому при изменении количества узловых точек N и степени многочлена n (n=N-1) интерполяционный многочлен Лагранжа требуется строить заново. В многочлене Ньютона при изменении количества узловых точек N и степени многочлена n требуется только добавить или отбросить соответствующее число стандартных слагаемых в формуле Ньютона (2.2). Это удобно на практике и ускоряет процесс вычислений.

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

На рисунке 2.1 представлена схема алгоритма решения задачи №2.

На рисунке 2.2 представлена схема алгоритма ввода исходных данных (подпрограмма-процедура Vvod).

На рисунке 2.3 представлена схема алгоритма интерполяции функции по методу Ньютона с разделенными разностями (newt)

На рисунке 2.4 представлена схема алгоритма записи данных и результата в файл (подпрограмма-процедура zapisb_v_fail).

На рисунке 2.5 представлена схема алгоритма вывода содержимого записанного файла на экран (подпрограмма-процедура outputtoscreen).

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

program newton;

uses crt,graph;

const c=10;

type matr=array[0..c,0..c] of real;

mas=array[0..c] of real;

var x,y,koef_polinoma:mas;

a:matr;

b:mas;

d1:real;

n:integer;

fail,fail1,ekran:text;

procedure Vvod(var kolvo:integer; var uzel,fun:mas);

{Процедура осуществляет ввод данных:пользователь вводит с клавиатуры

узлы интерполяции и значения функции в них. Также определяется количество узлов.}

var code,i:integer; s:string;

begin

writeln(введите количество узлов);

readln(kolvo);

kolvo:=kolvo-1;

for i:=0 to kolvo do

begin

repeat

writeln(введите ,i,-й узел интерполирования);

readln(s);

val(s,uzel[i],code);

until code=0;

repeat

writeln(введите значение функции, соответствующее данному узлу);

readln(s);

val(s,fun[i],code);

until code=0;

end;

end;

procedure newt(var kolvo:integer; D:real; var koef,uzel,fun:mas)

var L,P:real;

begin

L:=fun[0];

P:=1;

for i:=1 to kolvo do

begin

P:=P*(D-uzel[i-1]);

for j:=1 to kolvo-i do

begin

fun[j]:=(fun[j-1]-fun[j])/(uzel[j+i]-uzel[i])

end;

koef[i]:=fun[0];

L:=L+P*fun[0];

end;

end;

procedure newt(var kolvo:integer; D:real; var koef,uzel,fun:mas) {процедура интерполяции функции методом Ньютона}

var L,P:real;

begin

L:=fun[0];

P:=1;

for i:=1 to kolvo do

begin

P:=P*(D-uzel[i-1]);

for j:=1 to kolvo-i do

begin

fun[j]:=(fun[j-1]-fun[j])/(uzel[j+i]-uzel[i])

end;

koef[i]:=fun[0];

L:=L+P*fun[0];

end;

end;

procedure zapisb(koef:mas; uzel,fun:mas; kolvo:integer; var f:text);

{В данной процедуре осуществляется запись в файл данных и результата}

var i:integer;

begin

assign(f,interpol.txt);

rewrite(f);

for i:=0 to kolvo do writeln(f,x= ,uzel[i]:8:4, f(x)=,fun[i]:8:4);

writeln(f,Интерполяционный полином);

write(f,p(x)=,koef[0]:8:4);

for i:=1 to kolvo do if i>1 then write (f,+(,koef[i]:8:4,)*x^,i)

else write (f,+(,koef[i]:8:4,)*x);

close(f);

end;

procedure vblvod(var f1,f2:text);

{Вывод содержимого записанного файла на экран}

var s1:string;

begin

clrscr;

assign(f1,interpol.txt);

reset(f1);

assigncrt(f2);

rewrite(f2);

while not eof(f1) do

assigncrt(f2);

rewrite(f2);

while not eof(f1) do

begin

Readln(f1,s1);

writeln(f2,s1);

end;

close(f2);

close(f1);

end;

procedure grafik(kolvo:integer; uzlbl,funktsiya:mas; c:mas);

{Построение графика полученной функции}

var driver,mode,Err,a1,b1,z,i,j:integer; s:string; xt,yt:real;

begin

driver:=detect;

InitGraph(driver,mode,d: p7pgi);

err:=graphresult;

if err<>grok then writeln(Ошибка при инициализации графического режима)

else

begin

Setcolor(9);

line(320,0,320,480);

line(0,240,640,240);

settextstyle(smallfont,horizdir,3);

setcolor(10);

outtextxy(320,245,0);

a1:=0;

b1:=480;

z:=-10;

for i:=0 to 20 do

begin

if z<>0 then

begin

str(z,s);

setcolor(10);

outtextxy(a1,245,s);

outtextxy(325,b1,s);

setcolor(8);

line(0,b1,640,b1);

line(a1,0,a1,480);

end;

outtextxy(325,b1,s);

setcolor(8);

line(0,b1,640,b1);

line(a1,0,a1,480);

end;

a1:=a1+32;

b1:=b1-24;

z:=succ(z);

end;

setcolor(5);

for i:=0 to kolvo do

begin

xt:=uzlbl[i];

yt:=funktsiya[i];

putpixel(round(320+xt*32),round(240-yt*24),5);

end;

moveto(round(320+uzlbl[0]*32),round(240-funktsiya[0]*24));

setcolor(11);

for i:=0 to kolvo do

begin

xt:=uzlbl[i];

yt:=0;

for j:=0 to kolvo do yt:=yt+c[j]*vozvedenie_v_stepenb(uzlbl[i],j);

lineto(round(320+xt*32),round(240-yt*24));

moveto(round(320+xt*32),round(240-yt*24));

end;

readln;

closegraph;

end;

end;

{Основная часть программы}

BEGIN

CLRSCR;

Writeln(Программа осуществляет интерполирование функции, заданной в узлах);

Vvod(N,X,Y);

writeln(`введите значение промежуточной точки);

readln(d1);

2.7 Тестовый пример

writeln(Нажмите Enter);

readln;

newt(N,d1,X,Y,koef_polinoma);

zapisb(koef_polinoma,x,y,n,fail);

vblvod(fail,fail1);

writeln(Нажмите Enter для просмотра графика функции, затем еще раз для выхода из программы);

readln;

grafik(N,X,Y,koef_polinoma);

END.

readln(d1);

writeln(Нажмите Enter);

readln;

newt(N,d1,X,Y,koef_polinoma);

zapisb(koef_polinoma,x,y,n,fail);

vblvod(fail,fail1);

writeln(Нажмите Enter для просмотра графика функции, затем еще раз для выхода из программы);

readln;

grafik(N,X,Y,koef_polinoma);

END.

Дана табличная функция:

Вычислить разделенные разности 1-го, 2-го, 3-го порядков (n=3) и занести их в диагональную таблицу.

Разделенные разности первого порядка:

Разделенные разности второго порядка:

Разделенная разность третьего порядка:

Интерполяционный многочлен Ньютона для заданной табличной функции имеет вид:

График интерполяционного многочлена будет таким:

procedure zapisb(koef:mas; uzel,fun:mas; kolvo:integer; var f:text);

{В данной процедуре осуществляется запись в файл данных и результата}

var i:integer;

begin

assign(f,interpol.txt);

rewrite(f);

for i:=0 to kolvo do writeln(f,x= ,uzel[i]:8:4, f(x)=,fun[i]:8:4);

writeln(f,Интерполяционный полином);

write(f,p(x)=,koef[0]:8:4);

for i:=1 to kolvo do if i>1 then write (f,+(,koef[i]:8:4,)*x^,i)

else write (f,+(,koef[i]:8:4,)*x);

close(f);

end;

procedure vblvod(var f1,f2:text);

{Вывод содержимого записанного файла на экран}

var s1:string;

begin

clrscr;

assign(f1,interpol.txt);

reset(f1);

assigncrt(f2);

rewrite(f2);

while not eof(f1) do

begin

Readln(f1,s1);

writeln(f2,s1);

end;

close(f2);

close(f1);

end;

procedure grafik(kolvo:integer; uzlbl,funktsiya:mas; c:mas);

{Построение графика полученной функции}

var driver,mode,Err,a1,b1,z,i,j:integer; s:string; xt,yt:real;

begin

driver:=detect;

InitGraph(driver,mode,d: p7pgi);

err:=graphresult;

if err<>grok then writeln(Ошибка при инициализации графического режима)

else

begin

Setcolor(9);

line(320,0,320,480);

line(0,240,640,240);

settextstyle(smallfont,horizdir,3);

setcolor(10);

outtextxy(320,245,0);

a1:=0;

b1:=480;

z:=-10;

for i:=0 to 20 do

begin

if z<>0 then

3. Среднеквадратическое приближение функции

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

Разработать схему алгоритма и написать программу на языке Turbo Pascal 7.0 для выполнения среднеквадратического приближения функции, заданной в узлах.

3.2 Математическая формулировка задачи

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

, (3.1)

где - весовая функция.

Такое приближение называют среднеквадратичным.

Делись добром ;)