Помогите разобраться и переделать программу.
Решение ДУ методом Адамса
Есть код программы для:
А мне нужно переделать на;
Прошу помощи, сам не могу разобраться, не получается.
Код:
#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
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 = (double*)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 порядка - //
///////////////////////////////////////////////////////////////////////
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; // Следующий шаг }