Решение систем линейных дифференциальных уравнений пятиточечным методом Адамса Башфорта

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

Работу выполнил студент гр.И-29 Уханов Е.В.

Кафедра “Системы и Процессы Управления”

“ХАРЬКОВСКИЙ ПОЛИТЕХНИЧЕСКИЙ ИНСТИТУТ”

Харьков 2001

ВВЕДЕНИЕ

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

Решение такого рода задач связано с необходимостью использования численных методов , таких как : метод прогноза и коррекции , метод Адамса-Башфорта , метод Эйлера , метод Рунге-Кута , и др.  При этом , стоит задача решения системы линейных дифференциальных уравнений первого порядка одним из методов интегрирования , на произвольном промежутке времени . Одним из оптимальных методов дающих высокую точность результатов – является пяти точечный метод прогноза и коррекции Адамса-Башфорта . Для повышения точности метода используется трех точечный метод прогноза и коррекции с автоматическим выбором шага , что приводит к универсальному методу интегрирования систем дифференциальных уравнений произвольного вида на любом промежутке интегрирования  .

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

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

1. ПОСТАНОВКА ЗАДАЧИ

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

                          (1.1)



тогда как :



               А =                                                          (1.2)



 где  А заданная матрица размером  N x N .

  - вектор с N координатами , который подлежит определению ;

N – произвольное целое число ;



заданные вектора правых частей с N  координатами .

    С  использованием метода прогноза и коррекции Адамса-Башфорта пятого порядка  , необходимо получить значения  неизвестных для заданных  временных интервалов . Для  стартования  метода необходимо использовать метод прогноза и коррекции   третьего порядка с переменным шагом  , на заданных временных промежутках ..

2. МЕТОДЫ РЕШЕНИЯ

2.1. Метод прогноза и коррекции

Метод прогноза и коррекции относится к задачам класса Коши , а именно к численным решениям многошаговыми методами .

Рассмотрим задачу Коши :

            ,                                                   (2.1.1)

Подставим в (2.1.1) точное решение  y(x)  , и проинтегрируем это уравнение на отрезке   , тогда получим :

                  (2.1.2)


      где в последнем член предполагаем , что   p(x)   полином , аппроксимирующий   f(x,y(x))  . Чтобы построить этот полином , предположим , что    - приближения к решению в точках              . Будем считать для начала , что узлы Xi  расположены равномерно с шагом  h . тогда  fi  = f(xi,yi), ( i=k,k-1,k-2,…,k-N) есть приближения к   f (x,y(x))  в точках   и мы в качестве P  возьмем интерполяционный полином для выбора данных (xi,fi) , 

( i =k,k-1,k-2,…,k-N) . Таким образом , P – полином степени N , удовлетворяющий условиям  P(xi)=fi  , ( i = k,k-1,k-2,…,k-N) .  В принципе , можем проинтегрировать этот полином явно , что ведет к следующему методу  :

                                                                           (2.1.3)

         В простейшем случае , когда  N=0  , полином  P  есть константа  , равная fk  ,  и (2.1.3) превращается в обычный метод Эйлера :

                                                                                       (2.1.4)

         Если  N=1 , то P есть линейная функция , проходящая через точки  

(xk-1,fk-1)  и (xk,fk) , т.е.  

                                              (2.1.5)

     интегрируя этот полином  от Xk  до Xk+1  , получим следующий метод :

                                                                  (2.1.6)

         который является двухшаговым , поскольку использует информацию в двух точках  xk  и   xk-1  . Аналогично , если N=2  , то P  - есть кубический интерполяционный полином , а соответствующий метод определяется формулой  :

                   (2.1.7)

                Отметим , что метод (2.1.6) – есть метод Адамса-Башфорта второго порядка , (2.1.7) – метод Адамса-Башфорта четвертого порядка  .

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

                                    

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

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

	        Рассуждая также , как для метода Адамса-Башфорта , который излагается в работах : [1],[2],[3]  ,  мы мы приходим к формулам :

Прогноз :

                                                                       (2.1.8)

Коррекция :

                                                        (2.1.9)

где  h  - шаг интегрирования , изменяющийся на малом промежутке времени  в соответствии с условиями Рунге :

                     ,

где в свою очередь  - малое конкретное значение , при невыполнении условия которого увеличивается шаг  h=h*N  а - малое конкретное значение , при невыполнении условия шаг соответственно уменьшается  h=h/N , где N  - некоторое целое число больше единицы .

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

                                                  (2.1.10)

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

