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

ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

Задаче Коши в разделе численных методов решения дифференциальных уравнений уделяется основное внимание.

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

Среди них особое место занимает семейство методов Рунге-Кутта, относительно которых сегодня мы с полным правом можем сказать, что это - классика.

Однако при решении конкретных задач с помощью этих методов бывают непредвиденные неприятности типа аварийного останова программы с выдачей сообщения о переполнении, или же процесс решения зацикливается.

При попытке найти полезный совет в тех же учебниках, там либо «молчание», либо вместо совета читаем: «Для эффективного решения задачи большое значение имеет опыт, интуиция и квалификация исследователя как при постановке задачи, так и в процессе выбора метода, разработки алгоритма и программы решения задачи на ЭВМ» ([1], с.102).

Согласитесь, что эта фраза мало может помочь студенту или начинающему специалисту.

Итак, предметом исследования настоящей работы являются численные методы решения таких задач Коши, когда классика не дает ожидаемого результата.

Формулировка задачи Коши в этом случае ничем не отличается от классической: требуется решить с заданной точностью на отрезке [х, х+b] уравнение:

y'=F(x,y) при y(х)=y (1)

Пожалуй, наиболее популярным методом ее численного решения считается обычный метод Рунге-Кутта, которому вместе с высокой степенью наглядности присущи такие положительные черты, как возможность применения переменного шага и достаточно малая погрешность получаемого решения (порядка h, где h- шаг, см., например, [2], с.102).

При решении конкретных задач этим классическим методом (как, впрочем, и другими) бывают непредвиденные неприятности, о которых мы говорили выше. Если мы сталкиваемся с зацикливанием, то это, как правило, связано с недопустимо большим дроблением шага h, впрочем есть и другие причины. Иногда (при переполнении) помогает совет стартовать с более мелким шагом, но не всегда. Анализ причин получаемых неприятностей при численном решении конкретных задач приводит к формулировке следующих четырех вопросов:

1. Насколько корректной является общепринятая трактовка определения точности численного решения дифференциального уравнения?

2. Если эта трактовка не всегда корректна, то что можно предложить взамен ей?

3. Всегда ли мы можем надеяться на успешный поиск численного решения задачи Коши на интервале [х, х+b]?

4. Как же на плоскости xOy с помощью численного метода строить интегральные кривые, которые при различных начальных условиях будут являться решением рассматриваемой задачи Коши?

До того как рассматривать эти вопросы и возможные ответы на них, с помощью простого примера покажем, что указанные выше проблемы действительно существуют. Рассмотрим следующую задачу Коши: для уравнения (2):

y'= - (2)

с начальным условием y(0)=1 попытаемся найти решение на интервале [0, 2] с заданной точностью (например, =0,0001).

Написать программу для реализации обычного метода Рунге-Кутта нетрудно. Однако если мы пойдем прямой дорожкой (как нам советует метод), то решения от компьютера мы не получим. Нетрудно видеть (уравнение (2) легко решается аналитически путем разделения переменных), что неявное задание функции y от х в виде равенства хІ+yІ=1 и определяет искомую интегральную кривую для задачи Коши (2). Это - окружность, и если она проходит через точку с координатами (0,1), что соответствует нашему начальному условию, то точек с координатами х>1 на ней просто нет. Потому искать решение мы можем только на отрезке [0,1], а никак не на отрезке [0,2]. Но и в таком случае нам компьютер не выдаст решение, т.к. при х=1 имеем y=0 и F(x,y) = -, а символ компьютер, как известно, “не переваривает”, останавливаясь и выдавая сообщение о переполнении.

Для задачи (2) мы можем легко определить ту длину интервала (величину b), где можно искать решение, но сколько нам может встретиться задач, когда это сделать будет не так легко?! Если же задавать длину интервала «от потолка», то это для многих задач может повлечь зацикливание программы, если ранее не встретится другой не менее приятный вариант переполнения.

Искомую интегральную кривую мы находим в виде системы точек {М} c координатами {X,Y}, при этом традиционно сравниваем значения Y и , полученные при шаге h и шаге h/2, считая, что достигнута требуемая точность, как только:

|Y-| < (3).

Когда мы говорим о точности полученного решения, то подразумеваем, что речь идет о близости двух кривых - точного и полученного решений задачи Коши. Однако условие (3) далеко не всегда может быть инструментом, адекватно соответствующим поставленной цели, т.е. не всегда может корректно отражать близость тех или иных кривых. Очевидно, что в случаях, когда интегральные кривые резко «подымаются» вверх или, наоборот, резко «падают» вниз, иными словами, при выполнении условия |y'| >>1, условие (3) как инструмент, с помощью которого делается заключение о близости двух указанных выше кривых, применять нельзя.

