Файл: ПОСОБИЕ_ВычМат.pdf

47

Замечание. В практических задачах возможна ситуация, когда Dkyi = 0, для некоторого 1< k £ n. Тогда степень полинома равна k–1 < n.

3.2.2.1 Алгоритмы интерполирования многочленами Ньютона

Рассматриваемые алгоритмы интерполирования основаны на формулах (7) и

(8). Для вычисления разностей различных порядков удобно использование следующего рекуррентного соотношения:

в i-ом слагаемом формулы (7) может быть вычислен с использованием предыдущего множителя, а именно

q( q − 1)( q − 2 )...( q −( i − 2 )) × q −( i − 1) . ( i - 1)! i

1444442444443

предыдущий множитель

Аналогичное соотношение имеет место в формуле (8).

Алгоритм расчета по первой интерполяционной формуле Ньютона

1.Ввести значения: n (n+1 – число узлов интерполяции), x0 – начальное значение аргумента, h – шаг интерполирования, yi (i = 0, 1, …, n) –значение функции в i-ом узле интерполяции, x – значение аргумента для интерполирования функции.

2.Вычислить q = (x – x0)/h.

3.Вычислить Pn(x) = y0 – учитываем первое слагаемое в формуле (7).

4.Определить p = 1.

5.Организовать цикл интерполирования, изменяя параметр i =1, 2, … n.

5.1. Вычислить p = p× q −( i − 1) – множитель (10). i

5.2.Вычислить Diy0 по формуле (9).

5.3.Вычислить Pn(x) = Pn(x) + p×Diy0 – суммирование i-го слагаемого в фор-

муле (7).

6.Конец цикла интерполирования по параметру i.

7.Вывод x – значение аргумента, Pn(x) – результат интерполирования функции.

8.Конец вычислений.

Алгоритм расчета по второй интерполяционной формуле Ньютона можно получить модификацией вышеописанного алгоритма, в частности:

1.Вычислить q = (x – xn)/h, где xn = x0 + nh.

2.Вычислить Pn(x) = yn.

48

2.1. Вычислить p = p× q + ( i − 1) . i

Пример 2. Составить программы на языке Turbo Pascal интерполирования функции f(x) на отрезке многочленами Ньютона по системе 5 узлов в точках: M1(1.05, 2.8577), M2(1.09, 2.9743), M3(1.13, 3.0957), M4(1.17, 3.2220), M5(1.21, 3.3535).

{первая интерполяционная формула Ньютона} program pNuton1;

uses crt; const

n = 4; {число узлов интерполяции, счет от нуля} type

Tmas = array[0..n] of real;

{*******функция вычисляет k-ю разность yi} function fdelta( var y : Tmas; i, k : byte) : real; begin

if k = 0 then fdelta := y[i] else fdelta := fdelta(y, i+1, k-1) - fdelta(y,i,k-1) end;

{основная программа} const

y : Tmas = (2.8577, 2.9743, 3.0957, 3.2220, 3.3535); var

i : byte; x0 , x , h , p , pn : real; delta, q : real; begin

x0:=1.05; h := 0.04;

writeln('Введите значение аргумента x для выполнения интерполяции'); readln(x);

pn:=y[0]; q:=(x-x0)/h; p:=1;

for i:=1 to n do {цикл интерполирования}

begin p:=p*(q-(i-1))/i; delta:=fdelta(y,0,i); pn:=pn+p*delta; end; writeln('Для значения x = ', x:12:5, ' получен результат Pn =', pn:12:5 ); writeln('Точное значение: ', exp(x):12:5);

readkey;

end.

{вторая интерполяционная формула Ньютона} program pNuton2;

uses crt; const n=4;

type Tmas=array[0..n] of real; {*******функция вычисляет k-ю разность yi}

49

function fdelta( var y:Tmas;i,k:byte):real; begin

if k = 0 then fdelta := y[i] else fdelta := fdelta(y, i+1, k-1) - fdelta(y, i, k-1) end;

{основная программа} const

y : Tmas = (2.8577, 2.9743, 3.0957, 3.2220, 3.3535); var

i : byte; x0, x, h, p, pn, delta, xn, q : real; begin

h:=0.04; x0:=1.05;

writeln('Введите значение аргумента x для выполнения интерполирования'); readln(x);

pn := y[n]; xn := x0 + n*h; q := (x-xn)/h; p:=1; for i:=1 to n do

begin

p:=p*(q+(i-1))/i; delta:=fdelta(y,n-i,i);{вычисляет i-ю разность y[n-i]} pn:=pn+p*delta;

end;

writeln('x=',x:12:5,' Pn =',pn:12:5); writeln('Точное значение = ',exp(x):12:5); readkey;

end.

3.2.3 Интерполирование функций средствами MachCAD

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

Для осуществления сплайновой интерполяции предлагается четыре встроенные функции. Три из них служат для получения векторов вторых производных сплайн-функций при различном виде интерполяции:

§cspline (X, Y) – возвращает вектор S вторых производных при приближении в опорных точках к кубическому полиному;

§pspline (X, Y) – возвращает вектор S вторых производных при приближении в опорных точках к параболической кривой;

§lspline (X, Y) – возвращает вектор S вторых производных при приближении в опорных точках к прямой.

Четвертая функция - interp (S, X, Y, z) – возвращает значение у(х) для заданных векторов S, X, Y и заданного значения z.

Выполните в системе MachCAD следующие действия:

1. Организуйте данные в виде отсортированных соразмерных векторов:

2. Задайте вектор вторых производных при помощи одной из приведенных выше функций семейства *spline, например lspline:

s1 := lspline(X ,Y)

3.Задайте функцию interp(s1,x,y, z). В нашем случае:

I1(Z) := interp(lspline(X ,Y) , X ,Y ,Z)

4.Постройте график, выбрав пиктограмму на панели Graph. Пропишите в соответствующих маркерах имя функции интерполяции и переменную. Получили сплайн-интерполяцию с линейным продолжением:

5.Аналогичным образом постройте сплайн-интерполяцию, используя функции cspline, pspline. Сравните результаты.

Задания к данной теме приведены в приложении В (1).