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

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

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

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).

Порядок решения задачи.

Для решения задачи вводятся стандартная и рабочая информативы и директива.

В рабочей информативе после метки Ц программа вычисления правых частей системы. Здесь Z[1]=…; Z[2]=…; …;Z[n]=…; - правые части исходной системы обыкновенных дифференциальных уравнений как функции от X1 и Y[1], Y[2], …, Y[n], X1 — соответствует аргументу, Y[I] - соответствует функциям. I=1, 2, …, N. Операторная часть рабочей информативы заканчивается оператором перехода «НА» Ф.

В описательной части рабочей информативы задаются X0, XK — соответственно начало и конец отрезка интегрирования, Q -шаг интегрирования методом Хемминга, J — число, определяющее, во сколько раз следует уменьшить шаг интегрирования методом Рунге-Кутта на участке «разгона» для получения решения того же порядка точности, что и в методе Хемминга,

N=n — порядок системы;

Y[n] - вектор начальных условий,

W[n] - вектор коэффициентов для вычисления невязки

W[I]=1, и описаны

A[n], B[n], C[n] - массивы значений функций в точках i-3,

i-2, i-1 соответственно,

Я[n], Б[n], Г[n], D[n] - массивы значений производных в точках i-3, i-2, i-1, i соответственно, Z[n] - массив правых частей,

П[n], P[n] - рабочие массивы.

В директиве задаются: R — разрядность вычислений по методу Хемминга («разгон» происходит с увеличенной разрядностью), Ю — число, определяющее период печати (количество шагов). Директива должна оканчиваться оператором «НА» HMG.

Описание работы программы

Данная расчетно-графическая работа (далее РГР) составлена на языке PC MathLab (PC-MATLAB © Copyright The MathWorks, Inc. 1984−1989 Version 3.5f 17-July-89 Serial Number 22 961) и выполнена в виде двух модулей (третий — контрольный пример), распечатка которых приведена в приложении.

1. Hemming. m

«Стандартный» головной модуль.

Входные данные: отсутствуют.

Выходные данные: отсутствуют.

Язык реализации: PC MathLab.

Операционная система: MS-DOS 3.30 or higher

Пояснения к тексту модуля:

Структура данного модуля элементарна. Вначале очищается экран, задаются исходные данные для второго модуля, как X0, XK — начальное и конечное значение, Q — шаг, J — число, определяющее во сколько раз нужно уменьшать шаг интегрирования методом Рунге-Кутта (далее Р-К) на участке «разгона» для получения того же порядка точности, что и в методе Хемминга, N — порядок системы, Y — вектор начальных значений, W — вектор коэффициентов для вычисления невязки и т. д. Затем вызывается модуль решения системы в формате:

[x, y, dg]=hem ('ours', x0, xk, q, j, n, y, w, ur), где

x, y — точки решения

dg — ошибка остальные параметры описаны выше.

Необходимо отметить, что несмотря на отсутствие входных и выходных данных, внутри данного модуля задаются начальные значения и выводятся результаты вычислений в числовом виде и графиков, а также оценка по быстродействию (TIME) и количеству выполненных операций (FLOPS), однако эти данные нельзя охарактеризовать как входные и выходные.

2. Hem. m

Модуль, которые непосредственно и решает систему ОДУ ме

тодом Хемминга.

Входные данные:

FunFcn — имя функции, вычисляющей левые части

X0,XK — начальное и конечное значение для счета

Q — шаг интегрирования

J — число, определяющее во сколько раз нужно уменьшать шаг интегрирования методом Рунге-Кутта (далее Р-К) на участке «разгона» для получения того же порядка точности, что и в методе Хемминга

N — порядок системы

Y — вектор начальных значений

W — вектор коэффициентов для вычисления невязки

UR — число, определяющее период печати

Выходные данные:

x — матрица точек, для которых вычислено решение

y — матрица решений

dg — ошибка интегрирования

Язык реализации: PC MathLab.

Операционная система: MS-DOS 3.30 or higher

Пояснения к тексту модуля:

Данный модуль содержит в своем теле всего одну функцию, входные и выходные данные которой являются входными и выходными данными текущего модуля. Они описаны выше. Мы же займемся описанием данной функции: