Программа для численного решения системы обыкновенных дифференциальных уравнений

Теоретическая часть

Программа предназначена для численного решения системы обыкновенных дифференциальных уравнений вида:

Y'=F (X, Y), с начальными условиями Y (X0)=Yo на отрезке [X, X] методом Хемминга с постоянным шагом интегрирования. В каждой i+1 точке находим начальное приближение Р к решению Y по предсказывающей формуле:

Pi+1=Yi-3+4*h*(2*Y'i-Y'i-1+2*Y'i-2)/3, где

Yi-3 решение в i-3 точке,

Y’i, Y'i-1,Y'i-2 — значения производных в точках i, i-1,i-2 соответственно.

Для улучшения решения используется корректирующая формула

Ci+1=[9*Yi-Yi-2+3*h*(M'i+1+2*Y'i-Y'i-1)]/8, где

Mi+1=Pi+1−112*(P-Ci)/121; M’i+1=F (Xi+1,Mi+1).

Решение системы в i+1 точке находится по формуле

G=Wj*|Pi+1,j-Ci+1|, где

Wj=1

j- номер компоненты вектора.

На участке «разгона» значения Yi-k и Y’i-k (k=0, 1, 2) вычисляются методом Рунге-Кутта по формуле Yi=Ui (2)-(Ui (i)-Ui (2))/15, где i- номер точки, в которой ищется решение, Ui- решение системы в i-ой точке, полученное с шагом h/l;

U (l)i-m+1/l=A (l)i-m+(K (l)1+2*K (l)2+2*K (l)3+K (l)4)i--m+1/l/6, где

j=1, 2, …, n,

K (l)1=h*F (Xi-m, A (l)i-m)/l;

K (l)2=h*F (Xi-m+h/2*l, A (l)i-m+K (l)½)/l;

K (l)3=h*F (Xi-m+h/2*l, A (l)i-m+K (l)2/2)/l;

K (l)4=h*F (Xi-m+h/l, A (l)i-m+K (l)3)/l.

A, U, K — векторы n-го порядка; l=1, 2; m=1 при l=1; m=1,½ при l=2;

A (l)i-1=Y (l)i-1; A (2)i-½=U (2)i-½.

Характеристика программы.

Программа состоит из стандартной информативы, реализующей описанный метод, рабочей информативы, задающей правые части уравнений системы и директивы.

Длина стандартной информативы 1600 символов. Объем исходных данных: 7 чисел, 2 массива, n функций. В результате работы программы на печать выводится на участке «разгона» X, значения функций и производных, далее X, G и Y[n] на всем отрезке интегрирования через Ю шагов и в конце отрезка.

Программа рекомендуется для решения систем обыкновенных дифференциальных уравнений на больших отрезках, так как считает быстрее одноточечных методов. Для контроля постоянно выводится погрешность вычислений G, которая позволяет следить за точностью решения.

«Разгон» (нахождение значений функций и производных в точках X0, X0+Q, X0+2*Q, X0+3*Q, где Q — шаг интегрирования) осуществляется методом Рунге-Кутта с увеличенной разрядностью.

В программе предусмотрена возможность при получении большой погрешности вычисления в точка «разгона» уменьшить шаг интегрирования в этих точках (см. способ задания J), а при быстром возрастании погрешности вычислений G уменьшить шаг интегрирования методом Хемминга или увеличить разрядность вычислений.

Программа позволяет производить интегрирование как с положительным, так и с отрицательным шагом (соответственно меняются X0, Xк и Q).