Таким образом, пришли к формулировке первого из указанных выше вопросов и, более того, к утверждению того, что общепринятая трактовка определения точности численного решения дифференциального уравнения не всегда корректна. Там, где решения получаются достаточно гладкими с ограниченной (по модулю) величиной y', все нормально, а иначе мы с большой вероятностью столкнемся с проблемами при отыскании решения нашей задачи. Коль так, то естественна формулировка второго из приведенных выше вопросов.

С тем чтобы предложить альтернативу условию (3) вместе с задачей (1) (при ее решении), будем рассматривать родственную ей задачу:

x'=(F(x,y))-1 при x(y)=x (4)

при этом пока ?y<= ?x, где ?y = | Y-Y| и ?x = | X-X|, решение ищем для задачи (1). Как только выполнится условие ?y>?x, переходим к дальнейшему поиску интегральной кривой при помощи задачи (4) с надлежащей корректировкой начального условия (ниже мы этот алгоритм выпишем более строго) и продолжаем этот поиск, пока это условие (?y>?x) выполняется. В противном случае снова переходим к поиску с помощью формулировки задачи (1), опять таки с соответствующей корректировкой начального условия.

Далее все повторяется аналогично. Подчеркнем: мы пользуемся условием (3) лишь при том, что углы наклона рассматриваемых кривых (по отношению к оси Ох) меньше р/4, иначе мы рассматриваем условие (5):

|Х-| < (5)

естественно до тех пор, пока угол наклона кривых по отношению к оси Oy будет меньше чем р/4. Если у вас возник вопрос, почему именно углы наклона не должны превышать значения р/4, то ответ прост: потому что при этом максимальная погрешность при определении близости интегральных кривых как с использованием критерия (3), так и с использованием критерия (5) будет одинакова.

В противном случае при одном из этих критериев нами будет допускаться максимальная погрешность, гораздо больше той, которая допускается нами с углом наклона, равным р/4.

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

Впрочем, благодаря предложенному здесь алгоритму решения рассматриваемой задачи Коши, мы можем надеяться, что тем самым практически избежали неприятности, связанной с переполнением, а что касается зацикливания, то заметим, что недопустимо большое дробление шага h, приводящее к зацикливанию, часто напрямую связано с ростом величины производной y'. Рассмотренный здесь подход также позволяет избежать и зацикливания (если оно возникает по указанной причине).

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

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

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

Итак, сформулируем основные положения предложенного метода решения задачи Коши:

1. Метод может опираться на любой известный пошаговый метод, но рекомендуется все же обычный метод Рунге-Кутта.

2. Вместе с задачей Коши (1) (решить с заданной точностью на отрезке [х, х+bx] уравнение y'=F(x,y) при y(x)=y) рассматривается родственная задача (4) (решить с заданной точностью на отрезке [y, y+by] уравнение x'=(F(x,y))-1 при x(y)=x). При этом задачу (1) избранным пошаговым методом решаем только при условии, что |y'|<1, иначе переходим к задаче (4), которую решаем до тех пор, пока |x'|<1, иначе переходим к решению задачи (1).

3. Стартовый момент (какую из задач (1) или (4) решать первой) определяем путем вычисления значения производной y' в начальной точке, т.е. y'=F(x,y) и сравнения этого значения с единицей. Если |y'|<1, то первой решается задача (1), иначе - первой решается задача (4).

4. В соответствии с формулировкой задач (1) и (4), данными в п.2 эти задачи решаются только начиная со стартового момента (точнее, только одна из этих задач решается в указанной формулировке), а дальше формулировки меняются: для задачи (1): решить с заданной точностью на отрезке [х, х+bx] уравнение y'=F(x,y) при y(x)=y; для задачи (4): решить с заданной точностью на отрезке [y, y+by] уравнение x'=(F(x,y))-1 при x(y)=x, где x,y - координаты последней точки искомой интегральной кривой, вычисленные по предыдущей задаче. Значения bx и by определяются очевидным образом из условий:

х+bx =х+bx, y+by = y+by (6)

В приведенной в начале данной статьи цитате говорилось, что эффективность решения задачи зависит не только от того, как разработан алгоритм ее решения, но и от того, как составлена программа решения задачи на ЭВМ. Увы, тонкости этого вопроса выходят за рамки данной статьи. Вместе с тем считаем необходимым здесь уделить этому важному вопросу хоть какое-то внимание.

