Курсова робота з інформатики
Тема: «Чисельні методи при вирішенні завдань»
Автор: студент групи ПС-146
Нечаєв Л. В.
Перевірив: Альошин Є. А.

Зміст

"1-3" Зміст. 2
Програми та опису. 3
Програма для вирішення завдання 17. 3
Умова задачі 17. 3
Рішення задачі за методом Адамса. 3
Блок-схема функції main з програми 17.c. 4
Блок-схема функції Adams з програми 17.c. 5
Лістинг програми 17.c. 6
Результат розв'язання задачі 17 на ЕОМ .. 9
Висновок: 9
Програма для вирішення завдання 30. 10
Умова задачі 30. 10
Рішення задачі за методом найменших квадратів. 10
Блок-схема функції main з програми 30.c. 11
Блок-схема функції MMinor з програми 30.c. 11
Блок-схема функції MatrixMultiply з програми 30.c. 12
Блок-схема функції Determinant з програми 30.c. 12
Лістинг програми 30.c. 12
Результат розв'язання задачі 30 на ЕОМ .. 17
Висновок: 17

Програми та опису

Програма для вирішення завдання 17

Умова задачі 17.

Розробити функцію чисельного інтегрування системи диференціальних рівнянь методом Адамса. Прототип функції:
void Adams (
void f (double * y, double * ys, double t),
double * y,
int n,
double tn,
double tk,
int m,
double eps);
де:
f - Функція обчислення правих частин системи диференціальних рівнянь:
y - Масив розміру n значень залежних змінних;
ys - Масив розміру n значень залежних похідних;
n - Порядок системи диференціальних рівнянь;
t - Незалежна змінна;
tn - Початкове значення інтервалу інтегрування;
tk - Кінцеве значення інтервалу інтегрування;
m - Початкове число розбиття відрізка інтегрування [tn; tk]
eps - відносна похибка інтегрування. Обчислення припиняються, коли

, Де

- Значення i-ї компоненти вектора залежних змінних при t = tk для кількості розбиття відрізка інтегрування m.
Початкові кроки робляться за методом Рунге-Кутта.
Застосувати цю функцію для інтегрування диференціального рівняння 3-його порядку y ⁽³⁾ +2 y''+3 y '+ y = 5 + x ² в інтервалі xÎ [0; 2] з кроком Δ x = 0, і початковими умовами x = 0; y (0) = 1; y '(0) = 0.1; y''(0) = 0.

Рішення задачі за методом Адамса

Для запуску екстраполяціонного методу Адамса потрібно 4 початкових значення функції. Одне значення вже поставлено, а інші виходять за методом Рунге-Кутта 4 порядку. Після обчислення значення в кінці відрізка відбувається обчислення відносної похибки (з поточних і раніше отриманих з кроком h значень функції) і порівняння її із заданим значенням. Якщо отримана погрішність менше, ніж задана, то вважається, що завдання виконано, і відбувається повернення в зухвалу програму з отриманим значенням функції. Якщо ж ні - то зменшується в 2 рази крок і весь процес, починаючи з методу Рунге-Кутта, повторюється знову (для обчислення нових значень функції). Так продовжується до тих пір, поки отримане значення похибки не стане менше ніж задане.
Для роботи програми необхідна функція обчислення правих частин системи диференціальних рівнянь. Це функція func (double * y, double * ys, double x). Так як в задачі потрібно вирішити рівняння y ⁽³⁾ +2 y''+3 y '+ y = 5 + x ^2, складаємо систему диференціальних рівнянь першого порядку. Виглядає вона так:

При кожному обчисленні лівих частин цієї системи відбувається диференціювання y, y 'і y'', тобто обчислення відповідно нових значень y ', y'', y'''.
Ну, а якщо перекласти це все в програму на Сі, то вийде функція func (дивися лістинг 17 завдання).

Блок-схема функції main з програми 17.c

SHAPE \ * MERGEFORMAT

Виклик функції Adams
Adams (func, y, 3, xs, xe, 100, 0.001);

xs + = 0.1; xe + = 0.1;

Висновок значення залежної змінної і результату інтегрування

I + +

Початок

Завдання початкових значень:
y '= 1; y''= 0.1; y''' = 0;
xs = 0; xe = 0.1;

