// вычисления 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
{ // Вычислим наконец-то 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
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])
{ // Все приращения
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
// по отношению к текущему
{ 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;
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
0
-1
-2
2
3
29
4
54
Рассчитываемая модель линейна относительно своих коэффициентов ai. Задана матрицы и , а также функция для получения матрицы F. F - Специальная матрица, которая вычисляется по алгоритму, приведённому ниже. Функция представляет собой мою собственную разработку, но вполне возможно её вводить вручную. Алгоритм составления матрицы F (учитывая разложение ):
, где - функции из модели y, а .- n-й элемент матрицы .
Исходя из этих формул строится функция f (смотри листинг программы 30.c).
Далее, по формуле находится матрица с коэффициентами ai и выводится на экран.
Министерство образования Российской Федерации
Южно-Уральский государственный университет
Кафедра Систем управления
Челябинск, 2004
Страницы: 1, 2, 3
