А. А. Литвиненко, доц icon

А. А. Литвиненко, доц




Скачати 74.67 Kb.
НазваА. А. Литвиненко, доц
Дата14.07.2012
Розмір74.67 Kb.
ТипДокументи

УДК 519.622


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

А.А.Литвиненко, доц.

Сумский государственный университет



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


Стандартные численные методы решения обыкновенных дифференциальных уравнений (ОДУ) и их систем в ряде случаев бессильны выполнить поставленную задачу, и тогда применяют методы нестандартные, которые обычно намного сложнее стандартных (см., например, [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

# include

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

# 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

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)
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)
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)
if (fabs(z)>=bz) return 0;

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

}

ВЫВОДЫ


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


SUMMARY


Studied are the problems which might arise while solving simple differential equations and systems (Koshi problem). There is proposed an algorithm, which nearly always guarantees normal completion of the necessary solution search.


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


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

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

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


Поступила в редакцию 12 декабря 2005 г.

Схожі:

А. А. Литвиненко, доц iconА. С. Литвиненко, О. М. Ляшенко програма та робоча програма навчальної дисципліни електричні апарати
«Електротехніка та електротехнології» спеціальності «Світлотехніка І джерела світла»./ Укл. А. С. Литвиненко, О. М. Ляшенко – Харків:...
А. А. Литвиненко, доц icon«Затверджую» Перший проректор Г. В. Стадник
Програма розроблена колективом кафедри "Електропостачання міст" у складі: д т н., к т н., доц. Абраменка І. Г., к т н., доц. Швеця...
А. А. Литвиненко, доц iconП л а н роботи лабораторії
Денисенко В. В., доц. Манько Н. В., доц. Нагрибельна І. А., доц. Гриценко І. В., доц. Полєвікова О. Б., доц. Саган О. В., ст викл....
А. А. Литвиненко, доц iconКурс п’ятий
Доц Лотоцька 436 доц. Потятинник 404 доц. Курпіль 402 доц. Бехта 423 ас. Рядська 403 доц. Федорчук 401
А. А. Литвиненко, доц icon«Затверджую» Перший проректор Г. В. Стадник
Світлотехніка І джерела світла у складі: проф., д т н. Назаренко Л. А., проф., д т н. Овчинников С. С., проф., к т н. Салтиков В....
А. А. Литвиненко, доц iconПолтавський національний педагогічний університет імені В. Г. Короленка
«молодший спеціаліст дошкільної освіти» / Укладачі: доц. А. В. Пасічніченко, доц. Н. В. Ковалевська, доц. О.І. Гришко, доц. Л. В. Зімакова,...
А. А. Литвиненко, доц iconПолтавський національний педагогічний університет імені В. Г. Короленка
«молодший спеціаліст дошкільної освіти» / Укладачі: доц. А. В. Пасічніченко, доц. Н. В. Ковалевська, доц. О.І. Гришко, доц. Л. В....
А. А. Литвиненко, доц iconФлу-41 флу-42
Науковий семінар доц. Білоус 226; доц. Захлюпана 230; доц. Ощипко 324 Сучасна українська літ мова (л) доц. Асіїв 308 сулм (п) доц....
А. А. Литвиненко, доц iconНавчальні дисципліни напряму «Психологія» що можуть бути запропоновані для вільного вибору студентів
Вовк А. О., доц. Волошок О. В., доц. Карковська Р.І., доц. Кліманська М. Б., доц. Крупська О.І., доц. Пилат Н.І, асист. Чолій С....
А. А. Литвиненко, доц iconІнститут гуманітарних І соціальних наук
Робочу навчальну програму склали: доц. Вознюк Г. Л., доц. Наконечна Г. В., доц. Микитюк О. Р
Додайте кнопку на своєму сайті:
Документи


База даних захищена авторським правом ©zavantag.com 2000-2013
При копіюванні матеріалу обов'язкове зазначення активного посилання відкритою для індексації.
звернутися до адміністрації
Документи