О решении жестких уравнений с помощью явных схем

О РЕШЕНИИ ЖЕСТКИХ УРАВНЕНИЙ С ПОМОЩЬЮ ЯВНЫХ СХЕМ

А.А. Литвиненко, доц. Сумской государственный университет.

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

Стандартные численные методы решения обыкновенных дифференциальных уравнений (ОДУ) и их систем в ряде случаев бессильны выполнить поставленную задачу, и тогда применяют методы нестандартные, которые обычно намного сложнее стандартных (см., например, [1, с.141] или [2, с.70]). В работе [3] автором предложен алгоритм, основанный на произвольно выбранном стандартном методе, который для ОДУ во многих, если не во всех, случаях гарантирует нормальное завершение поиска искомого решения. Предлагаемый алгоритм решения ОДУ, который без особого труда распространяется и на системы ОДУ, позволяет находить решение с помощью стандартных методов как для уравнений с особенностями, так и для жестких уравнений. Ниже излагается идея предлагаемого алгоритма для численного решения одного уравнения и для систем ОДУ.

Итак, идея предлагаемого алгоритма для ОДУ заключается в том, что наряду с уравнением y'=f(x,y) (назовем это задачей 1) рассматривается родственное уравнение x'=(f(x,y))О№ (задача 2). Предлагается вместо решения одной задачи 1 поочередно решать обе задачи, при этом рассматриваемый отрезок, на котором отыскивается решение ОДУ, автоматически разбивается на множество отдельных отрезков так, что на одной их части решается задача 1, а на другой - решается задача 2. При этом задачу 1 избранным пошаговым методом решаем только при условии, когда | y' | < 1 (до тех пор, пока…), иначе переходим к решению задачи 2 и решаем ее до тех пор, пока | x' | <= 1, иначе снова переходим к задаче 1. Стартовый момент определяем путем вычисления значения производной y' в начальной точке, далее, если | y' | < 1, то решаем задачу 1, иначе - решаем задачу 2.

Для системы ОДУ , где идея предлагаемого алгоритма заключается в следующем.

, а

Помимо исходной системы уравнений (задача 1), рассматривается еще n родственных задач (систем), которые получаются из задачи 1 заменой i-го уравнения на уравнение . Таким образом, рассматривается n + 1 задача, при этом в каждый конкретный момент решается только одна из них. Какая именно - определяем следующим образом: находим значение j, для которого в “начальной” точке (начальная точка изменяется при переходе от одной задачи к другой) . Далее, если |, то решается задача 1, иначе решается (j+1)-я задача. В последнем случае начинаем с j-го уравнения, определяя (и = h), после чего для этого с остальных уравнений находим . Решение любой из этих задач продолжается до тех пор, пока либо достигается конец интересующего нас рассматриваемого отрезка, либо получаемое решение таково, что наступает необходимость смены текущей задачи на другую. Подробно об этом ниже в приводимом примере.

Для демонстрации результативности рассматриваемого алгоритма решения возьмем пример жестких уравнений из работы [1, с.140]. С начальными условиями y(0)=z(0)=1 требуется найти решение для следующей системы ОДУ:

 (1)

Величина отрезка, на котором ищется решение, здесь особой роли не играет, если не брать ее слишком малой. Можно рассмотреть отрезок [0,9], заведомо включающий всю «интересную» область. Что касается величины заданной точности , то ее можно положить равной 0.0001.

Как и в работе [1] за основу возьмем метод Эйлера, но будем применять его не стандартно, а в соответствии с идеей, изложенной в работе [3]. Описание алгоритма для решения систем ОДУ, которое здесь приводится для конкретного примера, очевидным образом распространяется на общий случай.

Итак, вместе с системой (1) рассмотрим еще две родственные системы (2) и (3):

(2)

(3)

где, очевидно, в системе (2) x'=dx/dy, а в системе (3) x'=dx/dz.

Уточним изложенную выше идею решения системы ОДУ, для чего сформулируем основные положения предложенного метода для рассматриваемой задачи:

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

