МІНІСТЕРСТВО ОСВІТИ І НАУКИ УКРАЇНИ
СУМСЬКИЙ ДЕРЖАВНИЙ УНІВЕРСИТЕТ
КАФЕДРА ІНФОРМАТИКИ
До У Р З Про У А Я Р А Б Про Т А
З чисельних методів
на тему:
Рішення диференціального рівняння
ПЕРШОГО ПОРЯДКУ
Суми, 2005 р.
1. Метод Адамса
Цей метод чисельного інтегрування розроблений Адамсом в 1855р. Надалі цей метод був забутий і знову відкритий на початку століття. Популяризація методу Адамса і подальше його удосконалення пов'язані з ім'ям А.Н. Крилова.
Викладемо метод Адамса стосовно до рівняння першого порядку
(1)
з початковою умовою
(2).
Нехай x (I = 0,1,2, ....) - Система равностоящими значень з кроком h і
=
. Очевидно, маємо
(3).
В чинності другий інтерполяційної формули Ньютона з точністю до різниць четвертого порядку отримуємо
(4)
де .
Підставляючи вираз (4) в формулу (3) і враховуючи, що dx = hdq, матимемо
Звідси отримуємо екстраполяційних формулу Адамса
. (5)
Для початку процесу потрібні чотири початкові значення , Так званий початковий відрізок, який визначають виходячи з початкової умови (2), яким-небудь чисельним методом. Можна, наприклад, використовувати метод Рунге-Кутта. Знаючи ці значення, з рівняння (1) можна знайти значення похідних
і скласти таблицю різниць.
(6)
Подальші значення (I = 4,5, ...) шуканого рішення можна крок за кроком обчислювати за формулою Адамса, поп про лняя в міру необхідності таблицю різниць (6).
Для контролю рекомендується обчисливши перше наближення для за формулою
визначити , Підрахувати кінцеві різниці.
,
,
(7)
і потім знайти друге наближення за більш точною формулою
(8)
Якщо і
відрізняються лише на кілька од і ниць останнього зберігається десяткового розряду, то можна покласти
, А потім, знайшовши
, Перевичіслів кінцеві різниці (7). Після цього, строго кажучи, слід знову знайти
за формулою (8). Тому крок h повинен бути таким, щоб цей перерахунок був зайвим.
На практиці крок h вибирають настільки малим, щоб можна було знехтувати членом у формулі (8).
Якщо ж розбіжність величин і
значно, то слід зменшити крок h.
Зазвичай крок h зменшують в два рази. Покажемо, як у цьому випадку, маючи до деякого значення i таблицю величин і
,
(J
i) c кроком
, Можна просто побудувати таблицю величин
(K = 0,1,2 ...) з кроком
. Для кратності введення скорочені позначення:
(K = 0,1,2 ...).
На основі формули (4) будемо мати
, (9)
де . Звідси, вважаючи j = i-2 і q = 1 / 2 та враховуючи, що
, Знаходимо
. (10)
Аналогічно при j = i-1, q = 1 / 2 з формули (9) отримуємо, що аргументу відповідає значення
. (11)
Що стосується значень і
, То вони є в старій таблиці. Після цього складаємо початковий відрізок для нової таблиці.
і знаходимо кінцеві різниці:
(K =- 3, -2, -1),
(K =- 3, -2),
(K =- 3,).
Далі таблиця триває звичайним шляхом, за допомогою відповідної модифікації формули (5):
,
(J = 0,1,2, ...).
Для роботи на електронних лічильниках машинах формулу Адамса (5) вигідно застосовувати в розкритому вигляді. Враховуючи, що
після приведення подібних членів маємо
,
причому .
2. Методи, засновані на застосуванні похідних вищих порядків
До цих пір для чисельного інтегрування диференціального рівняння першого порядку
(1)
з початковою умовою
(2)
ми застосовували формули, в яких явно використовується лише перша похідна шуканого рішення.
Проте якщо використовувати формули, явно містять похідні вищих порядків від шуканого рішення, то можна вказати методи, що дають більш точний результат на даному проміжку без збільшення числа кроків.
Виведемо відповідні формули, припускаючи, що права частина рівняння (1) дифференцируема достат про чное число разів.
Нехай - Значення шуканого рішення y = y (x) і, відповідно, значення його похідних першого і другого порядків в точках
. Маючи величини
в ряди за степенями h, знаходимо:
З отриманих формул виключимо члени, що містять і
.
Для цього другу формулу помножимо на , А третю - на
і складемо з першою. Будемо мати:
Таким чином, з точністю до маємо наближену формулу
(3)
Можна показати, що залишковий член формули (3) дорівнює де
Аналогічно маємо:
і
Звідси
З іншого боку
Тому
Таким чином, з точністю до h 5 маємо прибл і женную формулу
(4)
Можна довести, що залишковий член формули (4) є
де
До формулами (3) і (4) прис е д і ним вирази для похідних:
(5)
(6)
Процес чисельного диференціювання рівняння (1) при наявності початкової умови (2), іспользющій формули (3) і (4), відбувається наступним чином. Яким-яким методом обчислюємо три початкові рядки (початкова таблиця):
З формули (4) при i = 2 отримуємо перше наближення для :
(7)
і, користуючись формулами (5) і (6), знаходимо для відповідних похідних і
їх перші наближення:
і
.
Друге наближення для визначаємо при i = 2 з формули (3):
(8)
Після цього виправляємо значення похідних і
, Під рах ітивая їх другі наближення:
і
.
Для контролю ще раз обчислюємо за формулою (3) третє наближення значення
, Використовуючи знайдені значення
і
.
Якщо крок h обраний підходящою кандидатурою, то перещет не дає нового результату, і в цьому випадку можна покласти
В іншому випадку слід зменшити крок. Аналогічно знаходяться подальші значення при i> 3.
Для отримання початкових значень і
зазвичай використовують метод послідовних наближень або метод Рунге-Кутта, після чого потрібні похідні
і
(I = 0,1,2) визначаються за формулами (5) і (6).
Можна також застосувати наступний прийом: спочатку, використовуючи дане початкове значення , Безпосередньо обчислюємо
і
.
Тим самим буде заповнена перший рядок початковій таблиць и.
Далі на підставі формули Тейлером наближено отримуємо
і, отже, можна будить знайти
і
.
Користуючись цими даними, уточнюємо значення за формулою (3):
і потім Переобчислювати значення і
. Тим самим заповнюємо другий рядок початковій таблиці. Аналогічно, виходячи з другого рядка, знаходимо елементи
,
і
останньої, третьої рядки початкової таблиці.
Відзначимо, що якщо перерахунки елементів рядків дають значні розбіжності, то цей прийом не є надійним. У такому випадку слід або зменшити крок h обчислень, або ж звернутися до більш точним методам.
На закінчення наведемо формули, що забезпечують більш високу ступінь точності, але вимагають обчислення, крім другої, ще й третьої похідної шуканого рішення. А саме, використовуючи Формулу Тейлером і вживаючи прийом, аналогічний вказаному вище, отримуємо формули
, (11)
де
, І
, (12)
де .
Формула (11) вживається для знаходження першого наближення ; Формула (12) дає уточнене значення
. Само собою зрозуміло, що до останніх двох формулах доцільно вдаватися тоді, коли форма диференціального рівняння дозволяє порівняно просто знаходити другу і третю похідні від шуканої функції y.
Додаток
program proizw_w_p;
uses crt;
const epsilon = 0.05;
type mas = array [1 .. 100] of real;
nabl = array [1 .. 3] of real;
var i: integer;
x, y, y1, y2: mas;
nabl1, nabl2, nabl3: nabl;
a, h: real;
n: integer;
function f (x, y: real): real;
begin
f: = sqr (x)-sqr (y);
end;
procedure metod (xi, yi, step: real; var rez: real);
var k1, k2, k3, k4: real;
begin
k1: = F (xi, yi);
k2: = F (xi + step / 2, yi + k1 * step / 2);
k3: = F (xi + step / 2, yi + k2 * step / 2);
k4: = F (xi + step, yi + k3 * step);
rez: = yi + (step / 6) * (k1 +2 * k2 +2 * k3 + k4)
end;
procedure osn_metod (xi, yi, step: real; var yh22: real; var h: real);
var yh, yh2: real;
begin
repeat
metod (xi, yi, step, yh);
metod (xi, yi, step / 2, yh2);
metod (xi, yh2, step / 2, yh22);
if abs (yh-yh22) / 15> epsilon then step: = h / 2;
until abs (yh-yh22) / 15 <epsilon;
end;
procedure iteraziya (j: integer; xi, h: real);
begin
{Перше наближення}
nabl1 [1]: = y [j-3] +3 * (y [j-1]-y [j-2]) + sqr (h) * y2 [j-1]-y2 [j-2];
{Похідна першого наближення}
nabl1 [2]: = sqr (xi)-sqr (nabl1 [1]);
{Друга похідна першого наближення}
nabl1 [3]: = 2 * (xi-nabl1 [1] * nabl1 [2]);
{Друге наближення}
nabl2 [1]: = y [j-1] + (h / 2) * (y1 [j-1] + nabl1 [2 ])+(( sqr (h)) / 12) * (nabl1 [3] - y2 [j-1]);
{Похідних другого наближення}
nabl2 [2]: = sqr (xi)-sqr (nabl2 [1]);
{Друга похідна Другого наближення}
nabl2 [3]: = 2 * (xi-nabl2 [1] * nabl2 [2]);
{Третього наближення}
nabl3 [1]: = y [j-1] + (h / 2) * (y1 [j-1] + nabl2 [2]) - (sqr (h) / 12) * (nabl2 [3]-y2 [ j-1]);
{Похідна третього наближення}
nabl3 [2]: = sqr (xi)-sqr (nabl3 [1]);
{Друга похідна третє наближення}
nabl3 [3]: = 2 * (xi-nabl2 [1] * nabl2 [2]);
end;
procedure solution (h: real);
begin
{============== Метод Рунге-Кута =============================== ==}
a: = 0;
i: = 1;
y [1]: = 1;
while i <4 do
begin
x [i +1]: = a + i * h;
osn_metod (x [i], y [i], h, y [i +1], h);
inc (i);
end;
{====== Кінець методу Рунге-Кута =================================}
{============ Знайдемо перші і другі похідні ===============}
for i: = 1 to 3 do
begin
y1 [i]: = sqr (x [i])-sqr (y [i]);
y2 [i]: = 2 * (x [i]-y [i] * y1 [i]);
end;
{================================================= ================}
for i: = 4 to n do
begin
iteraziya (i, x [i], h);
if abs (nabl3 [1]-nabl2 [1]) <epsilon
then
begin
y [i]: = nabl3 [1];
y1 [i]: = nabl3 [2];
y2 [i]: = nabl3 [3];
end
else
begin
h: = h / 2;
if keypressed then halt;
solution (h);
end;
end;
end;
BEGIN
{===================== Init =========================== ===============}
clrscr;
write ('введіть кількість значень, які необхідно обчислити n =');
readln (n);
h: = 0.1;
{================== End of init ============================ =============}
for i: = 1 to n do
begin
x [i]: = (i-1) * h;
end;
solution (h);
for i: = 1 to n do
begin
write ('y [', i, '] =', y [i], 'y "[', i, '] =', y1 [i], 'y ""[', i,'] = ' , y2 [i]);
writeln;
end;
writeln ('');
writeln ('');
write ('Press <enter> to exit ....');
readln;
END.