Южно-Уральский государственный университет
Кафедра Систем управления
Челябинск, 2004
// Задача 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Вычисления успешны, получен следующие коэффициенты:");
printf ("\na%d = %lg", i, *(REV + i));
// Освободить память
free (F);
printf ("\nНажмите any key");
getch ();
printf ("\nDone.\n");
return 0;
После запуска программа сразу же начинает расчёт коэффициентов. На экран выводится следующее:
Программа аппроксимации функции методом наименьших квадратов для модели
y = a1x1 + a2x2^2 + a3x1x2
по заданной таблице эксперимента.
Разработчик: студент группы ПС-146
Нечаев Леонид Владимирович
25.02.2004
Известны результаты наблюдений:
x1 x2 y
1 1 0
-1 -1 -2
2 2 -2
3 -2 29
-2 4 54
Начинаем аппроксимацию...
Вычисления успешны, получены следующие коэффициенты:
a0 = 1
a1 = 2
a2 = -3
Нажмите any key
Done.
Результат верен, так как при подстановке переменных в модель получается верное равенство:
0 = 1 * 1 + 2 * 1 - 3 * 1 * 1
-2 = 1 * (-1) + 2 * (-1) - 3 * (-1) * (-1)
-2 = 1 * 2 + 2 * 2 - 3 * 2 * 2
29 = 1 * 3 + 2 * (-2) - 3 * 3 * (-2)
54 = 1 * (-2) + 2 * 4 - 3 * (-2) * 4
Задача решена верно.
Страницы: 1, 2, 3
