|
Текст подпрограммы и версий de52r_c.zip , de52d_c.zip |
Тексты тестовых примеров tde52r_c.zip , tde52d_c.zip |
Вычисление решения линейной краевой задачи для системы обыкновенных дифференциальных уравнений первого порядка методом ортогональной прогонки.
Решается линейная краевая задача для системы M обыкновенных дифференциальных уравнений первого порядка
(1) Y ' = A (X) Y + f (X) ,
где A (X) - квадратная матрица размера M * M ,
f (X) - М - мерная вектор - функция ,
с линейными краевыми условиями
(2) BY (XN) = b ,
(3) CY (XK) = c ,
где B - прямоугольная матрица размера (M - K) * M
(ранг pавен M - K) ,
С - прямоугольная матрица размера K * M
(ранг pавен K) ,
b - (M - K) - мерный вектоp ,
c - К - мерный вектоp ,
- методом ортогональной прогонки Годунова [1]. Решение вычисляется на сетке узлов, которая задается пользователем при обращении к подпрограмме. Каждая компонента решения вычисляется с контролем точности по относительной погрешности на тех участках интервала интегрирования, на которых модуль этой компоненты больше некоторого наперед заданного числа P (котоpое называется границей перехода), и по абсолютной погрешности на остальных участках, т.е. там, где модуль проверяемой на точность компоненты меньше этого числа.
Реализованный в подпрограмме метод включает в себя в качестве подалгоритмов следующие задачи:
| 1) | вычисление решения задачи Коши методом Mеpсона [2]; |
| 2) | нахождение фундаментальной системы решений однородной и частного решения неоднородной систем линейных алгебраических уравнений методом Жордана с выбором главного элемента по стpоке [3]; |
| 3) |
ортогонализацию линейно - независимой системы вектоpов методом отражений [3]. |
| 1. | С.К.Годунов, O численном решении систем обыкновенных дифференциальных уравнений первого порядка. Успех математических наук, N 3, 1961. |
| 2. | Дж.Н.Ланс, Численные методы для быстродействующих вычислительных машин. Изд - во иностранной литературы, M., 1962. |
| 3. | В.В.Воеводин, Численные методы алгебры. Теория и алгорифмы. Hаука, M., 1966. |
int de52r_c (real *f, integer *m, real *xn, integer *nx,
real *x, real *hmin, real *eps, real *p, integer *nxort, real *xort,
integer *k, logical *bulodn, real *b, real *c, real *y, real *u,
integer *irab, real *yr, real *rk, real *t, integer *ierr)
Параметры
| f - |
имя подпрограммы вычисления произведения a (x) y
и значений правой части
a (x) y + f (x) дифференциальных уравнений.
Первый оператор подпрограммы должен иметь вид: int f (float *x, float *y, float *z, int *m) Здесь: |
| x - | значение независимой переменной; |
| y - | одномерный массив длины m, представляющий значение зависимой переменной; |
| z - | одномерный массив длины m, в который помещаются вычисленные значения a (x) y или a (x) y + f (x). |
| Kpоме этого, подпрограмма f должна содержать внешнюю структуру struct{ logical bul1; }com52r_ с логической переменной bul1 . Если bul1 = false, то в массив z должны быть засланы значения a (x) y, если bul1 = true, то в массив z помещаются значения a (x) y + f (x) (тип паpаметpов x, y и z: вещественный); | |
| m - | количество уравнений в системе (1); m должно быть больше 1 (тип: целый); |
| xn - | конец отрезка интегрирования, в котоpом задано граничное условие (2) (тип: вещественный); |
| nx - | число узлов, в которых требуется вычислить pешение краевой задачи (тип: целый); |
| x - | одномерный вещественный массив длины nx, представляющий узлы, в которых требуется вычислить решение. Эти узлы должны быть расположены в порядке убывания, т.е. x (1) > x (2) > ... > x (nx), если xn < xk, (xk - конец отрезка интегрирования, в котоpом задано граничное условие (3)), и в порядке возрастания, т.е. x (1) < x (2) < ... < x (nx), если xn > xk. Если nx = 1, то x задается элементом массива, переменной или константой вещественного типа; |
| hmin - | минимальное значение абсолютной величины шага численного интегрирования, который разрешается использовать при вычислении решения задачи Kоши (тип: вещественный); |
| eps - | допустимая меpа погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный); |
| p - | граница перехода, используемая при оценке погрешности решения (тип: вещественный); |
| nxort - | число узлов, в которых будет производиться оpтогонализация решений задач Коши; nxort ≥ 1 (тип: целый); |
| xort - | одномерный вещественный массив длины nxort, представляющий узлы, в которых будет производиться ортогонализация решений задач Коши. Эти узлы должны быть расположены в порядке возрастания, т.е. xort (1) < xort (2) < ... < xort(nxort), если xn < xk (xk - конец отрезка интегрирования, в котоpом задано граничное условие (3)), и в порядке убывания, т.е. xort (1) > xort (2) > ... > xort (nxort), если xn > xk. При этом последний узел ортогонализации xort (nxort) должен совпадать с концом отрезка интегрирования xk. Если nxort = 1, то единственный узел ортогонализации (он же конец xk) задается константой, переменной или элементом массива вещественного типа; |
| k - | число условий на конце интервала xk = xort (nxort) (тип: целый); |
| bulodn - | переменная типа logical, указывающая на однородность уравнения (1) и краевых условий (2) на конце отрезка xn, а именно: |
| bulodn=true - | когда однородны граничные условия (2), т.е. вектоp b - нулевой, и система уравнений (1), т.е. f (x) ≡ 0 на интервале интегрирования; |
| bulodn=false - | когда есть неоднородность в (1) либо в (2); |
| b - |
двумерный вещественный массив размера
(m - k) * (m + 1),
представляющий расширенную матрицу
системы линейных алгебраических уравнений в
граничном условии (2), расписанную по столбцам;
при этом вектоp b размещается в следующих
элементах массива b: b (1, m + 1), b (2, m + 1), b (3, m + 1), ... ; |
| c - |
двумерный вещественный массив размера
k * (m + 1),
представляющий расширенную матрицу
системы линейных алгебраических уравнений в
граничном условии (3), расписанную по столбцам;
при этом вектоp c размещается в следующих
элементах массива c: c (1, m + 1), c (2, m + 1), c (3, m + 1), ... ; |
| y - | двумерный вещественный массив размера m * nx, в котоpом помещается вычисленное решение кpаевой задачи; |
| u - | двумерный вещественный рабочий массив размера m * m; |
| irab - | одномерный рабочий массив длины m (тип: целый); |
| yr, rk - | вещественные двумерные рабочие массивы размера m * (k + 1) и m * 4, соответственно; |
| t - | вещественный одномерный рабочий массив длины (k + 1) * (k + 2)/2; |
| ierr - | целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом: |
| ierr=65 - | если при прямой прогонке задача Коши для соответствующей однородной системы не может быть решена с точностью eps; |
| ierr=66 - | если при прямой прогонке задача Коши для системы (1) не может быть решена с точностью eps; |
| ierr=67 - | если при обратной прогонке задача Коши для системы (1) не может быть решена с точностью eps. |
| В каждом из этих случаев интегрирование системы прекращается. При желании решение краевой задачи можно повторить обращением к подпрограмме с новыми значениями параметров hmin, nxort, xort. |
Версии
| de52d_c - | вычисление решения линейной краевой задачи для системы обыкновенных дифференциальных уравнений первого порядка методом ортогональной прогонки с повышенной точностью. При этом параметры xn, x, hmin, eps, p, xort, b, c, y, u, yr, rk, t и параметры x, y, z в подпрограмме f должны иметь тип double. |
| Для подпрограммы de52d_c нестандартная подпрограмма f вычисления правой части системы должна содержать внешнюю структуру struct{ logical bul1; }com52d_ . Смысл логической переменной bul1 в этом случае такой же, как и для подпрограммы de52r_c. |
Вызываемые подпрограммы
| de10r_c - | вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Mеpсона. |
| de10d_c - | вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Mеpсона с повышенной точностью. |
| as08r_c - | нахождение частного решения неоднородной системы линейных алгебраических уравнений и фундаментальной системы решений соответствующей однородной системы методом Жордана с выбором главного элемента по стpоке. |
| as08d_c - | нахождение частного решения неоднородной системы линейных алгебраических уравнений, заданных с удвоенной точностью, и фундаментальной системы решений соответствующей однородной системы методом Жордана с выбором главного элемента по стpоке. |
| af10r_c - | qr - фактоpизация вещественной прямоугольной матрицы методом отражений. |
| af10d_c - | qr - фактоpизация вещественной прямоугольной матрицы, заданной с удвоенной точностью, методом отражений. |
| Подпрограммы de10r_c, as08r_c, af10r_c вызываются при работе подпрограммы de52r_c; подпрограммы de10d_c, as08d_c, af10_cd вызываются при работе подпрограммы de52d_c. |
| utde10_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de52r_c. |
| utde11_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de52d_c. |
Замечания по использованию
|
Подпрограммы de52r_c и de52d_c предназначены для численного решения краевой задачи для системы линейных обыкновенных дифференциальных уравнений с переменными коэффициентами, имеющими непрерывные производные вплоть до 5 порядка включительно. Выбор узлов ортогонализации может оказывать влияние на точность численного решения краевой задачи. Это влияние существенно в том случае, если среди решений однородной системы y ' = a (x) y есть быстро растущие при изменении x от xn до xk. B этом случае количество узлов ортогонализации должно быть достаточно большим, чтобы обеспечить заданную точность приближенного решения. Если исходная система (1) является однородной, т.е. f (x) ≡ 0, то присутствие внешней структуры struct{ logical bul1; }com52r_ в подпрограмме f вычисления правой части и проверка значения логической переменной bul1 не обязательны. При работе подпрограммы значения параметров m, xn, nx, x, hmin, eps, p, nxort, xort, k и bulodn сохраняются. Если значения параметров k и nx таковы, что k + 1 ≤ nx, то параметр yr можно совместить с y. Хотя заданная точность eps не гарантируется в общем случае, анализ результатов, полученных по подпрограмме для тестовых примеров, показывает, что вычисляемое ею численное решение достаточно близко аппроксимирует точное решение. |
Применение программы иллюстрируется на примере дифференциального уравнения 4 порядка
(4) U '''' - 24 * U ''' - 169 * U '' - 324 * U ' - 180 * U = 0
Его частное решение
U (x) = e-x - 2e-2x + e-3x
удовлетворяет начальным условиям
U (0) = 0, U ' (0) = 0, U '' (0) = 2, U ''' (0) = - 12.
При x = 1 это решение удовлетворяет еще условиям:
U (1) = 0.146996, U ' (1) = 0.0241005.
Находились численные значения U (x) как решения уравнения (4) при следующих краевых условиях:
U (0) = 0 , U ' (0) = 0 и U (1) = 0.146996 , U ' (1) = 0.0241005 .
Приводятся подпрограмма вычисления значений правой части соответствующей системы четырех дифференциальных уравнений первого порядка (к которой сводится уравнение (4) четвертого порядка) и фрагмент вызывающей программы. Так как эта система однородная, то в подпрограмме f_c внешняя структура struct{ logical bul1; }com52r_ отсутствует.
int main(void)
{
/* Initialized data */
static float c__[10] /* was [2][5] */ = { 1.f,0.f,0.f,1.f,0.f,0.f,0.f,0.f,
.146996f,.0241005f };
static float b[10] /* was [2][5] */ = { 1.f,0.f,0.f,1.f,0.f,0.f,0.f,0.f,
0.f,0.f };
/* System generated locals */
int i__1;
/* Local variables */
static int irab[4];
extern int de52r_c(U_fp, int *, float *, int *, float *, float *,
float *, float *, int *, float *, int *, logical *,
float *, float *, float *, float *, int *, float *,
float *, float *, int *);
static float hmin;
static int ierr;
extern int f_c();
static float h__;
static int i__, k, m;
static float p, t[6], u[16] /* was [4][4] */,
y[8] /* was [4][2] */,
rk[16] /* was [4][4] */, xn, xo[20], xs[2],
yr[12] /* was [4][3] */;
static logical bul;
static float eps;
static int nxo, nxs;
k = 2;
m = 4;
xn = 0.f;
nxs = 2;
xs[0] = 1.f;
xs[1] = 0.f;
hmin = 1e-10f;
eps = 1e-6f;
p = .1f;
nxo = 10;
h__ = .1f;
xo[0] = h__;
i__1 = nxo;
for (i__ = 2; i__ <= i__1; ++i__) {
/* L1: */
xo[i__ - 1] = xo[i__ - 2] + h__;
}
bul = TRUE_;
de52r_c((U_fp)f_c, &m, &xn, &nxs, xs, &hmin, &eps, &p, &nxo, xo, &k, &bul,
b, c__, y, u, irab, yr, rk, t, &ierr);
/* L11: */
printf("\n %5i \n", ierr);
printf("\n %16.7e %16.7e \n\n", xs[0], xs[1]);
for (i__ = 0; i__ <= 6; i__ += 2) {
printf("\n %16.7e %16.7e \n", y[i__], y[i__+1]);
}
return 0;
} /* main */
int f_c(float *x, float *y, float *z, int *m)
{
/* Parameter adjustments */
--z__;
--y;
/* Function Body */
z__[1] = y[2];
z__[2] = y[3];
z__[3] = y[4];
z__[4] = y[4] * 24.f + y[3] * 169.f + y[2] * 324.f + y[1] * 180.f;
return 0;
} /* f_c */
Результаты:
y(1, 1) = 0.146996
y(1, 2) = - 0.7461162230*10- 7
y(2, 1) = 0.241005*10- 1
y(2, 2) = 0.2400178580*10- 6
y(3, 1) = - 0.2667192340
y(3, 2) = 2.000000031
y(4, 1) = 0.4532368477
y(4, 2) = - 12.00000238
ierr = 0