Метод Адамса-Башфорта использует уже посчитанные значения в точке Xk  и в предыдущих точках . В принципе , при построении интерполяционного полинома , мы можем использовать и точки Xk+1,Xk+2,… . Простейший случай при этом состаит в использовании точек Xk+1,Xk,…,Xk-N

 и построения интерполяционного полинома степени  N+1  , удовлетворяющего условиям  P(Xi)=fi , (I=k+1,k,…,k-N) .  При этом возникает класс методов , известных как методы Адамса-Моултона . Если  N=0  , то  p – линейная функция , проходящая через точки  (Xk,fk) и  (Xk+1,f k+1) , и соответствующий метод :

                                                                   (2.1.11)

является методом Адаиса-Моултона  [2]  , именно им мы воспользовались в формуле  (2.1.9) – коррекции спрогнозированной точки в трех шаговом методе . Если  N=2  , то p – кубический полином , построенный по точкам  и соответствующий метод :

                (2.1.12)

            является методом Адамса-Моултона четвертого порядка . В силу того , что по сути fk+1 – неизвестная , то методы Адамса-Моултона   (2.1.11),(2.1.12) называют неявными . В тоже время методы Адамса-Башфорта – называют явными .

Теперь воспользовавшись явной формулой  (2.1.7)  , и неявной формулой  (2.1.12)  , используя их совместно , мы приходим к методу Адамса-Башфорта четвертого порядка :

       (2.1.13)









        Стоит обратить внимание , что в целом этод метод является явным   . Сначало по формуле Адамса-Башфорта вычисляется значение , являющееся  “прогнозом”   . Затем      используется для вычисления приближенного значения  , которое в свою очередь используется в формуле Адамса-Моултона . Таким образом формула  Адамса-Моултона “корректирует” корректирует приближение , называемое формулой Адамса-Башфорта .

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

	

 где

                              A =   

Заданная матрица размером  NxN ;   - вектор с N координатами , который подлежит определению . В связи с тем , что связь между искомыми неизвестными  определяется матрицей коэффициентов A , на каждом шаге по времени , необходимо решить систему относительно неизвестных скоростей , для её решения воспользуемся модифицированным методом Гаусса , который описан в разделе 2.2  .

   Далее, интегрируя сначала  ранее описанными методами  : методом Эйлера  на первом шаге , трех точечным методом прогноза и коррекции с авто подбором шага , на малом промежутке времени и с малым начальным шагом  , для повышения точности стартующих методов на оставшемся промежутке времени производим интегрирование с постоянным шагом – пяти точечным методом прогноза и коррекции Адамса-Башфорта (2.1.13) ,  [2] , [3]  . 

2.2 Модифицированный метод Гаусса

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

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



Составить систему :       (2.2.1)



1) Каждое уравнение делиться на коэффициент   при X1



                 





 2) Теперь образуем нули в первом столбце матрицы системы : вычитаем 2-ое

из 1-ого , 3-е из 2-ого , 4-ое из 3-его :



                                                                          



                                                                                                                          (2.2.2)





3) Повторив еще раз эти операции   получим систему двух уравнений с двумя неизвестными , решение которой можно получить по формулам Крамера :



                                                                                 (2.2.3)



         Решение же X1  и  X2   можно получить , подставив в какое-либо из уравнений систем (2.2.1) и (2.2.2) и разрешив эти уравнения относительно соответствующей переменной .	

3. ОПИСАНИЕ АЛГОРИТМА

Программа начинается с вывода сообщения о программе . После происходит считывание необходимых исходных данных из файла , для дальнейшей работоспособности алгоритма , а именно – начальных условий и матрицы коэффициентов системы линейных дифференциальных уравнений первого рода , начального шага интегрирования , левого и правого условий Рунге , время интегрирования по  трех шаговому методу прогноза и коррекции ,  время интегрирования по пяти точечному методу Адамса-Башфорта .

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

Далее составляем цикл , для реализации алгоритма нахождения всех Yk+1   точек на заданном малом промежутке времени , и проверкой на условия Рунге , по трех шаговому методу прогноза и коррекции с авто подбором шага  . После чего мы организовываем цикл , реализующий алгоритм нахождения точек по методу Адамса-Башфота  , на заданном большом промежутке времени и с шагом автоматически подобранным предыдущим методом .

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

