--PAGE_BREAK-- // Для временного массива сбора данных
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
f (y0, q0, xi); // Заполняем q0, основываясь на значениях из y0
for (j = 0; j
xi += h; // Следующий шаг
//… а остальные 3 добываем с помощью метода Рунге-Кутта 4 порядка.
for (i = 0; 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
{
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
{ // Вычислим наконец-то 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
{
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
for (j = 0; 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
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
{ // Все приращения
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
// Очередное значение функции вычислено. Следующиий шаг
xi += h;
// Продолжить интегрирование?
if (xi
// Если первый раз здесь, то просчитать ещё раз с шагом h/2
if (flag == 0)
flag = 1; // Сравнивать уже будет с чем
else
{
// Не первый раз — оценить погрешность
// Сейчас в y3 — значение только что вычисленной функции ,
// а в y2 — занчение функции, вычисленной с шагом h * 2
// по отношению к текущему
for (j = 0; j
{ eps2 = fabs (((y3[j] — y2[j]) / y2[j]));
if (eps2 > eps) break; // Если погрешность слишком великА
}
if (j == n) // Если всё ОК
{ // Копируем результат
for (j = 0; 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
{
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 SHAPE \* MERGEFORMAT
Блок-схема функции MMinor из программы 30.c SHAPE \* MERGEFORMAT
Блок-схема функции MatrixMultiply из программы 30.c SHAPE \* MERGEFORMAT продолжение
--PAGE_BREAK--