// Для временного массива сбора данных 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; } Для работы программу необходимо скомпилировать в модели не ниже 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 знаков после запятой не отличаются от полученных. Задача решена верно. Разработать программу аппроксимации функции методом наименьших квадратов для модели по таблице результатов эксперимента: 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 и выводится на экран. SHAPE \* MERGEFORMAT Удалось выделить память для всех матриц F и FT? | Задание начальных значений в матрице F | Транспонирование матрицы F (TF=FT) | REV (aji) = Determinant(AC2) /Dt * ((i + j) % 2 == 0) ? 1:-1 | Вывод получившихся коэффициентов | SHAPE \* MERGEFORMAT *(res + j * (siz - 1) + i) = *(m + (j+kj) * siz + (i+ki)); | SHAPE \* MERGEFORMAT |