Блок-схема приведена в  Приложении 1 .

4.ОПИСАНИЕ ПРОГРАММЫ

Программа реализующая универсальный алгоритм для решения  систем линейных дифференциальных уравнений первого порядка произвольного вида , - построена по принципам  объектно-ориентированного программирования .Основная программа построена на объектной библиотеке VFH  , реализующей возможности реализации гибкого интерфейса между программой и пользователем .

Основная программа включает в себя только один модуль PACM , и использует всего два метода объекта TApplPandC   , -  метод  Application  - рабочий цикл программы  ; деструктор Done – реализует разрушение таблицы виртуальных методов  , и операций , связанных с завершением программы .

Модуль PACM включает в себя модули библиотек  - реализующих построение интерфейса  . Модуль реализующий алгоритм метода Адамса-Башфорта , и по вычесленным данным строящий график , есть – PACMBtn .

Главным родителем всех объектов есть объект – Tobject  . Основным рабочим объектом библиотеки VFH  есть объект  Tform  .  Рассмотрим потомка являющегося типичным представителем родителя  TForm  -  TApplPandC .  Он имеет два виртуалых метода  :  MouseHandler : Boolean  Б – выходным параметром которого  есть  признак закрытия формы  , и метод  FormCreate  - реализующий построение интерфейса формы  . Не виртуальный метод  Application  - предназначен для создания формы , конфигурирования программной среды , и дальнейшего управления программой .

Модуль реализующий создание и управления главного и субменю ,    есть – PACMMenu    , позволяющий пользователю изменять параметры и настройки системы , предоставляющий справку о разработчике , а также дает доступ к справочной системе  PrandCo M Help System  . Данные свойства меню реализуют объекты  TMenu  , и  THelpForm  ,  объектной библиотеки  VFH  . 

Теперь  рассмотрим модуль PACMBtn – рреализующий алгоритм построения вычисленных данных . Процедура реализующая алгоритм пяти точечного метода прогноза и коррекции Адамса-Башфорта  , -  MethodAdamsaBashforta ( h,tp,ta : real ; NU : array[1..N] of real ) – параметры которой представляют : h -  начальный шаг интегрирования ; tp – время интегрирования трех точечным методом прогноза и коррекции  , ta – время интегрирования по методу Адамса-Башфорта , NU – массив начальных условий . Данная процедура способна производить решения систем линейных дифференциальных уравнений произвольного размера , на произвольном промежутке времени интегрирования  . Вычисленные данные записываются в файлы  prandcom*.df  .  Метод реализующий алгоритм построения вычисленных данных произвольной степени сложности  , с возможностью построения графиков с не линейно изменяющимся шагом  ,  построения одновременно любого количества графиков , - есть объект TCartFile  , обладающего всеми свойствами родителей   Tform , Tchart  .

К заключению стоит заметить , что программа   PrandCo M version 2.41 -  разработана на языке Borland Pascal  под защищенный режим работы процессора и  имеет доступ ко всей оперативной памяти компьютера  . Реализует гибкий интерфейс , облегчающим работу с программным обеспечением .  Позволяет решить систему линейных дифференциальных уравнений первого порядка методом Адамса-Башфорта , с возможность просмотра результатов вычисления в виде графиков . 

Как показали тестовые программы – разработанный алгоритм предоставляет точность вычислений , погрешность которых не превышает  1% .

Тексты  программной оболочки PrandCo M  version 2.41 приведены в приложении 4 .

5. ПРИМЕРЫ РАСЧЕТОВ

Для анализа достоверности получаемых результатов рассмотрим следующие примеры :

5.1.Решение одного дифференциального уравнения

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

Пусть требуется  решить уравнение  :

                                        

при начальном условии  y(0)=1 , 0<=x<=1 ,  и шаге интегрирования  h=0.1 . Это линейное уравнение , имеющее следующее точное решение :



                  



которое поможет нам сравнить точность численного решения для случая с постоянным шагом  , т.к. точность решений  с переменным шагом выше  . Результаты расчета представлены в Таблице 1 .Как видно из таблицы, отличие между численными и аналитическими решениями удовлетворительное  даже для такого большого шага , и не превышает 2% . Теперь решим этот же пример тем же методом , но с переменным шагом . Получаем любопытные зависимости точности от выбора шага , а также шага сходимости , - которые носят периодический характер . Результаты исследования приведены в таблице 2 . Как мы видим, погрешность резко уменьшается  с использованием метода с переменным шагом , и показывает очень высокую точность решения  для численных методов , не превышающею 1% .



