Чисельні методи: розв`язування типових математичних задач, розв`язування систем лінійних рівнянь алгебри методом простих ітерацій

Зміст

• Введення
• 1. Розв`язування систем лінійних рівнянь алгебри методом простих ітерацій
• 1.1 Постановка завдання
• 1.2 Математичне формулювання завдання
• 1.3 Огляд існуючих чисельних методів рішення задачі
• 1.4 Чисельний метод рішення задачі
• 1.5 Схема алгоритму
• 1.6 Текст програми
• 1.7 Тестовий приклад
• 2. Поліноміальна інтерполяція функції методом Ньютона з розділеними різницями
• 2.1 Постановка завдання
• 2.2 Математичне формулювання завдання
• 2.3 Огляд існуючих чисельних методів рішення задачі
• 2.4 Чисельний метод рішення задачі
• 2.5 Текст програми
• 2.6 Тестовий приклад
• 3. Середньоквадратичне наближення функції
• 3.1 Постановка завдання
• 3.2 Математичне формулювання завдання
• 3.3 Огляд існуючих чисельних методів рішення задачі
• 3.4 Чисельний метод рішення задачі
• 3.5 Схема алгоритму
• 3.6 Текст програми
• 3.7 Тестовий приклад
• Висновок
• Список використаної літератури
• Введення
• Поява і безперервне вдосконалення швидкодіючих електронних обчислювальних машин (ЕОМ) привела до достовірно революційного перетворення павуки взагалі і математики особливо. Змінилася технологія наукових досліджень, колосально збільшилися можливості теоретичного вивчення, прогнозу складних процесів, проектування інженерних конструкцій. Рішення крупних науково-технічних проблем, прикладами яких можуть служити проблеми оволодіння ядерною енергією і освоєння космосу, стало можливим лише завдяки застосуванню математичного моделювання і нових чисельних методів, призначених для ЕОМ.
• В даний час можна говорити, що з'явився новий спосіб теоретичного дослідження складних процесів, що допускають математичний опис, - обчислювальний експеримент, тобто дослідження природничонаукових проблем засобами обчислювальної математики. Пояснимо істоту цього способу дослідження на прикладі рішення якої-небудь фізичної проблеми. Хай потрібно вивчити деякий фізичний процес. Математичному дослідженню передує вибір фізичного наближення, тобто рішення питання про те, які чинники треба врахувати, а якими можна нехтувати. Після цього проводиться дослідження проблеми методом обчислювального експерименту.
• Розробка і дослідження обчислювальних алгоритмів і їх застосування до рішення конкретних задач складає зміст величезного розділу сучасної математики -- обчислювальної математики.
• Обчислювальну математику визначають в широкому сенсі цього терміну як розділ математики, що включає круг питань, зв'язаних з використанням ЕОМ, і у вузькому сенсі -- як теорію чисельних методів і алгоритмів рішення поставлених математичних задач.
• Загальним для всіх чисельних методів є зведення математичного завдання до конечномерной. Це найчастіше досягається дискретизацією початкового завдання, тобто переходом від функцій безперервного аргументу до функцій дискретного аргументу. Після дискретизації початкового завдання треба побудувати обчислювальний алгоритм, тобто вказати послідовність арифметичних і логічних дій, виконуваних па ЕОМ і дій, що дають за кінцеве число, рішення дискретної задачі. Одержане рішення дискретної задачі береться за наближене рішення початкової математичної задачі.
• При рішенні задачі па ЕОМ ми завжди одержуємо не точне рішення початкової задачі, а деяке наближене рішення. Чим же обумовлена виникаюча погрішність? Можна виділити три основні причини виникнення погрішності при чисельному рішенні початкової математичної задачі. Перш за все, вхідні дані початкового завдання (початкові і граничні умови, коефіцієнти і праві частини рівнянь) завжди задаються з деякою погрішністю. Погрішність чисельного методу, обумовлену неточним завданням вхідних даних, прийнято називати неусувною погрішністю. Далі, при заміні початкового завдання дискретним завданням виникає погрішність, звана погрішністю дискретизації або, інакше, погрішністю методу. Нарешті, кінцева розрядність чисел, що представляються в ЕОМ, приводить до помилок округлення, які можуть наростати в процесі обчислень
• Чисельні методи дають наближене рішення задачі. Це означає, що замість точного рішення і (функції або функціонала) деякого завдання ми знаходимо рішення у іншого завдання, близьке в деякому розумінні (наприклад, по нормі) до шуканого. Основна ідея всіх методів -- дискретизація або апроксимація (заміна, наближення) початкового завдання іншим завданням, зручнішим для вирішення на ЕОМ, причому рішення апроксимуючої задачі залежить від деяких параметрів, управляючи якими, можна визначити рішення з необхідною точністю. Наприклад, в завданні чисельної інтеграції такими параметрами є вузли і ваги формули, квадратури. Далі, рішення дискретної задачі є елементом конечномерного простору.

1. Рішення систем лінійних рівнянь алгебри методом простої

ітерації

1.1 Постановка завдання

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

1.2 Математичне формулювання завдання

Хай А - невироджена матриця і потрібно вирішити систему

де діагональні елементи матриці А ненульові.

1.3 Огляд існуючих чисельних методів рішення задачі

Метод Гауса

У методі Гауса матриця СЛАР за допомогою рівносильних перетворень перетвориться у верхню трикутну матрицю, що виходить в результаті прямого ходу. У зворотному ході визначаються невідомі.

Хай дана СЛАР

Запишемо розширену матрицю системи:



На першому кроці алгоритму Гауса виберемо діагональний елемент (якщо він рівний 0, то перший рядок переставляємо з яким-небудь нижнім рядком) і оголошуємо його a₁₁_0 ведучим, а відповідний рядок і стовпець, на перетині яких він стоїть - ведучими. Обнулимо елементи провідного стовпця. Для цього сформуємо числа (-a₂₂/a₁₁) (-a₃₁/a₁₁) .., (a_n1/a₁₁).

LU-розкладання матриць

Комп'ютерна реалізація методу Гауса часто здійснюється з використанням LU-розкладання матриць.

LU - розкладання матриці A є розкладанням матриці A в твір нижньої і верхньої трикутних матриць, тобто

A=LU

де L - нижня трикутна матриця (матриця, у якої всі елементи, що знаходяться вище за головну діагональ рівні нулю lij=0 при i>j), U- верхня трикутна матриця (матриця, у якої всі елементи, що знаходяться нижче за головну діагональ рівні нулю uij=0 при i>j).

LU - розкладання може бути побудовано з використанням описаного вище методу Гауса. Розглянемо до - ий крок методу Гауса, на якому здійснюється обнулення піддіагональних елементів до - го стовпця матриці. Як було описано вище, з цією метою використовується наступна операція



В термінах матричних операцій така операція еквівалентна множенню A(k)=MkA(k-1), де елементи матриці визначаються таким чином

В термінах матричних операцій така операція еквівалентна множенню A(k)=MkA(k-1), де елементи матриці визначаються таким чином

В результаті прямого ходу методу Гауса одержимо, A(n-1)=U

де A(n-1)=U - верхня трикутна матриця, а - нижня трикутна матриця, що має вигляд .

Таким чином, шукане розкладання A=LU одержано.

Метод прогону

Метод прогону є одним з ефективних методів рішення СЛАР з трьох - діагональними матрицями, виникаючих при звичайно-різницевій апроксимації завдань для звичайних диференціальних рівнянь (ОДУ) і рівнянь в приватних похідних другого порядку і є окремим випадком методу Гауса. Розглянемо наступну СЛАР:

рішення якої шукатимемо у вигляді

де Qi,Pi,i=1,n - прогоночні коефіцієнти, що підлягають визначенню. Для їх визначення виразимо з першого рівняння СЛАР (1.1) x1 через x2, одержимо:

звідки

З другого рівняння СЛАР (1.1) з допомогою (1.3) виразимо x2 через x3, одержимо:

звідки

Продовжуючи цей процес, одержимо з i-го рівняння СЛАР (1.1):

Отже:

З останнього рівняння СЛАР маємо:

Метод простих ітерацій

При великому числі рівнянь прямі методи рішення СЛАР (за винятком методу прогону) стають важкореалізованими на ЕОМ перш за все із-за складності зберігання і обробки матриць великої розмірності. В той же час характерною особливістю ряду тих, що часто зустрічаються в прикладних завданнях СЛАР є розрідженість матриць. Число ненульових елементів таких матриць мало в порівнянні з їх розмірністю. Для вирішення СЛАР з розрідженими матрицями переважно використовувати ітераційні методи.

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

Метод Зейделя рішення СЛАР

Метод простих ітерацій досить поволі сходиться. Для його прискорення існує метод Зейделя, що полягає в тому, що при обчисленні компоненту xik+1 вектора невідомих на (k+1) -ій ітерації використовуються x1k+1, x2k+1, .., xik+1 вже обчислені на (k+1) -ій ітерації. Значення інших компонент беруться з попередньої ітерації. Так само як і в методі простих ітерацій будується еквівалентна СЛАР і за початкове наближення береться вектор правих частин . Тоді метод Зейделя для відомого вектора на к-ой ітерації має вигляд попередньої ітерації. Так само як і в методі простих ітерацій будується еквівалентна СЛАР і за початкове наближення береться вектор правих частин . Тоді метод Зейделя для відомого вектора на к-ой ітерації має вигляд:

З цієї системи видно, що, де В - нижня трикутна матриця з діагональними елементами, рівними нулю, а З - верхня трикутна матриця з діагональними елементами, відмінними від нуля, б=B+C. Отже

При такому способі приведення результатної СЛАР до еквівалентного вигляду метод простих ітерацій носить назву методу Якобі.

звідки

Таким чином, метод Зейделя є методом простих ітерацій з матрицею правих частин б=(E-B)-1C і вектором правих частин (E-B)-1в.

Вирішимо систему відносно невідомих при ненульових діагональних елементах aii_0, i= 1,n (якщо який-небудь коефіцієнт на головній діагоналі рівний нулю, достатньо відповідне рівняння поміняти місцями з будь-яким іншим рівнянням). Одержимо наступні вирази для компонентів вектора в і матриці б еквівалентної системи

Як нульове наближення x(0) вектора невідомих приймемо вектор правих частин x(0)=в. Тоді метод простих ітерацій прийме вигляд:

З (1.19) видно перевага ітераційних методів по порівнянню, наприклад, з розглянутим вище методом Гауса. У обчислювальному процесі беруть участь тільки твори матриці на вектор, що дозволяє працювати тільки з ненульовими елементами матриці, значно спрощуючи процес зберігання і обробки матриць.

Має місце наступна достатня умова збіжності методу простих ітерацій.

Метод простих ітерацій (1.19) сходиться до єдиного рішення СЛАР при будь-якому початковому наближенні x(0), якщо яка-небудь норма матриці б еквівалентної системи менше одиниці

Якщо використовується метод Якобі (вирази (1.18) для еквівалентної СЛАР), то

достатньою умовою збіжності є діагональне переважання матриці A, тобто

(для кожного рядка матриці A модулі елементів, що стоять на головній діагоналі, більше суми модулів недіагональних елементів). Очевидно, що в цьому випадку ||б||_c менше одиниці і, отже, ітераційний процес (1.19) сходиться.

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

При виконанні достатньої умови збіжності оцінка погрішності рішення на k- ой ітерації дається виразом:

Слід підкреслити, що ця нерівність дає завищене число ітерацій до, тому рідко використовується на практиці

1.4 Чисельний метод рішення задачі

Вирішимо систему відносно невідомих при ненульових діагональних елементах aiiV0, i= 1,n (якщо який-небудь коефіцієнт на головній діагоналі рівний нулю, достатньо відповідне рівняння поміняти місцями з будь-яким іншим рівнянням). Одержимо наступні вирази для компонентів вектора в і матриці б еквівалентної системи:

При такому способі приведення результатної СЛАР до еквівалентного вигляду метод простих ітерацій носить назву методу Якобі.

Як нульове наближення x(0) вектора невідомих приймемо вектор правих частин x(0)=в. Тоді метод простих ітерацій прийме вигляд:

З (1.19) видно перевага ітераційних методів по порівнянню, наприклад, з розглянутим вище методом Гауса. У обчислювальному процесі беруть участь тільки твори матриці на вектор, що дозволяє працювати тільки з ненульовими елементами матриці, значно спрощуючи процес зберігання і обробки матриць.

Має місце наступна достатня умова збіжності методу простих ітерацій.

Метод простих ітерацій (1.19) сходиться до єдиного рішення СЛАР при будь-якому початковому наближенні x(0), якщо яка-небудь норма матриці б еквівалентної системи менше одиниці

Якщо використовується метод Якобі (вирази (1.18) для еквівалентної СЛАР), то

достатньою умовою збіжності є діагональне переважання матриці A, тобто

(для кожного рядка матриці A модулі елементів, що стоять на головній діагоналі, більше суми модулів недіагональних елементів). Очевидно, що в цьому випадку ||б||_c менше одиниці і, отже, ітераційний процес (1.19) сходиться.

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

При виконанні достатньої умови збіжності оцінка погрішності рішення на k- ой ітерації дається виразом:

де x·- точне рішення СЛАР.

Процес ітерацій зупиняється при виконанні умови, де ее_)(kе - точність, що задається обчислювачем.

Беручи до уваги, що з (1.20) слідує нерівність, можна одержати апріорну оцінку необхідного для досягнення заданої точності числа ітерацій. При використанні як початкове наближення вектора в така оцінка визначиться нерівністю:

звідки одержуємо апріорну оцінку числа ітерацій до при ||б||<1

Слід підкреслити, що ця нерівність дає завищене число ітерацій до, тому рідко використовується на практиці.

1.6 Текст програми

program Yakobi;

uses crt;

const

maxn=100;

type

matrix=array[1..maxn,1..maxn] of real;

vector=array[1..maxn] of real;

vector1=array[1..maxn] of real;

var

i,j,n,k1: integer;

e,norma:real;

а: matrix;

b: vector;

x2,x3: vector1;

imya,dannble_i_rezultat,ekran:string;

procedure input(var kolvo:integer; var pogreshnostb:real; var matr1:matrix; var matr2:vector);

{ Введення початкових даних}

var i,j,code:integer;

a_s: string;

b_s: string;

begin

writeln('введіть кількість рівнянь');

readln(kolvo);

writeln('введіть погрішність обчислення');

readln(pogreshnostb);

writeln('введіть матрицю коефіцієнтів при невідомих');

for i:=1 to kolvo do

for j:=1 to kolvo do

begin

repeat

begin

writeln('введіть елемент [',i,',',j,'] і натисніть Enter');

readln(a_s);

if (i=j) and (a_s='0') then

repeat

begin

writeln('Елементи на головній діагоналі повинні бути відмінні від нуля. Повторіть введення');

readln(a_s);

end;

until a_s<>'0';

val(a_s,matr1[i,j],code);

end;

until code=0;

end;

writeln('введіть вектор вільних коефіцієнтів');

for j:=1 to kolvo do

begin

repeat

writeln('введіть елемент ',j,' і натисніть Enter');

readln(b_s);

val(b_s,matr2[j],code);

until code=0;

end;

end;

1.7 Тестовий приклад

2. Поліноміальная інтерполяція функції методом Ньютона з

розділеними різницями

2.1 Постановка завдання

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

2.2 Математичне формулювання завдання

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

i¨xi¨yi¨1. _

0	x0	у
1	x1	0
2	x2	у
...	...	1
n	xn	у
		2
		...
		у
		n

або

Крапки з координатами (xi, yi) називаються вузловими точками або вузлами.

Кількість вузлів в табличній функції рівна N=n+1.

Необхідно знайти значення цієї функції в проміжній крапці, наприклад x=D, причому .

Для вирішення завдання будуємо інтерполяційний многочлен.

2.3 Огляд існуючих чисельних методів рішення задачі

Інтерполяція по Лагранжу

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

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

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

якщо x=x₀, то Ln(x₀)= y₀,

якщо x=x₁, то Ln(x₁)= y₁,

.....

якщо x=x_n, то Ln(xn)= yn.

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

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

Чим більше вузлів інтерполяції на відрізку [x₀,x_n], тим точніше інтерполяційний многочлен наближає задану табличну функцію, тобто тим точніше рівність:

Проте із збільшенням числа вузлів інтерполяції зростає ступінь інтерполяційного многочлена n і в результаті значно зростає об'єм обчислювальної роботи. Тому при великому числі вузлів необхідно застосовувати ЕОМ. В цьому випадку зручно знаходити значення функції в проміжних крапках, не одержуючи многочлен в явному вигляді.

При рішенні задачі екстраполювання функції за допомогою інтерполяційного многочлена обчислення значення функції за межами відрізка [x₀,x_n] звичайно проводять не далі, чим на один крок h, рівний найменшій величині

оскільки за межами відрізка [x₀,x_n] погрішності, як правило, збільшуються.

Інтерполяція по Ньютону

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

(2.2)

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

- розділені різниці 0-го, 1-го, 2-го.., n-го порядку, відповідно.

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

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

Використовуючи теорію вигину бруса при малих деформаціях, можна показати, що сплайн - це група кубічних многочленів, в місцях сполучення яких перша і друга похідні безперервні. Такі функції називаються кубічними сплайнами. Для їх побудови необхідно задати коефіцієнти, які єдиним чином визначають многочлен в проміжку між даними крапками.

Наприклад, для деяких функцій (рис.) необхідно задати всі кубічні функції q_{1(x), q2(x), .qn(x).}

У найбільш загальному випадку ці многочлени мають вигляд:

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

Перші 2n умов вимагають, щоб сплайни стикалися в заданих точках:

Наступні (2n-2) умов вимагають, щоб в місцях зіткнення сплайнів були рівні перші і другі похідні:

Система рівнянь, алгебри, має рішення, якщо число рівнянь відповідає числу невідомих. Для цього необхідно ввести ще два рівняння. Звичайно використовуються наступні умови:

При побудові алгоритму методу перші і другі похідні зручно апроксимувати розділеними різницями відповідних порядків.

Одержаний таким чином сплайн називається природним кубічним сплайном. Знайшовши коефіцієнти сплайна, використовують цю шматково-гладку полиноминальную функцію для представлення даних при інтерполяції.

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

Значення f(x₀), f(x₁) ., f(xn), тобто значення табличної функції у вузлах, називаються розділеними різницями нульового порядку (k=0).

Відношення називається розділеною різницею першого порядку (k=1) на ділянці [x₀, x₁] і рівно різниці розділених різниць нульового порядку на кінцях ділянки [x₀, x₁], розділеній на довжину цієї ділянки.

Для довільної ділянки [xi x_i+1] розділена різниця першого порядку (k=1) рівна

Відношення називається розділеною різницею другого порядку (k=2) на ділянці [x₀, x₂] і рівно різниці розділених різниць першого порядку, розділеній на довжину ділянки [x₀, x₂].

Для довільної ділянки [xi x_i+2] розділена різниця другого порядку (k=2) рівна

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

(2.3)

Де n - ступінь многочлена.

Максимальне значення до рівне n. Тоді i =0 і розділена різниця n-го порядку на ділянці [x₀,x_n] рівна

,

тобто рівна різниці розділених різниць (n-1)-го порядку, розділеній на довжину ділянки [x₀,x_n].

Розділені різниці

є цілком певними числами, тому вираз (2.2) дійсно є многочленом, алгебри n-й ступеня. При цьому в многочлені (2.2) всі розділені многочленом, алгебри, n-й ступені. При цьому в многочлені (2.2) всі розділені різниці визначені для ділянок [x₀, x_0+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 Текст програми

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:\tp7\bp\bgi');

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);

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:\tp7\bp\bgi');

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

2.6 Тестовий приклад

3. Середньоквадратичне наближення функції

3.1 Постановка завдання

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

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

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

(3.1)

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

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

3.3 Огляд існуючих чисельних методів рішення задачі

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

Коли рівень невизначеності в завданні функції f(xi), що наближається i=1..m, достатньо великий, що характерне для обробки експериментальних даних, безглуздо вимагати виконання умов інтерполяції; крім того, число точок завдання функції f(xi) часто вельми велике. Все це робить застосування інтерполяції мало перспективним з причин поганої обумовленості завдання високої розмірності і проблем збіжності процесу інтерполяції

Однієї з найбільш простих і, тому, широко використовуваних функцій, що наближають, є поліном, алгебри

Метод середньоквадратичного наближення забезпечує побудову полінома Pn(x), виходячи з мінімізації величини

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

у просторі параметрів a₀, a₁,...,an. Існують різні підходи до рішення задачі мінімізації функції D(a). Простий з них приводить до необхідності рішення нормальної системи лінійних рівнянь алгебри

Проте, вже при n > 5 матриця такої системи виявляється настільки погано обумовленою, що одержані з (3.4) значення aj виявляються мало придатними для обчислення Pn(x). Тому, при необхідності побудови поліномів якнайкращого середньоквадратичного наближення вищих ступенів застосовують інші алгоритми, наприклад, метод сингулярного розкладання.

3.4 Чисельний метод рішення задачі

Можна розглянути два завдання:

1 - підібрати функцію так, щоб виконувалася нерівність

; (3.5)

2 - знайти якнайкраще наближення, тобто таку функцію, щоб було справедливим співвідношення

. (3.6)

Далі займемося відшуканням якнайкращого наближення.

Розкладемо функцію за системою лінійно незалежних функцій :

. (3.7)

Надалі для скорочення запису користуватимемося визначенням скалярного твору в просторі функцій :

.

Підставляючи (3.7) в умову (3.6), одержимо

 .

Диференціюючи цей вираз по і прирівнюючи похідні нулю, одержимо

. (3.8)

Визначник цієї системи є визначник Грама функцій . Через їх лінійну незалежність цей визначник не рівний нулю. Отже, з системи (3.8) можна знайти коефіцієнти, що визначають функцію згідно (3.6) і що мінімізують інтеграл від погрішності . Таким чином, якнайкраще середньоквадратичне наближення існує і воно єдино.

При використанні ортонормованої системи функцій система (3.8) спрощується:

,

тобто є коефіцієнтами Фурьє, а якнайкраще наближення є ряд Фурьє, що обривається на якомусь члені.

Доведено, що в будь-якому лінійно нормованому просторі при лінійній апроксимації вигляду (3.4) якнайкраще наближення існує, хоча воно може бути не єдиним.

У тих випадках, коли функції не ортогональні, при визначник Грама зменшується, наближаючись до нуля. Тоді система стає погано обумовленою і її рішення дає велику погрішність. У цій ситуації звичайно беруть не більше п'яти-шести членів в сумі (3.7).

Як найчастіше використовують поліноми Лежандра, Чебишева, Лагерра, Ермита, ортогональні із заданою вагою.

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

(3.9)

де - число заданих вузлів.

Умова якнайкращого середньоквадратичного наближення записується таким чином:

. (3.10)

Вважаючи, де, і підставляючи цей многочлен в (3.10), дійдемо системи (3.8), в якій скалярні твори обчислюють згідно (3.9). Описана процедура апроксимації носить назву методу найменших квадратів.

Найбільш споживаний варіант методу найменших квадратів відповідає випадку статечного виду функцій, тобто, причому .

Система рівнянь (3.8) при цьому приймає вигляд

(3.11)

де .

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

Мал. 3.1 Основна програма



Мал. 3.2 Процедура введення даних



Рис 3.3 Процедура середньоквадратичного наближення

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

program srpribl; {$R+}

uses graph;

label 1,2,3,4;

const m=2;

type mas= array [1..21] of real;

mas1= array [1..m] of real;

mas2= array [1..m,1..m+1] of real;

var i,j:byte;

y1,x1:mas;

xx1:mas1;

a1:mas2;

procedure vvod (x,y: mas);

begin

x[1]:=0;y[1]:=0.290234387293458; x[2]:=0.25;y[2]:=0.201247759722173;

x[3]:=0.5;y[3]:=0.0712686786428094;x[4]:=0.75; у[4]:=0.044294935464859;

x[5]:=1;y[5]:=-0.0745576142333448; x[6]:=1.25;y[6]:=-0.0807833694852889;

x[7]:=1.5;y[7]:=-0.233371530473232;x[8]:=1.75;y[8]:=-0.283957027737051;

x[9]:=2;y[9]:=-0.308697660081089;x[10]:=2.25;y[10]:=-0.42091366359964;

x[11]:=2.5;y[11]:=-0.516991161741316;x[12]:=2.75;y[12]:=-0.427710095947851;

x[13]:=3;y[13]:=-0.374748698528856;x[14]:=3.25;y[14]:=-0.229905794281513;

x[15]:=3.5;y[15]:=-0.205365492492496639;x[16]:=3.75;y[16]:=-0.129155068378896;

x[17]:=4;y[17]:=-0.0438349825330079;x[18]:=4.25;y[18]:=0.0294586319476366;

x[19]:=4.5;y[19]:=0.132592592108995;x[20]:=4.75;y[20]:=0.206369574274868;

x[21]:=5;y[21]:=0.26959548862651;

end;

procedure srpribl (xx:mas1;a:mas2);

var l:real;

s,z,k1:integer;

begin

for s:=1 to m-1 do

for z:=s+1 to m do

begin

а[z,s]:=-a[z,s]/a[s,s];

for k1:=s+1 to m+1 do а[z,k1]:=a[z,k1]+a[z,s]*a[s,k1]

end;

xx[m]:=a[m,m+1]/a[m,m]; writeln(' xx[',m,']=',xx[m]:2:3);

for i:=m-1 downto 1 do

begin

l:=a[i,m+1];

for j:=i+1 to m do l:=l-xx[j]*a[i,j];

xx[i]:=l/a[i,i]; writeln(' xx[',i,']=',xx[i]:2:3)

end

end;

procedure Grafik (var x,y:mas;xx:mas1)

var ec,gd,gm:integer;

begin

gd:=detect;

initgraph (gd,gm,' ');

ec:=graphresult;

if ec<>grok then begin

writeln ('Oshibka v inicializacii grafika');

halt (1);

end;

setcolor(white);

line (0,420,620,420);

line (0,0,0,420);

setcolor (white);

for i:=1 to 21 do begin

circle (round (x[i]*150),round (у[i]*100),1);

end;

setcolor (yellow);

for i:=1 to m do begin

circle (round (x[i]*150),round (xx[i]*100),1);

end;

setcolor (green);

for i:=2 to m do begin

line (round (x[i-1]*150),round(xx[i-1]*100),round (x[i]*150)

round (xx[i]*100));

end;

end;

begin

vvod(x1,y1);

for i:=1 to 2 do

for j:=1 to 3 do а[i,j]:=0;

а[1,1]:=21;

for i:=1 to 21 do

begin

а[1,2]:=a[1,2]+x[i];

а[2,1]:=a[2,1]+x[i];

а[2,2]:=a[2,2]+x[i]*x[i];

а[1,3]:=a[1,3]+y[i];

а[2,3]:=a[2,3]+y[i]*x[i]

end;

srpribl(xx1,a1);

for i:=1 to 21 do

writeln(у[i]:2:3,' ',xx[1]+xx[2]*x[i]:2:3);

Grafik(x1,y1,xx1);

end.

3.7 Тестовий приклад

Висновок

У даній роботі описані і реалізовані за допомогою блок-схем і мови програмування Turbo Pascal базові завдання обчислювальної математики: рішення систем лінійних рівнянь алгебри, поліномінальна інтерполяція, середньоквадратичне наближення функції, чисельна інтеграція функцій. Представлені методи і реалізовані алгоритми достатньо прості, але в теж час ефективні для великої кількості завдань.

Список використаної літератури

1. Хвальків Н., Жідков Н., Псів Г. Численниє методи. М.: Лабораторія базових знань, 2001. 632 з.

2. Вержбіцкий В.М., Чисельних методи. Лінійна алгебра і нелінійні рівняння. М.: Вища школа, 2000. 266 з.

3. Вержбіцкий В.М., Основи чисельних методів. М.: Вища школа, 2002. 840 з.

4. Пірумов У.Г. Чисельні методи . М.: Дрохва, 2003. 224 з.

5. Буслов В.А., Яковльов С.Л. Чїсленниє методи і рішення рівнянь. Санкт-Петербург, 2001. 44 з.

6. Фаронов В.В. Турбо Паскаль 7.0. Початковий курс. Навчальний посібник. - М.: Нолідж, 1997. - 616 з.