// Для временного массива сбора данных
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 |