Таблица 2





Таблица 2

Начальный шаг

Максимальная погрешность

Сведение к шагу

 0.1

1.683 %

 0.0250

 0.01

1.163  %

 0.0100

 0.001

0.744  %

 0.0040

 0.0001

0.568  %

 0.0032

 0.00001

0.451  %

 0.0025

 0.000001

0.723  %

 0.0040

 0.0000001

0.578  %

 0.0032

 0.00000001

0.462  %

 0.0026

 0.000000001

0.740  %

 0.0041

 0.0000000001

0.592  %

 0.0033

 0.00000000001

0.473  %

 0.0026



Иллюстрация решения данного дифференциального уравнения в виде графика – приведена в Приложении 2 .

5.2.Решение системы дифференциальных уравнений

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

Рассмотрим следующую систему дифференциальных уравнений  , которую требуется решить методом Адамса-Башфорта :

                                  

		Начальными условиями здесь являются :

 .  Возьмем начальный шаг интегрирования   h=0.00001   ,  время  интегрирования по трех точечному методу прогноза и коррекции  tp=0.1  и время интегрирования по методу Адамса-Башфорта  ta=1 .

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

Иллюстрация решения данной системы дифференциальных уравнений приведены в виде графика в приложении 3 .

ЗАКЛЮЧЕНИЕ

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

Проведены тестовые расчеты  , подтвердившие высокую эффективность и точность метода Адамса-Башфорта со стартованием трех точечным методом прогноза и коррекции с  переменным шагом  .

Проведены  ряд исследований  решения систем как с постоянным шагом , так и с переменным шагом на сходимость к постоянному шагу .

Во всех случаях получены результаты высокой точности .

Список литературы 

1.Дж.Ортега , У.Пул  “Введение в численные методы решения дифференциальных уравнений ”. Пер.с англ.; под редакцией  А.А.Абрамова  - М.;Наука.Гл.ред.физ.мат.лит.1986.-288с.

2.Р.В.Хемминг “Численные методы для научных работников инженеров  ”: Пер с англ.:Под редакцией Р.С.Гутера.- Гл.ред.физ.мат.лит.1968.-203 с.

Т.Шуп.”Решение инженерных задач наЭВМ. Практическое пособие “

Пер.с англ.-М.Мир.1982.-238с.

Приложение 1 :

                       Блок схема Алгоритма





















                                                              -     

         

                                               +



















                                                               -

                                                  +











Приложение 2:









          Решение одного дифференциального уравнения









Приложение 3 :





          Решение системы линейных дифференциальных уравнений



	1-ое уравнение                                    2 –ое уравнение



                         3 – е уравнение                                     4 –ое уравнение







Приложение 4 : Тексты программ



{

 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

 |               PrandCoM version 2.41 Copiright ( c ) 2001                 |

 |                    Программа разработана студентом                       |

 |                Национального Технического Университета                   |

 |                " Харьковский Политехнический Институ "                   |

 |                            группы И - 29                                 |

 |              Кафедры Автоматического Управления Движением                |

 |                     ( Системы и процессы управления )                    |

 |                             Ухановым Е.В.                                |

 |                   NetMail ( FidoNet )   2:461/212.21                     |

 |                     E-Mail : JVUMailbox@rambler.ru                       |

 |                                                                          |

 | Программа разработана на основе объектной библиотеки  VFH version 4.XX   |

 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

}

{$M 10000,0,0}

(****************************************************************************)

(******  Дата последней разработки : 05.05.2001        **********************)

(****************************************************************************)

Program Prognoz_and_Correction_Modification;

(****************************************************************************)

 Uses PACM;

(****************************************************************************)

 var

   TPC  : TApplPandC;

(****************************************************************************)

(******************************) begin (*************************************)

                            TPC.Application;

                               TPC.Done;

(*******************************) end. (*************************************)

(****************************************************************************)







