p align="left"> // Очередное значение функции вычислено. Следующиий шаг 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.Разработать программу аппроксимации функции методом наименьших квадратов для модели по таблице результатов эксперимента: |
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Министерство образования Российской Федерации Южно-Уральский государственный университет Кафедра Систем управления Челябинск, 2004 Блок-схема функции MMinor из программы 30.cМинистерство образования Российской Федерации Южно-Уральский государственный университет Кафедра Систем управления Челябинск, 2004 Блок-схема функции MatrixMultiply из программы 30.cМинистерство образования Российской Федерации Южно-Уральский государственный университет Кафедра Систем управления Челябинск, 2004 Блок-схема функции Determinant из программы 30.cМинистерство образования Российской Федерации Южно-Уральский государственный университет Кафедра Систем управления Челябинск, 2004 Листинг программы 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; } Результат решения задачи 30 на ЭВМПосле запуска программа сразу же начинает расчёт коэффициентов. На экран выводится следующее: Программа аппроксимации функции методом наименьших квадратов для модели 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
|