додати матеріал


Чисельні методи при вирішенні завдань

[ виправити ] текст може містити помилки, будь ласка перевіряйте перш ніж використовувати.


Курсова робота з інформатики
Тема: «Чисельні методи при вирішенні завдань»
Автор: студент групи ПС-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 (3) +2 y''+3 y '+ y = 5 + x 2 в інтервалі 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 (3) +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, Δ 2 q 1, Δ 2 q 0, Δ 3 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 1
X 2
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?
Так
i = 0
K + +
K <cols1?
Так
* (Result + cols2 * j +)) + =
* (M1 + cols1 * j + k) * (* (m2 + cols2 * k + i))

Ні
j + +
i + +
* (Result + (cols2 * j + i)) = 0
,
k = 0
Визод MatrixMultiply

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

SHAPE \ * MERGEFORMAT
Повернути d зухвалому
dim = 1?
Ні
Так
i <dim?
i = 0
d =* m
Так
d + = k * (* (m + i)) *
Determinant (mm, dim - 1); k = 0 - k;

Ні
Виділити пам'ять для мінору матриці m
Викреслити з матриці m елемент a (I, 0)
Звільнити пам'ять з-під мінору mm
Визод Determinant, d = 0
I + +

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

/ / Завдання 30. Апроксимація функції методом найменших квадратів
/ / (C) 2004 REPNZ
/ / Файли, що включаються
# Include <stdio.h>
# Include <conio.h>
# Include <dos.h>
# Include <stdlib.h>
/ / -------------- Опис початкових значень ------------------
/ / Дано (Розміри матриць - (1 х висота):
/ / Xm - це матриця-стовпець незалежних змінних
/ / Xm = (x1, x2, ... xN) T висотою xr
/ / Вектор спостережень. ym - його матриця:
/ / Ym = (y1, y2, ..., yM) T висотою yr
/ / А також опису функцій при коефіцієнти a1, a2, ..., aK
/ / 1. Матриці з елементами типу double
/ / - Кількість елементів у столбцевих марітцах xm і ym
# Define xr 2
# Define yr 5
/ / - Дані значення х
static double xm1 [xr] = {1, 1};
static double xm2 [xr] = {-1, -1};
static double xm3 [xr] = {2, 2};
static double xm4 [xr] = {3, -2};
static double xm5 [xr] = {-2, 4};
/ / - Масив покажчиків на ці значення
static double * xmp [yr] = {xm1, xm2, xm3, xm4, xm5};
/ / - Матриця зі значеннями функції
static double ym [yr] = {0, -2, -2, 29, 54};
/ / 2. Функції з моделі
/ / - Скільки їх
# Define n 3
/ / І власне самі функції, записуються як тіло Сі-функції
double f (double xm [xr], int path)
/ / - Які саме (n штук шляхів, вибирається параметром path)
{
switch (path)
{
/ / Функція 1
case 1:
return xm [0]; / / x1
/ / Функція 2
case 2:
return xm [1] * xm [1]; / / x2 ^ 2
/ / Функція 3
case 3:
return xm [0] * xm [1]; / / x1 * x2
}
printf ("\ nНеправільная функція \ n");
abort ();
}
/ / Ну і модель відповідно вийшла: y = a1 * x1 + a2 * x2 ^ 2 + a3 * x1 * x2
char txtmodel [] = "y = a1x1 + a2x2 ^ 2 + a3x1x2";
/ / Коротше, n = K, xr = N, yr = M (!) ;-)
////////////////////////////////////////////////// /////////////////////////////
/ / =-=-=-=-=-=-=-=-=-=-=-=-=-= Функції і підпрограми =-=-=-=-=-=-=-=-=- =-=-=-=
////////////////////////////////////////////////// /////////////////////////////
/ / Друк матриці m. Розміри (x * y)
void mprint (double * m, int x, int y)
{
int i, j; / / Індекси для проходу
for (j = 0; j <y; j + +) / / За рядками
{
for (i = 0; i <x; i + +) / / За елементами рядка
{/ / Елемент
printf ("% 8.4lg", * (m + (j * x + i)));
}
printf ("\ n"); / / CR / LF
}
}
////////////////////////////////////////////////// /////////////////////////////
/ / Перемноження матриць m1 (розмір - rows1 * cols1) і m2 (розмір - cols1 * cols2)
/ / Результат поміщається в result
void MatrixMultiply (double * m1, int rows1, int cols1, double * m2, int cols2, double * result)
{
int i, j, k;
/ / Вийде матриця висотою rows1 і довжиною cols2
for (j = 0; j <rows1; j + +) / / Прохід по висоті
{
for (i = 0; i <cols2; i + +) / / Прохід по довжині
{
/ / Очищення елемента
* (Result + (cols2 * j + i)) = 0;
for (k = 0; k <cols1; k + +) / / Прохід по елементах
/ / Рядки першої матриці
/ / Обчислення чергового елемента результату
* (Result + (cols2 * j + i)) + =
* (M1 + (cols1 * j + k)) * (* (m2 + (cols2 * k + i)));
}
}
}
////////////////////////////////////////////////// /////////////////////////////
/ / Обчислює мінор матриці m, отриманий викреслюванням елемента (xel; yel)
/ / І кладе його в res
void MMinor (double * m, double * res, int siz, int xel, int yel)
{
int i, j, ki = 0, kj = 0; / / Початковий стан
for (j = 0; j <(siz - 1); j + +) / / Проходимо по рядках матриці res
{
if (j == yel) kj = 1; / / Пропустити поточний рядок
for (i = 0; i <(siz - 1); i ++)// Проходимо по стовпцях матриці res
{
if (i == xel) ki = 1; / / Пропустити поточний стовпець
* (Res + j * (siz - 1) + i) = * (m + (j + kj) * siz + (i + ki));
}
ki = 0; / / Для наступної строчки (yel рядок вже пропустили)
}
}
////////////////////////////////////////////////// /////////////////////////////
/ / Обчислення визначника матриці m розміром (dim * dim)
/ / (Рекурсивна функція)
double Determinant (double * m, int dim)
{
/ / Всі змінні - ОБОВ'ЯЗКОВО ЛОКАЛЬНІ!
double d = 0, k = 1; / / Визначник і прапорець
int ki, kj, di, dj, i; / / Коефіцієнти, індекси, зміщення
double * mm; / / Нова матриця з викресленою рядком і стовпцем
if (dim <1) {printf ("\ nНеправільние аргументи"); abort ();}
if (dim == 1) return * m; / / Якщо матриця 1х1
/ / Виділяємо пам'ять для мінору
if ((mm = malloc ((dim - 1) * (dim - 1) * sizeof (double))) == 0)
{Printf ("\ nОшібка розподілу пам'яті \ n"); abort ();}
/ / Якщо матриця 2х2
if (dim == 2) d = ((* m) * (* (m + 3)) - (* (m + 2) * (* (m + 1 ))));
else / / Розмір більше ніж 2
/ / Розкладаємо матрицю за нульовою рядку
for (i = 0; i <dim; i + +)
{
MMinor (m, mm, dim, i, 0); / / Викреслимо стовпець і
/ / Рядок у матріцк
d + = k * (* (m + i)) * Determinant (mm, (dim - 1));
k = 0 - k;
}
free (mm); / / Звільнити пам'ять під мінор
return d; / / Повернути значення визначника
}
////////////////////////////////////////////////// /////////////////////////////
/ / Основна частина програмии
int main (void)
{
/ / Апроксимація функції для моделі y
double * F; / / Спеціальна матриця F n * y
double * TF; / / Транспонована F y * n
double * REV; / / Зворотній матриця n * n
double * TMP; / / Тимчасова матриця n * n
double * AC2; / / Алгебраїчні доповнення (n-1) * (n-1)
double dt; / / Значення визначника матриці
double flag; / / Прапорець для оберненої матриці
int i, j; / / Індекси
/ / Уявімо програму користувачеві:)
printf ("\ nПрограмма апроксимації функції методом найменших квадратів для"
"Моделі \ n% s"
"\ Nпо заданої таблиці експерименту."
"\ N \ n Розробник: студент групи ПС-146"
"\ N Нечаєв Леонід Володимирович"
"\ N 25.02.2004"
, Txtmodel);
printf ("\ nІзвестни результати спостережень:"
"\ N x1 x2 y");
for (i = 0; i <yr; i + +)
printf ("\ n% 10.4lg% 8.4lg% 8.4lg", * (xmp [i]), * (xmp [i] + 1), ym [i]);
printf ("\ nНачінаем апроксимацію ... \ n");
/ / Требуется порахувати am. Так:
/ / Am - це матриця-стовпець шуканих коефіцієнтів. Представляє з себе
/ / Am = (a1, a2, ..., aK) T висотою n, а вважається так:
/ / Am = Inverse [Transpose [F]. F]. Transpose [F]. Ym, де
/ / F - мартицей, складена спеціальним чином (дивись нижче):
/ / Виділяємо пам'яті відразу на всі матриці - F, TF, REV, TMP, AC2
# Define memneed (((n * yr) + (yr * n) + (n * n) + (n * n) + ((n-1) * (n-1))) * eof (double))
if ((F = malloc (memneed)) == 0)
{
printf ("\ nОшібка розподілу пам'яті. Замініть комп'ютер");
abort (); / / Якщо не вдалося виділити для неї пам'ять
}
TF = F + (n * yr);
REV = TF + (yr * n);
TMP = REV + (n * n);
AC2 = TMP + (n * n);
/ / Заповнення значеннями матриці F
for (j = 0; j <yr; j + +) / / Цикл по рядках F
{
for (i = 0; i <n; i + +) / / І по стовпцях F
{
/ / Заповнюємо j-й рядок значеннями функцій fi
* (F + (j * n + i)) = f (xmp [j], (i + 1));
}
}
/ / Матриця F готова. Треба обчислити за формулою:
/ / Am = Inverse [Transpose [F]. F]. Transpose [F]. Ym значення
/ / Коефіцієнтів a1, a2, a3, ...
/ / Транспонуємо F
for (j = 0; j <n; j + +) / / Цикл по рядках TF
{
for (i = 0; i <yr; i + +) / / І по її стовпцями
{
* (TF + (j * yr + i)) = * (F + (i * n + j));
}
}
/ / Вважаємо TMP = TF * F
MatrixMultiply (TF, n, yr, F, n, TMP);
/ / Далі вважаємо оперделітель від TMP
if ((dt = Determinant (TMP, n)) == 0)
{
printf ("\ nТак, як визначник матриці TF * F дорівнює 0, \ n"
"Неможливо порахувати зворотний до них матрицю \ n");
free (F); abort ();
}
/ / Складаємо зворотну матрицю.
for (j = 0; j <n; j + +)
{
for (i = 0; i <n; i + +)
{
/ / Беремо Мінор елемента ij
MMinor (TMP, AC2, n, i, j);
/ / Знак елемента
flag = ((i + j)% 2 == 0)? 1. : -1.;
/ / Відразу транспонування
* (REV + (i * n) + j) = flag * Determinant (AC2, (n - 1)) / dt;
}
}
/ / Множимо обернену матрицю на транспоновану до F
/ / Тобто Inverse (TF * F) * TF
/ / Така матриця буде розміру yr * n, тому цілком вистачить пам'яті для F
MatrixMultiply (REV, n, n, TF, yr, F);
/ / І, нарешті, все це множимо на матрицю Y і отримуємо шукані
/ / Коефіцієнти a1, a2, ... aN
/ / Для такої матриці (розміром 1 * n) цілком вистачить пам'яті під REV
MatrixMultiply (F, n, yr, ym, 1, REV);
/ / Все, друкуємо відповідь
printf ("\ nВичісленія успішні, отриманий наступні коефіцієнти:");
for (i = 0; i <n; i + +)
printf ("\ na% d =% lg", i, * (REV + i));
/ / Звільнити пам'ять
free (F);
printf ("\ nНатисніть any key");
getch ();
printf ("\ nDone. \ n");
return 0;
}

Резуль

Додати в блог або на сайт

Цей текст може містити помилки.

Програмування, комп'ютери, інформатика і кібернетика | Курсова | 65.8кб. | скачати

Схожі роботи:
Розвиток логічного мислення учнів при вирішенні завдань на побудову
Актуалізація різного типу знань при вирішенні психологічних завдань
Метод моделювання розвитку психічної діяльності при вирішенні навчальних та ігрових завдань
Методи оптимізації при вирішенні рівнянь
Місцевий бюджет та його роль у вирішенні соціально-економічних завдань Адміністрації Красноармійського
Чисельні методи 6
Чисельні методи 5
Чисельні методи 3
Чисельні методи 4
© Усі права захищені
написати до нас
Рейтинг@Mail.ru