I = 0

I <20?

Так

Ні

Кінець

Блок-схема функції Adams з програми 17.c

SHAPE \ * MERGEFORMAT

Ні

Обчислення наступного значення функції YA = (Q3 / h) + (y3 + (1. / 2.) * Dq2 + (5./12.) * D2q1 + (3. / 8.) * D3q0);

Виклик функції Adams

Перевірка аргументів функції на правильність

Виділення пам'яті для всіх масивів

Не вдалося

Вдалося

Завершення програми

Розподіл пам'яті між всіма масивами

Ініціалізація змінних і масивів q0, xi початковими значеннями

i = 0

i <3?

Ні

Так

Обчислення K1 = f (xi, Y)

Корекція K1 для кожного рівняння
(K1 = K1 * h)

Обчислення аргументів для наступної функції (YA = Y + K1 / 2.)

Обчислення K2 = f (xi + h / 2., YA)

Корекція K2 для кожного рівняння
(K2 = K2 * h)

Обчислення аргументів для наступної функції (YA = Y + K2 / 2.)

Обчислення K3 = f (xi + h / 2., YA)

Корекція K3 для кожного рівняння
(K3 = K3 * h)

Обчислення аргументів для наступної функції (YA = Y + 3 * K3 / 2.)

Обчислення K4 = f (xi + h, YA)

Корекція K4 для кожного рівняння
(K4 = K4 * h)

Обчислення наступного значення функції (q [i +1] =
q [i] + (1. / 6.) * (K1 +2 * K2 +2 * K3 + K4))

xi + = h; i + +

Корекція: q [i] = q [i] * h

Обчислення
Δq _2, Δq _1, Δq _0, Δ ² q _1, Δ ² q _0, Δ ³ q0

x + = h