{

 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

 |                              Версия 2.XX                                 |

 |       Программа разработана студентом Национального Технического         |

 |    Университета " Харьковский Политехнический Институ " группы И - 29    |

 |       Кафедры Автоматического Управления Движением  - Ухановым Е.В.      |

 |                    NetMail ( FidoNet )   2:461/212.21                    |

 |                      E-Mail : jvumailbox@rambler.ru                      |

 |                                                                          |

 | Программа разработана на основе объектной библиотеки  VFH version 4.XX   |

 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

}

(****************************************************************************)

(****  Дата последней разработки модуля : 15.04.2001        *****************)

(****************************************************************************)



(****************************************************************************)

(*******************************) Unit PACM; (*******************************)

(****************************************************************************)

(*******************************) INTERFACE (********************************)

(****************************************************************************)

Uses FormObj,MouseObj,PACMEr,PACMMenu,PACMBtn,PACMPnl,PACMPC,PACMCnst;

(****************************************************************************)

type

    TApplPandC = object ( TForm )

    Function   MouseHandler : boolean;Virtual;

    Procedure  FormCreate;Virtual;

    Procedure  Application;

             end;

(****************************************************************************)

(******************************) IMPLEMENTATION (****************************)

(****************************************************************************)

 Procedure TApplPandC.FormCreate;

  var

    Pnl     : TPanel;

    Pnl1    : TPanel;

    TMenu1  : TCreateMenus;

   begin

     Pnl.Init(548,35,619,50,1,7,1,1,1,1,false,false);

     Pnl.Panel;

     Pnl1.Init(470,407,630,460,1,7,1,0,1,4,true,false);

     Pnl1.Panel;

     TPnl1.ToolBarCreate;

     TPnl1.PanelCreate;

     TPageControl1.PageControlCreater;

     TBitBtns.BitBtnCreaters;

     TMenu1.MenusCreate;

   end;

(********************************)

 Function TApplPandC.MouseHandler;

  var

    TMouse1     : TMouse;

    b,x,y       : word;

    TMenu1      : TCreateMenus;

    TSubMenu1   : TCreateMenus;

    ST1         : TSystemTime;

   begin

     MouseHandler:=false;

     TMouse1.GetMouseState(b,x,y);

     ST1.Init(549,36,618,49,1,15);

     ST1.SystemTime;

     TBitBtns.BitBtnHandlers(b,x,y);

     MouseHandler:=fExitBtn;

     TMenu1.MenusVisible(x,y);

     TMenu1.MenusHandlers(b,x,y);

     TPageControl1.PageControlHandlers(b,x,y);

   end;



 Procedure TApplPandC.Application;

  var

    TIEr      : TInitErrors;

   begin

     TIEr.FatalErrorVFH;

     TIEr.LoadFont('km_defj8.fnt');

     TIEr.FindImEr1('x.bi');

     InitObjGraph;

      if InitMouseJVU then

         begin

           TIEr.LfLoad('Lf.sys');

           TIEr.ErrorExec('x.bi');

           TIEr.FindFile('f1.dat');

           TIEr.FindFile('f2.dat');

           TIEr.FindFile('f3.dat');

           TIEr.FindFile('f4.dat');

           TIEr.FindFile('km_defj8.fnt');

           TIEr.FindFile('f_nfrj8.fnt');

           TIEr.FindFile('t_nfrj8.fnt');

           TIEr.FindFile('asdf.bi');

           TIEr.FindFile('pacm_n1.bi');

           TIEr.FindFile('pacm_n2.bi');

           TIEr.FindFile('pacm_n3.bi');

           TIEr.FindFile('pacm_n4.bi');

           TIEr.FindFile('PrandCoM.hlp');

           TIEr.FindFile('litj.chr');

           TIEr.FindFile('scri.chr');

           TIEr.FindFile('trip.chr');

           TIEr.FindFile('tscr.chr');

           TIEr.FindFile('initm.mtr');

           TIEr.FindFile('initnu.mtr');

        if not  fQuickHalt then

            begin

           TIEr.LoadCFG('PrandCom.cfg');

               With HT do

                  begin

                    hx1:=575;

                    hy1:=20;

                    hx2:=637;

                    hy2:=34;

                    hc:=true;

                    hs:='Закрыть';

                  end;

     Init(1,1,639,479,7,1,'Prognoz & Corrections Modifications');

     Form;

            end;

     end

       else

         begin

           TIEr.ErrorVFH;

         end;

   end;



(****************************************************************************)

(***********************************) END. (*********************************)

(****************************************************************************)