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

Пусть задана табличная функция своими (n+1) узловыми значениями (x₀,y₀ ) ; (x₁,y₁ ) ; (x₁,y₁ ) ; …, (x_n,y_n ) . Необходимо вычислить коэффициенты a_i , i=(0,1,2,…n) интерполирующего многочлена (1).

Так как по определению интерполирующий многочлен проходит через узлы интерполяции без погрешности, то значения каждого из них обращает в тождество уравнение (1), то есть

 a₀ + a₁x₀ + a₂x₀² +…+ a_nx₀ⁿ = y_0,

a₀ + a₁x₁ + a₂x₁² +…+ a_nx₁ⁿ = y_1,

a₀ + a₁x₂ + a₂x₂² +…+ a_nx₂ⁿ = y_2,

a₀ + a₁x_n + a₂x_n² +…+ a_nx_nⁿ = y_n .

Если функция однозначна и x_i различны, то в результате подстановки значений узловых точек интерполяции в уравнения (1) получаем систему (2) (n+1) уравнений с (n+1) неизвестным. Определитель этой системы называют определителем Вандермонда:

1 x₀ x₀² … x₀ⁿ

1 x₁ x₁² … x₁ⁿ

= 1 x₂ x₂² … x₂ⁿ

… … … …

1 x_n x_n² … x_nⁿ

Определитель Вандермонда отличен от нуля, если x_j все различные. Следовательно, данная система всегда имеет решение, и притом единственное. Таким образом, доказано существование и единственность интерполяционного многочлена. Значения коэффициентов интерполирующего многочлена получают, решив систему (2).

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

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

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

На рисунке 2.3 представлена схема алгоритма возведения действительного числа в степень (подпрограмма-функция vozvedenie_v_stepenb).

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

На рисунке 2.5 представлена схема алгоритма построения СЛАУ с неизвестными - коэффициентами интерполяционного полинома (подпрограмма-процедура sistema).

На рисунке 2.6 представлена схема алгоритма решения полученной CЛАУ (подпрограмма-процедура reshenie).

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

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

Рисунок 2.1 - Блок-схема алгоритма решения задачи №2

Рисунок 2.2 - Блок-схема алгоритма ввода данных

Рисунок 2.3 - Блок-схема алгоритма возведения действительного числа в степень

Рисунок 2.4 - Блок-схема алгоритма вычисления определителя квадратной матрицы

Рисунок 2.5 - Блок-схема алгоритма построения СЛАУ с неизвестными - коэффициентами интерполяционного полинома

Рисунок 2.6 - Блок-схема алгоритма решения полученной CЛАУ

Рисунок 2.7 - Блок-схема алгоритма записи в файл данных и результата

Рисунок 2.8 - Блок-схема алгоритма вывода содержимого записанного файла на экран

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

program vandermond;

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;

vand:real;

n:integer;

fail,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;

function vozvedenie_v_stepenb(osnovanie:real;pokazatelb:integer):real;

{Функция возводит в степень заданное действительное число}

begin

if pokazatelb=0 then vozvedenie_v_stepenb:=1

else

vozvedenie_v_stepenb:=osnovanie*vozvedenie_v_stepenb(osnovanie,pokazatelb-1);end;

Function Determinante(Mas:Matr; kolvo:Integer):Real;

{Функция вычисляет определитель заданной матрицы}

var

Arr:Matr;

i,j,k:Integer;

sum:Real;

begin

For i:=0 to kolvo do

For j:=0 to kolvo do Arr[i,j]:=Mas[i,j];

sum:= 0;

if kolvo>1 then begin

for i:=0 to kolvo do

begin

for j:=0 to kolvo do

begin

if i<>j then

begin

for k:=0 to kolvo-1 do

begin

if j>i then Arr[k,j-1]:=mas[k+1,j] else

arr[k,j]:=mas[k+1,j];

end;

end

end;

if (i+1) mod 2 = 1 then sum:=sum+mas[0,i]*Determinante(Arr,kolvo-1)