xi <tk? (Відрізок закінчився?

Так

Є з чим порівнювати значення функції в точці tk?

flag = 1

Так

Похибка менше заданої?

Так

Ні

Копіювати yt в y

Повернення в основну програму

Копіювати ya в yt

Зменшити крок у 2 рази (h = h / 2.)

Ні

xi = tn (з початку відрізка)

Лістинг програми 17. c

/ / Завдання 17. Чисельне інтегрування системи диференціальних рівнянь
/ / Методом Адамса. Програма розрахована на компіляцію в Micro $ oft C 6.00
/ / Або Borland C 3.1 +
/ / (C) 2004 REPNZ. All rights reserved. Release date is 2.04.2004
# Include <stdio.h>
# Include <stdlib.h>
# Include <math.h>
void func (double * y, double * ys, double t)
{/ / Функція обчислення правих частин рівнянь
ys [0] = y [1]; / / ys [1]-перша похідна; ys [2]-друга і т.д.
ys [1] = y [2]; / / t-незалежний аргумент
ys [2] = 5 + t * t - y [0] - 3. * Y [1] - 2. * Y [2];
}
void Adams (
void f (double * y, double * ys, double x),
/ / Функція вичіленія правих частин системи
double * y, / / Масив розміру n значень залежних змінних
int n, / / Масив розміру n значень похідних
double tn, / / Початок інтервалу інтегрування
double tk, / / Кінець інтервалу інтегрування
int m, / / Початкове число розбиття відрізка інтегрування
double eps) / / Відносна похибка інтегрування
{
double * k1, * k2, k3 *, * k4; / / Для методу Рунге-Кутта
double * q0, * q1, * q2, * q3; / / Значення похідних Для методу Адамса
double * ya; / / Тимчасовий масив
double * y0, * y1, y2 *, * y3; / / Значення функції для методу Адамса
double h; / / Крок інтегрування
double xi; / / Поточне значення незалежної змінної
double eps2; / / Для оцінки похибки
double dq2, dq1, dq0, d2q1, d2q0, d3q0; / / збільшення
int flag = 0; / / 0, поки йде перший прорахунок
int i, j; / / Індекси

if (m <4) m = 4; / / Мінімум 4 відрізки
if (tn> = tk)
{Printf ("\ nНеправільние аргументи \ n");
abort (); / / Неправильні аргументи
}
/ / Виділяємо пам'ять для масивів із змінними
if ((k1 = malloc ((4 + 4 + 4 + 1) * n * sizeof (double))) == 0)
{Printf ("\ nОшібка розподілу пам'яті \ n");
abort (); / / Перервати, якщо не вдалося
}
/ / Розподіляємо пам'ять між масивами:
/ / Для методу Рунге-Кутта 4 порядку
k2 = k1 + n; k3 = k2 + n; k4 = k3 + n;
/ / 4 пердидущіх значення функції
y0 = k4 + n; y1 = y0 + n; y2 = y1 + n; y3 = y2 + n;

/ / Для тимчасового масиву збору даних
ya = y3 + n;
/ / Для методу Адамса
q0 = ya + n; q1 = q0 + n; q2 = q1 + n; q3 = q2 + n;
h = (tk - tn) / m; / / Крок
eps = fabs (eps); / / Абсолютне значення похибки
start: / / Звідси починаються обчислення
xi = tn; / / Початок проміжку
/ / Обчислюємо значення функції y0 ... y3, тобто y [i-3] ... y [0]
/ / Перше значення системи рівнянь вже дано: y ...
////////////////////////////////////////////////// /////////////////////
/ / - Метод Рунге-Кутта 4 порядку - / /
////////////////////////////////////////////////// /////////////////////
for (j = 0; j <n; j + +) y0 [j] = y [j]; / / Копіюємо його в y0
f (y0, q0, xi); / / Заповнюємо q0, грунтуючись на значеннях з y0
for (j = 0; j <n; j + +) q0 [j] *= h; / / Робимо q0
xi + = h; / / Наступний крок
/ / ... а інші 3 видобуваємо за допомогою методу Рунге-Кутта 4 порядку.
for (i = 0; i <3; i + +) / / i - ЯКЕ ЗНАЧЕННЯ ВЖЕ Є
{/ / А Обчислюємо значення Y [i +1 ]!!!!
/ / Спочатку потрібні коефіцієнти k1
/ / Елемент y [i, j] = y0 + (i * n) + j = y0 [i * n + j]
f (& y0 [i * n], k1, xi); / / Обчислимо f (xi, yi) = k1 / h
/ / І для кожного диференціального рівняння системи проробляємо
/ / Операції обчислення k1, а також підготовки в ya аргументу для
/ / Обчислення k2
for (j = 0; j <n; j + +)
{
k1 [j] *= h; / / Обчислимо нарешті k1
ya [j] = y0 [i * n + j] + k1 [j] / 2.;
/ / І один з аргументів для функції
} / / Обчислення k2
f (ya, k2, xi + (h / 2.)); / / Обчислимо f (xi, yi) = k2 / h
for (j = 0; j <n; j + +)
{/ / Обчислимо нарешті k2
k2 [j] *= h;
ya [j] = y0 [i * n + j] + k2 [j] / 2.; / / І один з аргументів для функції
} / / Обчислення k3
f (ya, k3, xi + h / 2.); / / Обчислимо f (xi, yi) = k3 / h
for (j = 0; j <n; j + +)
{
k3 [j] *= h; / / Обчислимо нарешті k3
ya [j] = y0 [i * n + j] + k3 [j]; / / І один з аргументів для функції
} / / Обчислення k4
f (ya, k4, xi + h); / / Обчислимо f (xi, yi) = k4 / h
for (j = 0; j <n; j + +) k4 [j] *= h; / / Обчислимо нарешті k4
/ / Треба обчислити приріст кожної функції з n
for (j = 0; j <n; j + +) / / Обчислюємо таке значення
/ / Функції
/ / Y [i +1] = Yi + ...
y0 [(i +1) * n + j] = y0 [i * n + j] + (k1 [j] + 2. * k2 [j] + 2 * k3 [j] + k4 [j]) / 6 .;
/ / І нове значення q [i +1]
f (& y0 [(i +1) * n], & q0 [(i +1) * n], xi); / / qi = f (xi, yi);
for (j = 0; j <n; j + +) q0 [((i +1) * n) + j] *= h;
xi + = h; / / Наступний крок}
////////////////////////////////////////////////// /////////////////////
/ / - Метод Адамса - / /
////////////////////////////////////////////////// /////////////////////
/ / Отже, обчислені 4 перших значення. Цього достатньо для початку методу
/ / Адамса для кроку h.
/ / B y0 ... y3 лежать 4 значення функцій (_НЕ_ПРОІЗВОДНИХ!!!).
/ / A в q0 ... q3 лежать значення _проізводних_ цих функцій, помножених на h
/ / Q0 .. q3, а також y0 .. y3 представляють собою черги з 4 елементами
again: / / Обчислюємо нове значення функції Yi (Це Y [i +1])
for (j = 0; j <n; j + +)
{/ / Все збільшення
dq2 = q3 [j] - q2 [j]; dq1 = q2 [j] - q1 [j]; dq0 = q1 [j] - q0 [j];
d2q1 = dq2 - dq1; d2q0 = dq1 - dq0;
d3q0 = d2q1 - d2q0;
/ / Нове значення функції (в ya поки що)
ya [j] = y3 [j] + (q3 [j] + (dq2 / 2.) + (5. * d2q1 / 12.) + (3. * d3q0 / 8.));
/ / Зрушуємо всі масиви на 1 вперед і додаємо в чергу нове
/ / Значення функції
y0 [j] = y1 [j]; y1 [j] = y2 [j]; y2 [j] = y3 [j]; y3 [j] = ya [j];
/ / Просто зрушуємо q, нічого поки що не додаючи
q0 [j] = q1 [j]; q1 [j] = q2 [j]; q2 [j] = q3 [j];
}
/ / У чергу в якості q3 кладемо нове значення
f (y3, q3, xi); / / q3 = f (xi, y3);
for (j = 0; j <n; j + +) q3 [j] *= h; / / Обчислити q3
/ / Чергове значення функції обчислено. Следующиие крок
xi + = h;
/ / Продовжити інтегрування?
if (xi <tk) goto again; / / Так.
/ / Якщо перший раз тут, то прорахувати ще раз з кроком h / 2
if (flag == 0)
flag = 1; / / Порівнювати вже буде з чим
else
{
/ / Не перший раз - оцінити похибку
/ / Зараз у y3 - значення тільки що обчисленої функції,
/ / А в y2 - занчение функції, обчисленої з кроком h * 2
/ / По відношенню до поточного
for (j = 0; j <n; j + +)
{Eps2 = fabs (((y3 [j] - y2 [j]) / y2 [j]));
if (eps2> eps) break; / / Якщо похибка дуже велика
}
if (j == n) / / Якщо все ОК
{/ / Копіюємо результат
for (j = 0; j <n; j + +) y [j] = y3 [j];
free (k1); / / Звільняємо пам'ять
return; / / Повертаємося в main
}
}
/ / З якихось причин виходу з функції не сталося -
/ / Тоді зменшуємо крок в 2 рази і повторюємо
/ / Всі, починаючи з методу Рунге-Кутта
h / = 2.; / / Зменшити крок
goto start; / / Повторити розрахунок спочатку, з новими параметрами
}
int main ()
{
double y [3], xs, xe;
int i;
y [0] = 1.; y [1] = 0.1; y [2] = 0.; / / Початкові умови
xs = .0; xe = .1; / / Початок інтегрування
printf ("x =% 5.3lg, y (% 4.2lg) =% 10.3lg \ n", xs, xs, y [0]);
for (i = 0; i <20; i + +)
{
Adams (func, y, 3, xs, xe, 10, 1.e-3);
xs + = 0.1; xe + = 0.1;
printf ("x =% 5.3lg, y (% 4.2lg) =% 10.3lg \ n", xs, xs, y [0]);
}
return 0;
}

Результат розв'язання задачі 17 на ЕОМ

Для роботи програму необхідно скомпілювати в моделі не нижче SMALL. Використовувався компілятор Micro $ oft C 6.00 з однойменного пакета. Після запуску програма виводить наступне:
Програма чисельного інтегрування системи диференціальних
рівнянь екстраполяційних методом Адамса
Розробник: студент гр. ПС-146
Нечаєв Леонід Володимирович
17.03.2004
Диференціальне рівняння має вигляд y'''+ 2y''+ 3y' + y = x ^ 2 + 5
Отже, залежність y [x]:
x = 0, y (0) = 1
x = 0.1, y (0.1) = 1.01
x = 0.2, y (0.2) = 1.02
x = 0.3, y (0.3) = 1.04
x = 0.4, y (0.4) = 1.07
x = 0.5, y (0.5) = 1.11
x = 0.6, y (0.6) = 1.16
x = 0.7, y (0.7) = 1.22
x = 0.8, y (0.8) = 1.28
x = 0.9, y (0.9) = 1.37
x = 1, y (1) = 1.46
x = 1.1, y (1.1) = 1.56
x = 1.2, y (1.2) = 1.67
x = 1.3, y (1.3) = 1.79
x = 1.4, y (1.4) = 1.92
x = 1.5, y (1.5) = 2.06
x = 1.6, y (1.6) = 2.21
x = 1.7, y (1.7) = 2.36
x = 1.8, y (1.8) = 2.52
x = 1.9, y (1.9) = 2.69
x = 2, y (2) = 2.86

Висновок:

Перевіряємо рішення в програмі Mathematica 4.2. Результати, отримані з точністю до 2 знаків після коми не відрізняються від отриманих. Завдання вирішена вірно.

Програма для вирішення завдання 30.

Умова задачі 30.

Розробити програму апроксимації функції методом найменших квадратів для моделі

за таблицею результатів експерименту:

X ₁	X ₂	Y
1	1	0
-1	-1	-2
2	2	-2
3	-2	29
-2	4	54

Рішення задачі за методом найменших квадратів

Розраховується модель лінійна щодо своїх коефіцієнтів a _i. Задано матриці

, А також функція для отримання матриці F. F - Спеціальна матриця, яка обчислюється за алгоритмом, наведеним нижче. Функція являє собою мою власну розробку, але цілком можливо її вводити вручну. Алгоритм складання матриці F (враховуючи розкладання

, Де

- Функції з моделі y, а

.- N-й елемент матриці

.
Виходячи з цих формул будується функція f (дивися лістинг програми 30.c).
Далі, за формулою

знаходиться матриця з коефіцієнтами a _i і виводиться на екран.

Блок-схема функції main з програми 30.c

Ні

SHAPE \ * MERGEFORMAT

Dt дорівнює 0?

Кінець

Початок

Вдалося виділити пам'ять для всіх матриць F і FT?

Так

Завдання початкових значень в матриці F

Транспонування матриці F
(TF = F ^T)

Обчислення TMP = (F ^T * F)

Dt = Determinant (TMP)

Ні

Так

J = 0

j <3?

Ні

Так

i = 0

i <3?

Так

Ні

AC2 = M (TMP, a _ij)

REV (a _ji) = Determinant (AC2) / Dt * ((i + j)% 2 == 0)? 1: -1

I + +

F = (FT * F) -1 * F

REV = (F ^T * F) ^-1 * F * Y

Висновок одержані коефіцієнтів

j + +

Блок-схема функції MMinor з програми 30.c

SHAPE \ * MERGEFORMAT

Ні

Виклик MMinor

j = 0

j <(siz-1)?

Так

y = yel?

Ні

Так

ky = 1

Вихід

i = 0

i <(siz-1)?

Так

x = xel?

Ні

Так

kx = 1

Ні

j + +

* (Res + j * (siz - 1) + i) = * (m + (j + kj) * siz + (i + ki));

i + +

Блок-схема функції MatrixMultiply з програми 30.c

SHAPE \ * MERGEFORMAT

Ні

Вихід

j = 0

j <rows1?

Так

i <cols2?

Так