С этой целью проиллюстрируем предложенный метод решения задачи Коши текстом программы на языке С для решения рассмотренной выше задачи (2). Решение выдается в виде координат точек искомой кривой. Для того чтобы это не было слишком быстро, после каждой печати предусмотрена пауза, которая заканчивается тут же после ввода любого символа с клавиатуры, кроме символа `q', ввод которого связан с концом процесса расчета. Целочисленная величина ind, принимающая возможные значения 1, 0, -1, помогает при необходимости менять знак для шага h (иначе окружность не получить!).

#include<stdio.h>

#include<math.h>

#include<conio.h> // для использования небуферизованного ввода

# define hh 0.1f // начальное значение шага h

float f(float x, float y); // расчет F(x,y)

float f_1(float x, float y); // расчет (F (x,y))О№

float runge_(float x, float y); // расчет по методу Рунге на одном шаге задачи (1)

float runge_1(float x, float y); // расчет по методу Рунге на одном шаге задачи (4)

int runge_y(float * xx, float * yy, float bx, float e); // расчет по методу Рунге задачи (1)

int runge_x(float * xx, float * yy, float by, float e); // расчет по методу Рунге задачи (4)

float h; // шаг - глобальная переменная

int main(void)

{ int ind;

float x,y,x0,y0,bx,by,e;

printf("введите x0,y0,bx,by,e\n");

scanf("%f%f%f%f%f",&x0,&y0,&bx,&by,&e); // вводим: 0 1 1 1 0.0001 «Enter»

printf(" x y\n");

x=x0;

y=y0;

printf("1:%f %f\n",x,y);

ind=1;

ddd:

if (ind) {h=ind*hh; ind=runge_y(&x, &y, bx, e);} else return 0;

if (ind) {h=ind*hh; ind=runge_x(&x, &y, by, e);} else return 0;

goto ddd; }

int runge_y(float * xx, float * yy, float bx, float e)

{ float x,y,y1,y2,y3,x1;

char ch;

x=*xx;

y=*yy;

do {

do { y1=runge_(x,y);

h=h/2.0f;

y2=runge_(x,y);

x1=x+h;

y3=runge_(x1,y2);

} while (fabs(y3-y1)>e);

h=h*2.0f;

x=x+h;

y2=y; y=y3;

printf("%f %f\n",x,y);

*xx=x;

*yy=y;

ch=getch(); // вводим любой символ (для конца расчета - символ q)

if (ch= ='q') return 0;

} while (fabs(y3-y2)<fabs(h) && fabs(x)<bx);

if (fabs(y3-y2)>=fabs(h)) { if (y3-y2<0) {return (-1);} else {return (1);} } else return 0; }

float f(float x, float y)

{ return (-x/y); }

float runge_(float x, float y)

{ float k0,k1,k2,k3;

k0=h*f(x,y);

k1=h*f(x+h/2,y+k0/2);

k2=h*f(x+h/2,y+k1/2);

k3=h*f(x+h,y+k2);

return (y+(k0+2*k1+2*k2+k3)/6); }

float f_1(float x, float y)

{ return (-y/x); }

float runge_1(float x, float y)

{ float k0,k1,k2,k3;

k0=h*f_1(x,y);

k1=h*f_1(x+0.5*k0,y+0.5*h);

k2=h*f_1(x+0.5*k1,y+0.5*h);

k3=h*f_1(x+k2,y+h);

return (x+(k0+2*k1+2*k2+k3)/6); }

int runge_x(float * xx, float * yy, float by, float e)

{ float x,y,y1,x1,x2,x3;

char ch;

x=*xx;

y=*yy;

do {

do { x1=runge_1(x,y);

h=h/2;

x2=runge_1(x,y);

y1=y+h;

x3=runge_1(x2,y1);

} while (fabs(x3-x1)>e);

h=h*2.0f;

y=y+h;

x2=x; x=x3;

printf("%f %f\n",x,y);

*xx=x;

*yy=y;

ch=getch(); // вводим любой символ (для конца расчета - символ q )

if (ch= ='q') return 0;

} while (fabs(x3-x2)<fabs(h) && fabs(y)<by);

if (fabs(x3-x2)>=fabs(h)) { if (x3-x2<0) {return (-1);} else {return (1);}} else return 0; }

Выводы:

1. В ряде случаев классические численные методы решения обыкновенных дифференциальных уравнений (задачи Коши) не дают ожидаемого результата.

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

3. Одним из таких рецептов может служить приведенный в данной статье алгоритм.

СПИСОК ЛИТЕРАТУРЫ

1. Маликов В.Т., Кветный Р.Н. Вычислительные методы и применение ЭВМ. -К.: Головное издательство издательского объединения «Выща школа», 1989.

2. Методи обчислень. Практикум на ЕОМ: Навч. посібник / В.Л. Бурковська, С.О. Войцехівський, І.П. Гаврилюк та інші. -К.: Вища школа, 1995. -303 с.