else sum:=sum-mas[0,i]*Determinante(Arr,kolvo-1);

end;

determinante:=sum;

end

else begin

determinante:=mas[0,0]*mas[1,1]-mas[0,1]*mas[1,0];

end;

end;

procedure sistema(uzel,fun:mas; kolvo:integer; var a1:matr; var b1:mas; var opr:real);

{Процедура строит матрицы слау, необходимые для решения,

также в ней вычисляется определитель Вандермонда}

var i,j:integer;

begin

for i:=0 to kolvo do

for j:=0 to kolvo do

begin

a1[i,j]:=vozvedenie_v_stepenb(uzel[i],j);

end;

for i:=0 to kolvo do b1[i]:=fun[i];

opr:=determinante(a1,kolvo);

end;

procedure reshenie(a1:matr; b1:mas; opr:real; kolvo:integer; var koef:mas);

{В данной процедуре осуществляется поиск коэффициентов интерполяционного полинома }

var c1:matr;

i,j,k:integer;

begin

for k:=0 to kolvo do

begin

for i:=0 to kolvo do

for j:=0 to kolvo do

c1[i,j]:=0;

for i:=0 to kolvo do

for j:=0 to kolvo do

if j=k then c1[i,j]:=b1[i]

else c1[i,j]:=a1[i,j];

koef[k]:=determinante(c1,kolvo)/opr;

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

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;

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(Нажмите Enter);

readln;

sistema(X,Y,N,A,B,vand);

reshenie(a,b,vand,n,koef_polinoma);

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

vblvod(fail,ekran);

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

readln;

grafik(n,x,y,koef_polinoma);

END.

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

Построить интерполяционный полином по значениям функции в узлах

x0=-1 y0=-1; x1=0 y1=-1; x2=1 y2=-1; x3=2 y3=0

Решение. По данным узлам и значениям функции составим СЛАУ

c0-c1+c2-c3=-1;

c0=-1;

c0+c1+c2+c3=-1;

c0+2c1+4c2+8c3=0;

Решив данную систему, получим

c0=-1, с1=-0.1667, с2=0, с3=0.1667.

На рисунке 2.9 изображен тестовый пример работы программы для данной системы, на рисунке 2.10 - график полученной функции.

Рисунок 2.9 - Тестовый пример работы программы интерполирования функции, заданной в узлах, методом Вандермонда

Рисунок 2.10 - График функции, полученной в данном тестовом примере

2.8 Инструкция пользователю

Данная программа осуществляет интерполирование функции, заданной в узлах, методом Вандермонда. Вам будет предложено ввести количество узлов, узлы интерполяции и значения функции в них. Данные и результат записываются в файл interpol.txt, и его содержимое выводится на экран. После этого вы просмотрите график функции, для выхода из программы нажмите enter.

2.9 Инструкция программисту

Данная программа осуществляет интерполирование функции, заданной в узлах, методом Вандермонда.

Константы:

С - константа целого типа, обозначает верхнюю границу матрицы системы и свободных членов.

Типы:

matr=array[0..c,0..c] of real - двумерный массив с равным количеством строк и столбцов.

mas=array[0..c] of real - одномерный массив.

Переменные:

x,y:mas - узлы интерполяции и значения функции в них.

a:matr - матрица СЛАУ для нахождения коэффициентов полинома.

b:mas - вектор свободных членов СЛАУ.

vand:real - определитель Вандермонда.

n:integer - номер последнего узла (нумерация с нуля).

Fail:text-файловая переменная, соответствующая записываемому файлу;

Ekran:text-файл, в который копируется данный файл для вывода на экран.

Подпрограммы

1) procedure Vvod(var kolvo:integer; var uzel,fun:mas) - ввод исходных данных.

Формальные параметры:

-возвращаемые значения:

uzel,fun:mas-массивы узлов интерполяции и значений функции в них

kolvo:Integer-номер последнего узла(нумерация с нуля);

Локальные переменные:

code:integer - индикатор ошибки ввода данных(параметр процедуры перевода строковой переменной в число Val)

i:Integer-счетчик цикла;

s:string-переменная для контроля ввода данных .

2) Function Determinante(Mas:Matr; kolvo:Integer):Real - вычисление определителя квадратной матрицы.

Формальные параметры:

-передаваемые значения:

Mas:Matr-матрица, определитель которой надо вычислить;

kolvo:Integer-размерность матрицы;

Локальные переменные:

Arr:Matr-матрица, в которую копируются элементы данной (вспомогательная для вычисления определителя);

i,j,k:Integer-счетчики циклов;

sum:Real-результат(значение, возвращаемое функцией);

3) function vozvedenie_v_stepenb(osnovanie:real;pokazatelb:integer):real - возведение действительного числа в степень.

Формальные параметры:

-передаваемые значения:

Osnovanie - переменная типа Real, возвращающая значение основания степени.

Pokazatelb - переменная типа Integer, возвращающая значение показателя степени.

4) procedure sistema(uzel,fun:mas; kolvo:integer; var a1:matr; var b1:mas; var opr:real) - построение СЛАУ с неизвестными - коэффициентами интерполяционного полинома.

Формальные параметры:

-передаваемые значения:

uzel,fun:mas-массивы узлов интерполяции и значений функции в них

kolvo:Integer-номер последнего узла(нумерация с нуля);

-возвращаемые значения

a1:matr; b1:mas-матрица системы и вектор свободных членов;

opr:real-определитель Вандермонда.

Локальные переменные:

I,j:Integer-счетчики цикла;

5) procedure reshenie(a1:matr; b1:mas; opr:real; kolvo:integer; var koef:mas) - решение полученной СЛАУ.

Формальные параметры:

-передаваемые значения:

a1:matr; b1:mas-матрица системы и вектор свободных членов;

opr:real-определитель Вандермонда.

kolvo:Integer-номер последнего узла(нумерация с нуля);

-возвращаемые значения

koef:mas-искомые коэффициенты интерполяционного полинома

Локальные переменные:

c1:matr-матрица системы, в которой один из столбцов заменен на вектор свободных членов

I,j,k:Integer-счетчики цикла.

6) procedure zapisb(koef:mas; uzel,fun:mas; kolvo:integer; var f:text) - запись в файл данных и результата.

Формальные параметры:

-передаваемые значения:

uzel,fun:mas-массивы узлов интерполяции и значений функции в них

kolvo:Integer-номер последнего узла(нумерация с нуля);

koef:mas-коэффициенты интерполяционного полинома

-возвращаемые значения:

f:text-файловая переменная, соответствующая записываемому файлу

Локальные переменные:

i-счетчик в циклах.

7) procedure outputtoscreen(var imyafaila:string; var f1,f2:text) - вывод содержимого записанного файла на экран.

Формальные параметры:

-передаваемые значения:

imyafaila:string-имя файла, который выводится на экран;

-возвращаемые значения:

f1,f2:text - файловые переменные, соответствующие копируемому файлу и файлу, в который осуществляется копирование;

Локальные переменные:

s1:string-копируемая строка.

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

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

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

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

При интерполировании используется условие равенства значений интерполяционного многочлена и аппроксимируемой функции в узлах интерполяции. Для инженерной практики типична задача приближения функции с помощью многочлена по некоторому заданному базису {f _k(x), k=0,1,...,n}, когда значение приближаемой функции заданы в числе узлов, превышающих число базисных функций.

Будем нумеровать узлы от аргумента индексами х₀,...,х_N. При N>n практически всегда не найдется многочлена, значения которого совпадали бы со значениями аппроксимируемой функции во всех узлах. В этом случае задача аппроксимации функции полиномом ставится как задача поиска такого вектора коэффициентов многочлена

,

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

3.3 Обзор существующих методов решения

Среднеквадратичное приближение по тригонометрическому базису

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

(t)= (1)

