Комплект заданий по численным методам
я2Теоретическая часть. В данной расчетно-графической работе (далее РГР) требует-
ся составить программу для решения системы нелинейных уравне-
ний методом последовательной итерации обратной матрицы Якоби.
Суть метода в следующем: Пусть требуется решить систему нелинейных алгебраических
или трансцендентных уравнений: Fя41я0(Xя41я0,Xя42я0,...,Xя4nя0)=0; i=1,2,...,n, с начальным приближением к решению: Xя50я0=(xя41я50я0,xя42я50я0,...xя4nя50я0). Вычислительная схема реализованного метода состоит в сле-
дующем: В начале итерационного процесса матрица H полагается рав-
ной единичной: Hя50я0=E. Затем для k=0,1,... 1. Вычисляется Pя5k я0=я5 я0-я5 я0Hя5k я0*я5 я0F(Xя5kя0); 2. Находятся Xя5k+1 я0=я5 я0Xя5k я0+я5 я0tя5kя0*Pя5kя0. Первоначально tя5kя0=1. Затем путем последовательного деления - 2 -
tя5kя0 на 2 находим такое tя5kя0, чтобы выполнялось неравенство: є F(Xя5k+1я0) є < є F(Xя5kя0) є Итерационный процесс заканчивается при выполнении усло-
вия: є F(Xя5k+1я0) є < E,
где E - заданная точность. 3. Определяется Yя5kя0= F(Xя5k+1я0) - F(Xя5kя0) 4. Находится новое приближение матрицы:
Hя5k+1 я0=я5 я0Hя5k я0-я5 я0(Hя5kя0*Yя5k я0-я5 я0Pя5kя0*tя5kя0)я5 я0*я5 я0(Pя5kя0)я5T я0*я5 я0(Hя5kя0)я5T я0/я5 я0((Pя5kя0)я5T я0*я5 я0Hя5kя0*Yя5kя0) и снова повторяется вычислительный процесс с пункта 1. я2Порядок работы с программой Данная РГР представлена в виде 3 исполняемых модулей:
я1OBRJ.M, OBRF.M и FUN1.M.я0 Решением поставленной задачи занима-
ется модулья1 OBRF.Mя0, а два остальных являются вспомогательными:
я1OBRJ.M -я0 головной модуль, в котором вводятся входные данные и
выводятся результаты вычислений, ая1 FUN1.M -я0 модуль, который
пишет сам пользователь и который возвращает вычисленные левые
части для требуемого уравнения. В головной программе задаются начальные приближения, в - 3 -
виде вектора X0 а также запрашивается допустимая ошибка. Затем
вызывается модулья1 OBRJ.M,я0 который и реализует решение данной
системы уравнений методом последовательной итерации обратной
матрицы Якоби. Внутри себя данный модуль по мере необходимости
вызывает функциюя1 FUN1.Mя0, которую пишет сам пользователь. я2Описание работы программ В связи с тем, что данная РГР состоит из 3 частей, то
опишем их по одиночке (распечатки данных модулей приведены в
приложении): 1.я1 OBRJ.M Головной модуль Входные данные: отсутствуют. Выходные данные: отсутствуют. Язык реализации: PC MathLab. Операционная система: MS-DOS 3.30 or Higher. Пояснения к тексту модуля: "Стандартный" головной модуль. В данном модуле задаются
начальные значения в виде вектора, например: Xя40я0=[0.4 0.9] Также в данном модуле запрашивается допустимая ошибка,
очищается экран, а также производятся другие подготовительные
действия. Затем происходит вызов модуляя1 OBRF.Mя0 с полученными вход-
ными данными. Формат вызова данного модуля описан далее (в
описании самого модуля). - 4 - После вычислений в головную программу возвращаются ре-
зультаты вычислений на основе которых строятся графики а также
выводятся оценки по затратам машинного времени и быстродейс-
твия. 2.я1 OBRF.M Вычислительный модуль Входные данные: FunFcn - имя функции, написанной пользователем, которая
вычисляет левые части для требуемой системы в определенной
точке. X0 - вектор-строка, определяющий начальные значения (на-
чальное приближение). E - допустимая ошибка. Выходные данные: Tout - Столбец итераций ("Время") Xout - Столбцы значений вычисленных на каждом этапе для
каждой итерации DXout - Столбцы погрешностей по каждой компоненте, вычис-
ленные на определенном этапе Язык реализации: PC MathLab Операционная система: MS-DOS 3.30 or Higher Пояснения к тексту модуля: Данный "вычислительный" модуль реализует метод последова- - 5 -
тельной итерации обратной матрицы Якоби. Общая структура вызо-
ва данного модуля: [Tя4outя0,Xя4outя0,DXя4outя0]=OBRF(FunFcn,Xя40я0,E); Значения каждого из параметров были описаны выше. На начальном этапе в данном модуле инициализируются внут-
ренние переменные (например, задается единичная матрица H, в
соответствии с размерностью X0), формируются (на основе на-
чальных значений) первичные элементы матриц Tout,Xout,DXout.
Затем данная функция, как и многие другие в численных методах,
имеет вид: While ОШИБКА > ДОПУСТИМОЙ ОШИБКИ Оператор1 Оператор2 ......... ......... ОператорN End Внутри данного цикла происходят вычисления внутренней пе-
ременной Pя5kя0 на каждом шаге K и, вычисляется начальное прибли-
жение Xя5k+1я0. Первоначально t=1 (Не номер итерации, а внутренний
параметр!). Затем, в очередном цикле While...End в случае, ес-
ли єF(Xя5k+1я0)є < єF(Xя5kя0)є t=t/2 и снова вычисляется Xя5k+1я0. Когда
очередное Xя5k+1я0 найдено, вычисляется Yя5kя0, а затем и новое приб-
лижение матрицы H. Итерационный процесс заканчивается, если - 6 -
F(Xя5k+1я0)є < E. Если данное условие не выполняется - итерацион-
ный процесс продолжается заново. Формирование выходных значений-матриц происходит внутри
данного цикла и поэтому никаких дополнительных действий не
требуется, то есть с окончанием данного цикла заканчивается и
сама функция. 3.я1 FUN1.M Модуль, вычисляющий левые части Входные данные: X - вектор-строка, задающий точки для вычислений по каж-
дой компоненте. Выходные данные: FF - вектор-строка, возвращающий значения каждой компо-
ненты в определенной точке Язык реализации: PC MathLab Операционная система: MS-DOS 3.30 or Higher Пояснения к тексту модуля: В принципе, текст данного модуля не требует пояснений. В
нем пользователь реализует систему уравнений, которая подлежит
решению. То есть на входные значения X данная функция возвра-
щает левые части по каждому уравнению. Единственное требование
к данному модулю - соблюдение формата, то есть входные и вы-
ходные данные должны быть представлены в виде вектор-строк. я2Сравнительный анализ и я2оценка быстродействия. - 7 - Сравнительный анализ показал, что данный метод обладает
неплохой сходимостью, так как попробованный метод простой ите-
рации с параметром вообще отказался сходиться для данной сис-
темы. Однако хорошо подходит для сравнения дискретный метод
Ньютона, так как данные методы практически одинаковы что по
точности что по затратам. 1. Метод последовательной итерации обратной матрицы Якоби Число операций: порядка 682 Быстродействие: порядка 0.11 секунды 2. Метод Ньютона дискретный Число операций: порядка 990 Быстродействие: порядка 0.22 секунды Как видно из вышеприведенных данных, эти два метода очень
близки между собой, но метод Ньютона дискретный более сложен в
реализации, однако обладает лучшей сходимостью, например при
начальных значениях Xя50я0=[2.0 2.0]; метод последовательной ите-
рации обратной матрицы Якоби уже не справляется, в то время
как дискретный метод Ньютона продолжает неплохо работать. Од-
нако метод Ньютона требует больших затрат машинного времени и
поэтому при выборе метода необходимо исходить их конкретных
условий задачи и если известно довольно точное приближение и
требуется быстрота вычислений, то к таким условиям отлично
подходит разработанный метод последовательной итерации обрат-
ной матрицы Якоби. - 8 - я2Выводы В данной РГР был разработан и реализован метод последова-
тельной итерации обратной матрицы Якоби, предназначенный для
решения системы нелинейных уравнений. Программа, реализованная
на языке PC MathLab хотя и не является оптимальной, однако вы-
полняет поставленную задачу и решает системы уравнений. Реали-
зованный метод не отличается повышенной сходимостью и требует
довольно точного начального приближения, однако довольно быст-
ро сходится к точному решению, то есть его можно порекомендо-
вать для вычисления непростых систем нелинейных уравнений при
наличии довольно точного начального приближения и наличия вре-
менных ограничений. я2Список литературы 1. О.М.Сарычева. "Численные методы в экономике. Конспект
лекций", Новосибирский государственный технический универси-
тет, Новосибирск 1995г. 2. Д.Мак-Кракен, У.Дорн. "Численные методы и программиро-
вание на Фортране", Издательство "Мир", М. 1977г. 3. Н.С.Бахвалов. "Численные методы", Издательство "Нау-
ка", М. 1975г. я2Теоретическая часть Программа предназначена для численного решения системы
обыкновенных дифференциальных уравнений вида: Y'=F(X,Y), с начальными условиями Y(Xя40я0)=Yя4oя0 на отрезке
[X,X] методом Хемминга с постоянным шагом интегрирования. В
каждой i+1 точке находим начальное приближение Р к решению Y
по предсказывающей формуле: Pя4i+1я0=Yя4i-3я0+4*h*(2*Y'я4iя0-Y'я4i-1я0+2*Y'я4i-2я0)/3, где Yя4i-3я0 решение в i-3 точке, Y'я4iя0,Y'я4i-1я0,Y'я4i-2я0 - значения производных в точках i, i-1,
i-2 соответственно. Для улучшения решения используется корректирующая формула Cя4i+1я0=[9*Yя4iя0-Yя4i-2я0+3*h*(M'я4i+1я0+2*Y'я4iя0-Y'я4i-1я0)]/8, где Mя4i+1я0=Pя4i+1я0-112*(P-Cя4iя0)/121; M'я4i+1я0=F(Xя4i+1я0,Mя4i+1я0). Решение системы в i+1 точке находится по формуле G=Wя4jя0*|Pя4i+1,jя0-Cя4i+1я0|, где Wя4jя0=1 j- номер компоненты вектора. На участке "разгона" значения Yя4i-kя0 и Y'я4i-kя0 (k=0, 1, 2)
вычисляются методом Рунге-Кутта по формуле Yя4iя0=Uя4iя5(2)я0-(Uя4iя5(i)я0-Uя4iя5(2)я0)/15, где i- номер точки, в которой
ищется решение, Uя4iя0- решение системы в i-ой точке, полученное с
шагом h/l; Uя5(l)я4i-m+1/lя0=Aя5(l)я4i-mя0+(Kя5(l)я41я0+2*Kя5(l)я42я0+2*Kя5(l)я43я0+Kя5(l)я44я0)я4i--m+1/lя0/6я4, где j=1, 2, ..., n, Kя5(l)я41я0=h*F(Xя4i-mя0,Aя5(l)я4i-mя0)/l; Kя5(l)я42я0=h*F(Xя4i-mя0+h/2*l,Aя5(l)я4i-mя0+Kя5(l)я41я0/2)/l; Kя5(l)я43я0=h*F(Xя4i-mя0+h/2*l,Aя5(l)я4i-mя0+Kя5(l)я42я0/2)/l; Kя5(l)я44я0=h*F(Xя4i-mя0+h/l,Aя5(l)я4i-mя0+Kя5(l)я43я0)/l. A, U ,K - векторы n-го порядка; l=1, 2; m=1 при l=1; m=1,
1/2 при l=2; Aя5(l)я4i-1я0=Yя5(l)я4i-1я0; Aя5(2)я4i-1/2я0=Uя5(2)я4i-1/2я0. я2Характеристика программы.я0 Программа состоит из стандартной информативы, реализующей
описанный метод, рабочей информативы, задающей правые части
уравнений системы и директивы. Длина стандартной информативы 1600 символов. Объем исход-
ных данных : 7 чисел, 2 массива, n функций. В результате рабо-
ты программы на печать выводится на участке "разгона" X, зна-
чения функций и производных, далее X, G и Y[n] на всем отрезке
интегрирования через Ю шагов и в конце отрезка. Программа рекомендуется для решения систем обыкновенных
дифференциальных уравнений на больших отрезках, так как счита-
ет быстрее одноточечных методов. Для контроля постоянно выво-
дится погрешность вычислений G, которая позволяет следить за
точностью решения. "Разгон" (нахождение значений функций и производных в
точках Xя40я0, Xя40я0+Q, Xя40я0+2*Q , Xя40я0+3*Q, где Q - шаг интегрирования )
осуществляется методом Рунге-Кутта с увеличенной разрядностью. В программе предусмотрена возможность при получении боль-
шой погрешности вычисления в точка "разгона" уменьшить шаг ин-
тегрирования в этих точках (см. способ задания J), а при быст-
ром возрастании погрешности вычислений G уменьшить шаг интег-
рирования методом Хемминга или увеличить разрядность вычисле-
ний. Программа позволяет производить интегрирование как с по-
ложительным, так и с отрицательным шагом (соответственно меня-
ются Xя40я0, Xя4кя0 и Q). я2Порядок решения задачи.я0 Для решения задачи вводятся стандартная и рабочая инфор-
мативы и директива. В рабочей информативе после метки Ц программа вычисления
правых частей системы. Здесь 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. я2Описание работы программы Данная расчетно-графическая работа (далее РГР) составлена
на языке PC MathLab ( PC-MATLAB (c) Copyright The MathWorks,
Inc. 1984-1989 Version 3.5f 17-July-89 Serial Number 22961) и
выполнена в виде двух модулей (третий - контрольный пример),
распечатка которых приведена в приложении. 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.я1 я_я0Hem.m Модуль, которые непосредственно и решает систему ОДУ ме-
тодом Хемминга. Входные данные: FunFcn - имя функции, вычисляющей левые части X0,XK - начальное и конечное значение для счета Q - шаг интегрирования J - число, определяющее во сколько раз нужно уменьшать
шаг интегрирования методом Рунге-Кутта (далее Р-К) на участке
"разгона" для получения того же порядка точности, что и в ме-
тоде Хемминга N - порядок системы Y - вектор начальных значений W - вектор коэффициентов для вычисления невязки UR - число, определяющее период печати Выходные данные: x - матрица точек, для которых вычислено решение y - матрица решений dg - ошибка интегрирования Язык реализации: PC MathLab. Операционная система: MS-DOS 3.30 or higher Пояснения к тексту модуля: Данный модуль содержит в своем теле всего одну функцию,
входные и выходные данные которой являются входными и выходны-
ми данными текущего модуля. Они описаны выше. Мы же займемся
описанием данной функции: После описания функции HEM устанавливается формат выход-
ных данных LONG E, а также происходит инициализация рабочих
массивов, как массивы значений функции в точках i-3, i-2, i-1;
массивы значений производных в этих же точках, массивы правых
частей и т.д. Всвязи с отсутствием в языке MathLab конструкции
безусловного перехода, используется конструкции While 1 (бес-
конечный цикл), Break (переход к началу While) и IF (Если).
Из-за таких немного "странных" конструкций вся дальнейшая
часть программы может быть весьма условно представлена такой
схемой: While (не конец расчетов) While 1 ... IF ... ... END ... ... IF ... ... ... END ... end end В целом, в данных циклах вычисляется номер шага, и если
он меньше 5, то вычисления сводятся к вычислению методом Р-К,
а если больше 5, то производятся вычисления по методу Хемминга
со всеми своими дополнительными действиями, как вычисление по
корректирующей формуле и т.д. В обоих случаях происходит об-
новление рабочих и других промежуточных массивов и вывод ин-
формации на экран. В случае решения в точках "разгона" вычис-
ляются также коэффициенты K1, K2, K3, K4, используемые в мето-
де Р-К. Также функция сама проверяет точность вычислений и в
случае необходимости корректирует шаг. Если шаг "сделан", то
программа выводит результаты на экран и заносит их в массив,
который представлен в виде нескольких столбцов. Также в необ-
ходимых для функции случаях она обращается к подпрограмме
FunFcn, которая занимается вычислением левых частей, вызов и
возврат значений которой должен быть следующим: Z=feval(FunFcn,x,y), где Z - вектор вычисленных левых частей, X,Y - векторы точек, для которых производится вычисление. Для удобства отладки и описания, программа разбита на
части, обозначенные русскими заглавными буквами (Ш,Щ,Л и
т.д.), которые соответствуют блокам, обозначенным в примере
программы, приведенной в задании. Несмотря на то, что приве-
денная программа написана на условном языке, прокомментировать
текст нашей программы на языке MathLab довольно удобно с ис-
пользованием данных обозначений (Конечно, часть блоков опуще-
на, в связи с отсутствием принципиальной значимости. Кроме то-
го изменен порядок появления блоков в программе): "Э" - начальное вычисление левых частей. "Ф" - общий цикл, в котором и происходят собственно все
вычисления. Он начинается с конструкций: While (xk-x1)>0 While 1, то есть пока не будет достигнут конец, все вычисления
происходят внутри этого цикла. Также внутри блока "Ф" происходят вычисления корректирую-
щей формулы IAR(i) а также оценка погрешности вычислений G,
если они еще не были рассчитаны на предыдущем шаге. "Ц" - вычисление левых частей. "Щ" - на этом этапе происходит перемещение данных в рабо-
чих массивах и X=X+H, то есть происходит переход к следующему
шагу. Также на этом этапе происходит вывод на экран и формиро-
вание выходных массивов Yout, Xout, DGout. "Л" - в этом блоке происходит расчет самой предсказываю-
щей формулы метода Хемминга - P(i) и Y(i) и происходит расчет
левых частей, т.е. снова в программе появляется блок "Ц". Затем опять продолжается блок "Ф". Идет проверка на каком
шаге мы находимся и, если он (шаг) меньше 5, то идет подготов-
ка к расчету методом Р-К, то есть идет участок "разгона". Со-
ответственно идет расчет коэффициентов K1, K2, K3 и K4, необ-
ходимых для метода Р-К. Также внутри данного блока еще раз
встречается блок "Ц", в котором происходит расчет левых час-
тей, необходимых для метода Р-К. Далее происходит проверка шага и уменьшение или увеличе-
ние его в соответствии с заданной точностью. Для возможности
"отката" в случае большого или маленького шага используется
переменная Х1. Также еще раз встречается блок "Ц". Затем, в
случае если все коэффициенты К1-К4 вычислены и шаг удовлетво-
ряет требованиям точности, то происходит расчет методом Р-К,
а также, естественно происходит формирование выходных массивов
Yout, Xout и DGout, а также происходит переход к следующему
шагу (то есть X=X+H) и переход к блоку "Э". На этом кончается блок "Ф" и вся функция. В начале блока
"Ф" происходит проверка на конец вычислений и если расчеты за-
кончились, то есть мы достигли Xk то происходит возврат в го-
ловную программу. Все выходные данные формируются внутри блока
"Ф", поэтому никаких дополнительных действий не производится. я2Сравнительный анализ и оценка быстродействия Для сравнения полученных результатов с другими методами
используется метод Адамса, разработанный другой бригадой. Число операций в методе Хемминга: порядка 2200. Быстродействие: порядка 0,8 секунды. Число операций в методе Адамса: порядка 560. Быстродействие: порядка 0,55 секунды. (Вычисления проводились на компьютере i486DX4-100) Как видно из вышеприведенных данных, метод Хемминга проигрывает
по временным показателям и по затратам машинного времени методу Адам-
са, однако стоит заглянуть в приложение, где приведены распечатки гра-
фиков решений и ошибок обоих методов и сразу видно, что метод Адамса
не справляется с контрольным примером для нашей системы, так как ошиб-
ка у него к концу вычислений (Xk=1) возрастает, а в "нашем" методе -
стремится к 0. я2Выводы Данная РГР по предмету "Численные методы в экономике" ре-
ализует метод Хемминга, который предназначен для решения зада-
чи Коши для обыкновенных дифференциальных уравнений. Програм-
ма, написанная на языке MathLab, хотя и не является оптималь-
ной, но решает поставленную задачу и решает ОДУ довольно боль-
ших степеней сложности, с которыми не справляются другие мето-
ды (например метод Адамса). Это связано с тем, что метод Хем-
минга является многоточечным, а в связи с этим и повышается
точность вычислений, а также устойчивость метода. Однако дан-
ный метод требует больших вычислительных затрат, что связано с
довольно громоздкими формулами а также с большим объемом вы-
числений и поэтому для относительно простых систем целесооб-
разно использовать более простые методы решения. я2Список литературы 1. Д.Мак-Кракен, У.Дорн. "Численные методы и программиро-
вание на Фортране", Издательство "Мир", М. 1977г. 2. О.М.Сарычева. "Численные методы в экономике. Конспект
лекций", Новосибирский государственный технический универси-
тет, Новосибирск 1995г. 3. Н.С.Бахвалов. "Численные методы", Издательство "Нау-
ка", М. 1975г. я2ВВЕДЕНИЕ В экономике очень часто используется модель, называемая "черный
ящик", то есть система у которой известны входы и выходы, а то, что
происходит внутри - неизвестно. Законы по которым происходят изменения
выходных сигналов в зависимости от входных могут быть различными, в
том числе это могут быть и дифференциальные законы. Тогда встает проб-
лема решения систем дифференциальных уравнений, которые в зависимости
от своей структуры могут быть решены различными методами. Точное реше-
ние получить очень часто не удается, поэтому мы рассмотрим численные
методы решения таких систем. Далее мы представим две группы методов:
явные и неявные. Для решения систем ОДУ неявными методами придется ре-
шать системы нелинейных уравнений, поэтому придется ввести в рассмот-
рение группу методов решения систем нелинейных уравнений, которые в
свою очередь будут представлены двумя методами. Далее следуют теорети-
ческие аспекты описанных методов, а затем будут представлены описания
программ. Сами программы, а также их графики приведены в приложении. Также стоит отметить, что в принципе все численные методы так или
иначе сводятся к матричной алгебре, а в экономических задачах очень
часто матрицы имеют слабую заполненность и большие размеры и поэтому
неэффективно работать с полными матрицами. Одна из технологий, позво-
ляющая разрешить данную проблему - технология разреженных матриц.
В связи с этим, мы рассмотрим данную технологию и операции умножения и
транспонирования над такими матрицами. Таким образом мы рассмотрим весь спектр лабораторных работ. Опи-
сания всех программ приводятся после теоретической части. Все тексты
программ и распечатки графиков приведены в приложении. я2ТЕОРЕТИЧЕСКАЯ ЧАСТЬ 1. ОПИСАНИЕ МЕТОДОВ ИНТЕГРИРОВАНИЯ СИСТЕМ ОДУ И ИХ ХАРАКТЕРИСТИК ЯВНЫЙ МЕТОД ЭЙЛЕРА И ЕГО ХАРАКТЕРИСТИКИ Алгоритм этого метода определяется формулой: xя5m+1я0 = xя5mя0 + hя4mя0*F(xя5mя0, tя4mя0)я4,я0 (1) которая получается путём аппроксимации ряда Тейлора до членов перво- го порядка производной x'(tя4mя0),т.к. порядок точности равен 1 (s=1). Для аналитического исследования свойств метода Эйлера линеари- зуется исходная система ОДУ x' = F(x, t) в точке (xя5mя0,tя4mя0): x' = A*x, (2) где матрица А зависит от точки линеаризации (xя5mя0,tя4mя0). Входной сигнал при линеаризации является неизвестной функцией времени и при фиксированном tя4mя0 на шаге hя4mя0 может считаться констан- той. Ввиду того ,что для линейной системы свойство устойчивости за- висит лишь от А, то входной сигнал в системе (2) не показан. Элемен- ты матрицы А меняются с изменением точки линеаризации,т.е. с измене- нием m. Характеристики метода: 1. я_Численная устойчивостья.. Приведя матрицу А к диагональной форме : А = Р*я7lя0*Ря5-1я0 и перейдя к новым переменным y = Pя5-1я0*x , система (2) примет вид : y' = я7lя0*y (3) Тогда метод Эйлера для уравнения (3) будет иметь вид : yя5m+1я0 = yя5mя0 + h*я7lя0 * yя5mя0, (4) где величина h*я7lя0 изменяется от шага к шагу. В этом случае характеристическое уравнение для дискретной сис- темы (4) будет иметь вид : 1 + h*я7lя0 - r = 0. А корень характеристического уравнения равен: r = 1+ h*я7lя0, где я7lя0 =я7 aя0 я_+я.я7 bя0 . я_Случай 1я.. Исходная система (3) является асимптотически устой- чивой , т.е. нулевое состояние равновесия системы асимптотически ус- тойчиво и я7 aя0 < 0. Областью абсолютной устойчивости метода интегрирования системы (4) будет круг радиусом , равным 1 , и с центром в точке (0, -1). Таким образом , шаг h должен на каждом интервале интегрирования под- бираться таким образом, чтобы при этом не покидать область А . Но в таком случае должно выполняться требование : h < 2*я7tя0, (5) где я7tя0=1/я72a2я0 - постоянная времени системы (3) . Она определяет ско- рость затухания переходных процессов в ней. А время переходного про- цесса Tя4ппя0 = 3*я7tя0 , где я7tя0 = я72a2я5-1я0. Если иметь ввиду, что система (3) n-го порядка, то Tя4ппя0 > 3*я7tя4maxя0, где я7tя4maxя0 = я72aя4minя72я5-1я7 я0; я7aя4min я0= min я7aя4iя0 . Из всех я7aя4iя0 (i=[1;n]) выбирает- ся максимальное значение : maxя72aя4iя72я0=я7aя4maxя0 и тогда можно вычислить : я7tя4min я0= 1/я7aя4maxя0, а шаг должен выбираться следующим образом : h < 2/я7aя4maxя0 или h < 2*я7tя4minя0. я_Случай 2я.. Нулевое состояние равновесия системы (2) неустойчи- во, т.е. я7aя0 > 0 . Т.е. система тоже неустойчива , т.е. я72я0rя72я0>1. Поэтому областью относительной устойчивости явного метода Эйлера является вся правая полуплоскость . Выполняется требование : я72я0 1+h*я7lя0 я72 я0
|
