2020-01-15
136
Для работы программу необходимо скомпилировать в модели не ниже 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.
Разработать программу аппроксимации функции методом наименьших квадратов для модели 
по таблице результатов эксперимента:
| X1 | X2 | Y |
| 1 | 1 | 0 |
| -1 | -1 | -2 |
| 2 | 2 | -2 |
| 3 | -2 | 29 |
| -2 | 4 | 54 |
Решение задачи по методу наименьших квадратов
Рассчитываемая модель линейна относительно своих коэффициентов ai. Задана матрицы 
и 
, а также функция для получения матрицы F. F – Специальная матрица, которая вычисляется по алгоритму, приведённому ниже. Функция представляет собой мою собственную разработку, но вполне возможно её вводить вручную. Алгоритм составления матрицы F (учитывая разложение 
):
, где 
- функции из модели y, а 
.- n-й элемент матрицы 
.
Исходя из этих формул строится функция f (смотри листинг программы 30.c).
Далее, по формуле 
находится матрица с коэффициентами ai и выводится на экран.
Блок-схема функции main из программы 30.c
|
Блок-схема функции MMinor из программы 30.c
Блок-схема функции MatrixMultiply из программы 30.c
Блок-схема функции Determinant из программы 30.c
Листинг программы 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;
}
