.2004
# include
# include
# include
void func (double * y, double * ys, double t)
{//функція обчислення правих частин рівнянь
ys [0] = y [1];// ys [1]-перша похідна; ys [2]-друга і т.д.
ys [1] = y [2];// t-незалежний аргумент
ys [2] = 5 + t * t - y [0] - 3. * Y [1] - 2. * Y [2];
}
void Adams (
void f (double * y, double * ys, double x),
// Функція вичіленія правих частин системи
double * y,// ​​Масив розміру n значень залежних змінних
int n,// ​​Масив розміру n значень похідних
double tn,// ​​Початок інтервалу інтегрування
double tk,// ​​ʳнець інтервалу інтегрування
int m,// ​​Початкове число розбиття відрізка інтегрування
double eps)// Відносна похибка інтегрування
{
double * k1, * k2, * k3, * k4;// Для методу Рунге-Кутта
double * q0, * q1, * q2, * q3;// Значення похідних Для методу Адамса
double * ya;// Тимчасовий масив
double * y0, * y1, * y2, * y3;// Значення функції для методу Адамса
double h;// Крок інтегрування
double xi;// Поточне значення незалежної змінної
double eps2;// Для оцінки похибки
double dq2, dq1, dq0, d2q1, d2q0, d3q0;// прирощення
int flag = 0;// 0, поки йде перший прорахунок
int i, j;// Індекси
В
if (m <4) m = 4;// Мінімум 4 відрізка
if (tn> = tk)
{printf (" nНеправільние аргументи n");
abort ();// Неправильні аргументи
}
// Виділяємо пам'ять для масивів з змінними
if ((k1 = malloc ((4 + 4 + 4 + 1) * n * sizeof (double))) == 0)
{printf (" nОшібка розподілу пам'яті n");
abort ();// Перервати, якщо не вдалося
}
// Распределяем пам'ять між масивами
// Для методу Рунге-Кутта 4 порядку
k2 = k1 + n; k3 = k2 + n; k4 = k3 + n;
// 4 пердидущіх значення функції
y0 = k4 + n; y1 = y0 + n; y2 = y1 + n; y3 = y2 + n;
// Для тимчасового масиву збору даних
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 порядку -// p>///////////////////////////////////////////////////////////////////////
for (j = 0; j
f (y0, q0, xi);// Заповнюємо q0, грунтуючись на значеннях з y0
for (j = 0; j
xi + = h;// Наступний крок
// ... а інші 3 видобуваємо за допомогою методу Рунге-Кутта 4 порядку.
for (i = 0; i <3; i + +)// i - ЯКЕ ЗНАЧЕННЯ ВЖЕ Є
{//А Обчислюється значення Y [i +1]!! p>// Спочатку потрібні коефіцієнти 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.; p>// І один з аргументів для функції p>}// обчислення 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;// Наступний крок}
///////////////////////////////////////////////////////////////////////
// - Метод Адамса -// p>//////////////////////////////////...