2 Вместе с задачей (системой) (1) рассмотрим родственные задачи (2) и (3). При этом задачу (1) избранным пошаговым методом решаем только при условии, что и |y'|<1, и |z'|<1, иначе переходим к одной из задач (2) или (3). При этом выбор делаем в соответствии с тем, какая из величин на данный момент |y'| или |z'| окажется большей (если большей окажется |y'|, то переходим к задаче(2), иначе перейдем к задаче (3)).

3 Если была выбрана задача (2), то решаем ее до тех пор, пока |x'|<1, где x'=dx/dy, и |y'|>|z'|, иначе либо переходим к решению задачи (1), если на данный момент выполняются условия |y'|<1 и |z'|<1, либо переходим к решению задачи (3).

4 Если была выбрана задача (3), то решаем ее до тех пор, пока |x'|<1, где x'=dx/dz, и |z'|>|y'|, иначе либо переходим к решению задачи (1), если на данный момент выполняются условия |y'|<1 и |z'|<1, либо переходим к решению задачи (2).

5 Стартовый момент (какую из задач (1), (2) или (3) решать первой) определяем путем вычисления значений производных y' и z' в начальной точке, далее, если |y'|<1 и |z'|<1, то решаем задачу (1), иначе выбор делается так, как это указано в п.2.

6 В соответствии с формулировками задач (1), (2) и (3), данными выше, эти задачи решаем, только начиная со стартового момента (точнее, только одну из этих задач решаем в указанной формулировке), а дальше формулировки меняются, смысл изменений формулировок заключается в изменении начальных значений для этих задач.

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

ТЕКСТ ПРОГРАММЫ

# include <stdio.h>

# include <math.h>

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

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

float f1(float y, float z); // расчет F1(y,z)

float f1_1(float y, float z); // расчет (F1(y,z))О№

float f2(float y, float z); // расчет F2(y,z)

float f2_1(float y, float z); // расчет (F2(y,z))О№

int eiler_1( float * yy,float * zz); // расчет по методу Эйлера на // одном шаге задачи (1)

float eiler_2(float * yy,float * zz); // расчет по методу Эйлера на

// одном шаге задачи (2)

float eiler_3(float * yy,float * zz); // расчет по методу Эйлера на

// одном шаге задачи (3)

int eiler_yz(float * xx, float * yy, float * zz, float bx, float e);

// расчет по методу Эйлера задачи (1)

int eiler_xz(float * xx, float * yy, float * zz, float by, float e);

// расчет по методу Эйлера задачи (2)

int eiler_yx(float * xx, float * yy, float * zz, float bz, float e);

// расчет по методу Эйлера задачи (3)

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

int main(void)

{ int ind;

float x,y,z,x0,y0,z0,bx,by,bz,e;

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

scanf("%f%f%f%f%f%f%f",&x0,&y0,&z0,&bx,&by,&bz,&e);

// вводим: 0 1 1 5 5 5 0.0001 <Enter>

printf(" x y z\n");

x=x0; y=y0; z=z0;

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

ind=-3;

ddd: h=ind*hh/abs(ind);

switch (abs(ind))

{ case 0 : return 0;

case 1 : ind=eiler_yz(&x, &y, &z ,bx, e); break;

case 2 : ind=eiler_xz(&x, &y, &z ,by, e); break;

case 3 : ind=eiler_yx(&x, &y, &z ,bz, e); break;

default :; }

goto ddd;

}

int eiler_yz(float * xx, float * yy, float * zz, float bx, float e)

{ float x,y,y1,y2,y3,x1,z,z1,z2,z3;

char ch;

int ind;

x=*xx; y=*yy; z=*zz;

do { if (3.5f*h<=hh) h=h*3.5f;

do { y1=y; z1=z; eiler_1(&y1,&z1);

h=h/2.0f;

y2=y; z2=z; eiler_1(&y2,&z2);

y3=y2; z3=z2; eiler_1(&y3,&z3);

} while (fabs(y3-y1)>e || fabs(z3-z1)>e);

h=h*2.0f; x=x+h; y2=y; y=y3; z2=z; z=z3;

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

*xx=x; *yy=y; *zz=z;

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

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

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

if (fabs(x)>=bx) return 0;

if (fabs(y3-y2) > fabs(z3-z2)) { if ((y3-y2)<0) {ind=-2;} else {ind=2;}

return ind;} else { if ((z3-z2)<0) {ind=-3;} else {ind=3;} return ind;}

}

float f1(float y, float z)

{ return (998.0*y+1998.0*z); }

float f2(float y, float z)

{ return (-999.0*y-1999.0*z); }

float f_1(float y, float z)

{ return (1.0/(998.0*y+1998.0*z)); }

float f2_1(float y, float z)

{ return (1.0/(-999.0*y-1999.0*z)); }

int eiler_1( float * yy,float * zz)

{ float z,y;

y=*yy; z=*zz;

y=y+h*f1(y,z);

z=z+h*f2(y,z);

*yy=y; *zz=z; return 0; }

float eiler_2(float * yy,float * zz)

{ float h1,z,y;

y=*yy; z=*zz;

h1=h*f_1(y,z);

z=z+h1*f2(y,z);

*zz=z; return h1; }

float eiler_3(float * yy,float * zz)

{ float h1,z,y;

y=*yy; z=*zz;

h1=h*f2_1(y,z);

y=y+h1*f1(y,z); *yy=y; return h1; }

int eiler_xz(float * xx, float * yy, float * zz, float by, float e)

{ float x,y,y1,x1,x2,x3,z,z1,z2,z3;

char ch;

int ind;

x=*xx; y=*yy; z=*zz;

do { if (3.5f*h<=hh) h=h*3.5f;

do { z1=z; x1=x+eiler_2(&y, &z1);

h=h/2.0f; z2=z; x2=x +eiler_2(&y, &z2);

y1=y+h; z3=z2; x3=x2+eiler_2(&y1, &z3);

} while (fabs(x3-x1)>e || fabs(z3-z1)>e);

h=h*2.0f; y=y+h;

x2=x; x=x3; z2=z; z=z3;

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

*xx=x; *yy=y; *zz=z;

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

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

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

if (fabs(y)>=by) return 0;

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

{ if ((z3-z2)<0) {ind=-3;} else {ind=3;} return ind;}

}

int eiler_yx(float *xx, float *yy, float *zz, float bz, float e)

{ float x,y,y1,y2,y3,z,z1,x1,x2,x3;

char ch;

int ind;

x=*xx; y=*yy; z=*zz;

do { if (3.5f*h<=hh) h=h*3.5f;

do { y1=y; x1=x+eiler_3(&y1, &z);

h=h/2.0f; y2=y; x2=x +eiler_3(&y2, &z);

z1=z+h; y3=y2; x3=x2+eiler_3(&y3, &z1);

} while (fabs(x3-x1)>e || fabs(y3-y1)>e);

h=h*2.0f;

z=z+h;

x2=x; x=x3; y2=y; y=y3;

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

*xx=x; *yy=y; *zz=z;

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

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

} while (fabs(x3-x2)<fabs(h) && fabs(z)<bz && fabs(y3-y2)<fabs(h));

if (fabs(z)>=bz) return 0;

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

{ if ((y3-y2)<0) {ind=-2;} else {ind=2;} return ind;}

}

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

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

1 Форсайт Дж., Малькольм М., Моулер К. Машинные методы математических вычислений. - М.: Мир, 1980. - 280 с.

2 Деккер К., Вервер Я. Устойчивость методов Рунге-Кутты для жестких нелинейных дифференциальных уравнений / Пер. с англ.- М.: Мир, 1988. -332 с.

3 Литвиненко А.А. О численных методах решения обыкновенных дифференциальных уравнений // Вестник Сумского государственного университета. Серия. Технические науки. -2004. -№12(71). - С. 118-123.