елемента результату
* (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
{
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
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
{
for (i = 0; i
{
// Заповнюємо 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
{
for (i = 0; i
{