где б_р, в_р - произвольные числовые коэффициенты. Тригонометрическими многочленами естественно приближать периодические функции с периодом 2? . В теории рядов Фурье устанавливается, что коэффициенты тригонометрического многочлена наилучшего среднеквадратичного приближения непрерывной 2? -периодичной функции f(t), для которого величина

принимает минимальное значение, задаются формулами

, , .

Рассмотрим задачу нахождения тригонометрического многочлена наилучшего среднеквадратичного приближения на дискретном множестве точек. Пусть N,m - натуральные числа, m? N/2. Даны точки

i=0,1,...,N.

и значения функции в этих точках f(t). Требуется аппроксимировать функцию f(t) на интервале тригонометрическим полиномом (1).

Коэффициенты тригонометрического многочлена наилучшего среднеквадратичного приближения функции f(t) на дискретном множестве точек х_i, i=0,1,...,N, в данном случае находятся по формулам

.

Для приближения функции на произвольном отрезке необходимо произвести замену переменных.

Мы будем решать поставленную задачу среднеквадратичного приближения функции по степенному базису.

3.4 Численный метод решения

Критерий среднеквадратического приближения (метод наименьших квадратов):

Коэффициенты многочлена выбирают исходя из условия

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

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

,

,

Чтобы эта система имела единственное решение необходимо и достаточно, чтобы определитель матрицы А (определитель Грамма) был отличен от нуля.

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

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

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

На рисунке 3.2 представлена схема алгоритма чтения из файла исходных данных (подпрограмма-процедура chtenie_dannblh).

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

На рисунке 3.4 представлена схема алгоритма построения базиса степенных функций (подпрограмма-процедура bazis).

На рисунке 3.5 представлена схема алгоритма построения СЛАУ с неизвестными - коэффициентами искомого полинома (подпрограмма-процедура sistema).

На рисунке 3.6 представлена схема алгоритма решения полученной CЛАУ (подпрограмма-процедура reshenie_systembl).

Рисунок 3.1 - Блок-схема алгоритма решения задачи №3

Рисунок 3.2 - Блок-схема алгоритма чтения из файла исходных данных

Рисунок 3.3 - Блок-схема алгоритма чтения из файла исходных данных (продолжение)

Рисунок 3.4 - Блок-схема алгоритма построения базиса степенных функций

Рисунок 3.5 - Блок-схема алгоритма получения СЛАУ, вектор неизвестных которой - искомые коэффициенты полинома

Рисунок 3.6 - Блок-схема алгоритма решения полученной СЛАУ

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

program priblizhenie;

uses crt,graph;

type mas1=array [0..20] of real;

mas2=array [0..2] of real;

mas3=array [0..2,0..2] of real;

mas4=array [0..2,0..20] of real;

var fail:text;

x,y:mas1;

koef_polinoma,b:mas2;

n:integer; psi:mas4;

a:mas3;

procedure chtenie_dannblh(var f:text; var uzlbl,funktsiya:mas1; var kolvo:integer);

{Процедура осуществляет чтение данных из файла, получаем узлы и значения функции в них}

var code,j,tab,i,dlina_stroki:integer;

x_s,y_s,s:string;

begin

kolvo:=-1;

assign(f,e:202020.txt);

reset(f);

while not eof(f) do

begin

readln(f,s);

dlina_stroki:=length(s);

kolvo:=kolvo+1;

x_s:= ;

y_s:= ;

for i:=1 to dlina_stroki do if s[i]=chr(9) then

tab:=i;

i:=1;

while i<tab do

begin

x_s:=x_s+s[i];

i:=i+1;

end;

for i:=tab+1 to dlina_stroki do y_s:=y_s+s[i];

delete(x_s,1,1);

delete(y_s,1,1);

val(x_s,uzlbl[kolvo],code);

val(y_s,funktsiya[kolvo],code);

end;

close(f);

for j:=0 to kolvo do writeln(uzel=,uzlbl[j]:4:2, function=,funktsiya[j]:16:15);