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
Областью абсолютной устойчивости метода интегрирования системы (4) будет круг радиусом , равным 1 , и с центром в точке (0, -1). Таким образом , шаг h должен на каждом интервале интегрирования под бираться таким образом, чтобы при этом не покидать область А . Но в таком случае должно выполняться требование : h
где 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 ... Нулевое состояние равновесия системы (2) неустойчи во, т. е. 7a 0 > 0 . Т. е. система тоже неустойчива , т. е. 72 0r 72 0>1. Поэтому
областью относительной устойчивости явного метода Эйлера является вся правая полуплоскость . Выполняется требование : 72 0 1+h* 7l 0 72 0 2. _Точность метода ...
Так как формула интегрирования (1) аппроксимирует ряд Тейлора для функции x(t 4m 0+1) до линейного по h члена включительно. Существует такое значение t в интервале [t 4m 0 , t 4m+1 0], при котором Е 4i 5am 0 =1/2! * h 4m 52 0*x 4i 0''(t), где i=[1; n]. 3. _Выбор шага интегрирования ...
Должны соблюдаться условия абсолютной (5) или относительной (6) устойчивости в зависимости от характера интегрируемой системы. Для уравнения первого порядка шаг должен быть : h Для уравнений n-ого порядка : h 4i 0 А окончательно шаг выбирают по условиям устойчивости : h
Вначале задаётся допустимая ошибка аппроксимации , а в процессе ин тегрирования шаг подбирается следующим образом :
1) по формуле (1) определяется очередное значение x 5m+1 0; 2) определяется dx 4i 5m 0 = x 4i 5m+1 0 - x 4i 5m 0 ; 3) условие соблюдения точности имеет вид :
h 4i 5m 0
4) окончательно на m-м интервале времени выбирается в виде: h 4m 0 = min h 4i 5m 0. ЯВНЫЕ МЕТОДЫ РУНГЕ-КУТТА
Метод Эйлера является методом Рунге-Кутта 1-го порядка . Методы Рунге-Кутта 2-го и 4-го порядка являются одношаговыми , согласуются с рядом Тейлора до порядка точности s , который равен порядку метода . Эти методы не требуют вычисления производных функций , а только самой функции в нескольких точках на шаге h 4m 0. Алгоритм метода Рунге-Кутта 2-го порядка состоит в следующем : x 5m+1 0 = x 5m 0 + h 4m 0/2 (k 41 0+k 42 0),
где k 41 0=f(x 5m 0, t 4m 0) ; k 42 0=f(x 5m 0+h 4m 0*k 41 0, t 4m 0+h 4m 0). Ошибка аппроксимации Е 5a 0 = k*h 4m 53 0 . Алгоритм метода Рунге-Кутта 4-го порядка x 5m+1 0=x 5m 0+h 4m 0/6(k 41 0+2k 42 0+2k 43 0+k 44 0),
где k 41 0=f(x 5m 0, t 4m 0); k 42 0=f(x 5m 0+h 4m 0/2*k 41 0, t 4m 0+h 4m 0/2); k 43 0=f(x 5m 0+h 4m 0/2*k 42 0, t 4m 0+h 4m 0/2); k 44 0=f(x 5m 0+h 4m 0*k 43 0, t 4m 0+h 4m 0). Ошибка аппроксимации Е 5a 0 = k*h 4m 55 0. НЕЯВНЫЙ МЕТОД ЭЙЛЕРА И ЕГО ХАРАКТЕРИСТИКИ
Неявный метод Эйлера используется для интегрирования " жест ких " систем. "Жесткие" системы это такие системы, в которых 7 0 ( 7l 4max 0) и ( 7l 4min 0) сильно отключаются друг от друга , то в решениях системы x' = A*x (1)
будут присутствовать экспоненты, сильно отличаются друг от друга по скорости затухания . Шаг интегрирования для таких систем должен вы бираться по условиям устойчивости из неравенства h
где 7t 0=1/ 72a2 0 - постоянная времени системы y' = 7l 0*y . Она определяет скорость затухания переходных процессов в ней . Неравенство (2) должно выполняться на всем участке решения , что соответственно тре бует значительных затрат машинного времени. Алгоритм этого метода определяется формулой: x 5m+1 0 = x 5m 0 + h 4m 0*F(x 5m+1 0, t 4m+1 0) 4 0(3)
Если h 4m 0 мал, то x 5m 0 и х 5m+1 0 близки друг к другу. В качестве на чального приближения берется точка x 5m 0 , а следовательно , между x 5m 0 и x 5m+1 0 будет существовать итерационный процесс.
Для аналитического исследования свойств метода Эйлера линеа ризуется исходная система ОДУ x' = F(x, t) в точке (x 5m 0, t 4m 0): x' = A*x,
где матрица А зависит от точки линеаризации (x 5m 0, t 4m 0). Входной сигнал при линеаризации является неизвестной функцией времени и при фиксированном t 4m 0 на шаге h 4m 0 может считаться констан той. Ввиду того , что для линейной системы свойство устойчивости за висит лишь от А, то входной сигнал в системе (1) не показан. Элемен ты матрицы А меняются с изменением точки линеаризации, т. е. с измене нием m. Характеристики метода: 1. _Численная устойчивость ...
Приведя матрицу А к диагональной форме : А = Р* 7l 0*Р 5-1 0 и перейдя к новым переменным y = P 5-1 0*x , система (3) примет вид : y' = 7l 0*y (4) Тогда метод Эйлера для уравнения (4) будет иметь вид : y 5m+1 0 = y 5m 0 + h* 7l 0 * y 5m+1 0, (5) где величина h* 7l 0 изменяется от шага к шагу.
В этом случае характеристическое уравнение для дискретной сис темы (5) будет иметь вид : 1 - h* 7l 0*r - 1 = 0. А корень характеристического уравнения равен: r = 1/(1-h* 7l 0) , где 7l 0 = 7 a 0 _+ . 7 b 0 .
_Случай 1 ... Исходная система (4) является асимптотически устой чивой , т. е. нулевое состояние равновесия системы асимптотически ус тойчиво и 7 a 0
Областью абсолютной устойчивости метода интегрирования системы (5) будет вся левая полуплоскость. Таким образом , шаг h должен на каждом интервале интегрирования подбираться таким образом, чтобы при этом не покидать эту область. Но в таком случае должно выполняться требование : h
где 7t 0=1/ 72a2 0 - постоянная времени системы (4) . Она определяет ско рость затухания переходных процессов в ней. А время переходного про цесса T 4пп 0 = 3* 7t 0 , где 7t 0 = 72a2 5-1 0. Если иметь ввиду, что система (4) 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 ... Нулевое состояние равновесия системы (4) неустойчи во, т. е. 7a 0 > 0 . Т. е. система тоже неустойчива , т. е. 72 0r 72 0>1, а следовательно : 72 0 1/(1-h* 7l 0) 72 0 > 1.
С учетом ограничения на скорость изменения приближенного ре шения относительно точного : 72 0 1/(1-h* 7l 0) 72 0 > 7 2 0e 5hl 7 2 0. (7)
Из этого соотношения следует , что при 72 0h* 7l2 0 -> 1 левая часть стремится к бесконечности . Это означает , что в правой полуплоскос ти есть некоторый круг радиусом , равным 1 , и с центром в точке (0, 1), где неравенство (7) не выполняется. 2. _Точность метода ...
Ошибка аппроксимации по величине равна ошибке аппроксимации явного метода Эйлера , но она противоположна по знаку : Е 4i 5am 0 =-1/2! * h 4m 52 0*x 4i 0''(t), где t 4m 0 i=[1; n]. 3. _Выбор шага интегрирования ...
Должны соблюдаться условия абсолютной (6) или относительной (7) устойчивости в зависимости от характера интегрируемой системы. Для уравнения первого порядка шаг должен быть : h Для уравнений n-ого порядка : h 4i 0
Однако область абсолютной устойчивости - вся левая полуплос кость . Поэтому шаг с этой точки зрения может быть любым.
Условия относительной устойчивости жестче, так как есть об ласть , где они могут быть нарушены. Эти условия будут выполняться , если в процессе решения задачи контролировать ошибку аппроксимации , а шаг корректировать . Таким образом, шаг можно выбирать из условий точности, при этом условия устойчивости будут соблюдены автоматичес ки. Сначала задается допустимая погрешность аппроксимации : E 4i 5aдоп 0 где i=[1; n].
Шаг выбирается в процессе интегрирования следующим образом: 1. Решая систему (3) относительно x 5m+1 0 с шагом h 4m 0, получаем очередную точку решения системы x = F(x, t) ; 2. Оцениваем величину ошибки аппроксимации по формуле:
Е 4i 5am 0 = 72 0h 4m 7/ 0(h 4m 0+h 4m-1 0)*[(x 4i 5m+1 0 - x 4i 5m 0) h 4m 7/ 0h 4m-1 0*(x 4i 5m 0 -x 4i 5m-1 0)] 72
3. Действительное значение аппроксимации сравнивается с до пустимым. Если Е 4i 5am 0
тивном случае шаг уменьшается в два раза , и вычисления повторяются с п. 1. 4. Рассчитываем следующий шаг: h 4i 5m+1 0 = SQR(E 4i 5aдоп 7/2 0Е 4i 5am 72 0) * h 4m 0.
5. Шаг выбирается одинаковым для всех элементов вектора X : h 4m+1 0 = min h 4i 5m+1 0.
6. Вычисляется новый момент времени и алгоритм повторяется . НЕЯВНЫЕ МЕТОДЫ РУНГЕ-КУТТА
Метод Эйлера является методом Рунге-Кутта 1-го порядка . Методы Рунге-Кутта 2-го и 4-го порядка являются одношаговыми , согласуются с рядом Тейлора до порядка точности s , который равен порядку метода . Эти методы не требуют вычисления производных функций , а только самой функции в нескольких точках на шаге h 4m 0. Алгоритм метода Рунге-Кутта 2-го порядка состоит в следующем: x 5m+1 0 = x 5m 0 + h 4m 0/2 (k 41 5m+1 0+k 42 5m+1 0), где k 41 5m+1 0=f(x 5m+1 0, t 4m+1 0); k 42 5m+1 0=f(x 5m+1 0-h 4m 0*k 41 5m+1 0, t 4m+1 0). Ошибка аппроксимации Е 4m 5a 0 = k*h 4m 53 0 . Алгоритм метода Рунге-Кутта 4-го порядка
x 5m+1 0 = x 5m 0 + h 4m 0/6 (k 41 5m+1 0+2k 42 5m+1 0+2k 43 5m+1 0+k 44 5m+1 0), где k 41 0=f(x 5m+1 0, t 4m+1 0); k 42 0=f(x 5m+1 0-h 4m 0/2*k 41 5m+1 0, t 4m+1 0-h 4m 0/2); k 43 0=f(x 5m+1 0-h 4m 0/2*k 42 5m+1 0, t 4m+1 0-h 4m 0/2); k 44 0=f(x 5m+1 0-h 4m 0*k 43 5m+1 0, t 4m 0). Ошибка аппроксимации Е 5a 0 = k*h 4m 55 0. 2. МЕТОДЫ РЕШЕНИЯ НЕЛИНЕЙНЫХ САУ МЕТОД НЬЮТОНА Итерационная формула метода Ньютона имеет вид
X 5m+1 0=X 5m 0- 5 0J 5-1 0* 5 0(X 5m 0) 5 0* 5 0F(X 5m 0), где J(X)=F 4x 5| 0/ 4x=xm Характеристики метода: 1. Сходимость. Покажем, что в точке P(G 4x 5| 0(X 5* 0))=0
Здесь G(x)=X - J 5-1 0(x) * F(x) - изображение итерационного процес са. Условие сходимости метода: P(G 4x 5| 0(X))
2. Выбор начального приближения X 50 0 - выбирается достаточно близко к точному решению. Как именно близко - зависит от скорости изменения функции F(x) вблизи X 5* 0: чем выше скорость, тем меньше область, где соблюдается условие сходимости и следовательно необходимо ближе выби рать X 50 0 к X 5* 0. 3. Скорость сходимости определяется соотношением ¦E 5m+1 0¦ = C 4m 0 * ¦E 5m 0¦ 51+p 0, 0
При X -> X 5* 0 величина P -> 1. Это значит, что вблизи точного реше ния скорость сходимости близка к квадратичной. Известно, что метода Ньютона сходится на 6-7 итерации.
4. Критерий окончания итераций - здесь могут использоваться кри терии точности, то есть максимальная из норм изменений X и F(x). ДИСКРЕТНЫЙ МЕТОД НЬЮТОНА Трудность использования метода Ньютона состоит в 1. Необходимости вычисления на каждом этапе J=F 4x 5| 0. Возможно несколько путей решения:
- аналитический способ. Наиболее эффективен, однако точные форму лы могут быть слишком большими, а также функции могут быть заданы таб лично. - конечно-разностная аппроксимация:
dF 4i 0 F 4i 0(x 41 0, ...., x 4j 0+dx 4j 0, ...., x 4n 0) F 4i 0(x 41 0, ...., x 4j 0-dx 4j 0, ....x 4n 0) --- = ------------------------------------------------- dX 4j 0 2*dX 4j
В этом случае мы имеем уже дискретный метод Ньютона, который уже не обладает квадратичной сходимостью. Скорость сходимости можно увели чить, уменьшая dX 4j 0 по мере приближения к X 5* 0.
2. Вычисление матрицы J 5-1 0 на каждом шаге требует значительных вы числительных затрат, поэтому часто решают такую систему: dX 5 0= 5 0J 5-1 0(X 5m 0) 5 0* 5 0F(X 5m 0) или умножая правую и левую часть на J, получим: J(X 5m 0) 5 0* 5 0dX 5m 0= 5 0F(X 5m 0)
На каждом M-ом шаге матрицы J и F известны. Необходимо найти dX 5m 0, как решение вышеприведенной системы линейных АУ, тогда X 5m+1 0=X 5m 0+dX 5m
Основной недостаток метода Ньютона - его локальная сходимость, поэтому рассмотрим еще один метод. МЕТОД ПРОДОЛЖЕНИЯ РЕШЕНИЯ ПО ПАРАМЕТРУ
Пусть t - параметр, меняющийся от 0 до 1. Введем в рассмотрение некоторую систему H(X, t)=0, такую, что: 1. при t=0 система имеет решение X 50 2. при t=1 система имеет решение X 5* 3. вектор-функция H(X, t) непрерывна по t.
Вектор функция может быть выбрана несколькими способами, например H(X, t) = F(X) + (t-1) или H(X, t) = t * F(X)
Нетрудно проверить, что для этих вектор-функций выполняются усло вия 1-3. Идея метода состоит в следующем:
Полагаем t 41 0= 7D 0t и решаем систему H(X, t 41 0)=0 при выбранном X 50 0. Полу
чаем вектор X 5t1 0. Далее берем его в качестве начального приближения и решаем при новом t 42 0=t 41 0+ 7D 0t систему H(X, t 42 0)=0, получаем X 5t2 0 и так далее до тех пор, пока не будет достигнута заданная точность. 3. ТЕХНОЛОГИЯ РАЗРЕЖЕННЫХ МАТРИЦ ОСНОВНЫЕ ИДЕИ МЕТОДА.
Основные требования к этим методам можно свести в 3 утверждения: 1. Хранить в памяти только ненулевые элементы.
2. В процессе решения принимать меры, уменьшающие возможность по явления новых ненулевых элементов. 3. Вычисления производить только с ненулевыми элементами. Разреженный строчный формат
Это одна из широко используемых схем для хранения разреженных матриц, которая предъявляет минимальные требования к памяти и очень удобна для выполнения основных операций с матрицами.
Значения ненулевых элементов и соответствующие столбцовые индексы хранятся по строкам в двух массивах: NA и JA. Для разметки строк в этих массивах используется массив указателей IA, отмечающих позиции массивов AN и JA, с которых начинается описание очередной строки. Пос ледняя цифра в массиве IA содержит указатель первой свободной позиции в JA и AN. Рассмотрим конкретный пример, который будет также и далее применятся для демонстрации других операций и который был использован в качестве контрольного примера для программы, выполняющей основные операции с разреженными матрицами. - ¬ ¦ 3 0 0 5 0 ¦ Позиция: 1 2 3 4 5 6 7 8 9 10 ¦ 0 4 0 0 1 ¦ IA: 1 3 5 7 9 11 A = ¦ 0 0 8 2 0 ¦ JA: 1 4 2 5 3 4 1 4 2 5 ¦ 5 0 0 6 0 ¦ AN: 3 5 4 1 8 2 5 6 7 9 ¦ 0 7 0 0 9 ¦ L
Данный способ представления называется полным (так как представ лена вся матрица А) и упорядоченным (так как элементы каждой строки хранятся в соответствии с возрастанием столбцовых индексов). Обознача ется RR(c)O - Row-wise representation Complete and Ordered (англ. ). Массивы IA и JA представляют портрет матрицы А. Если в алгоритме разграничены этапы символической и численной обработки, то массивы IA и JA заполняются на первом этапе, а массив AN - на втором. Транспонирование разреженной матрицы
Пусть IA, JA, AN - некоторое представление разреженной матрицы. Нужно получить IAT, JAT, ANT - разреженное представление транспониро ванной матрицы. Алгоритм рассмотрим на нашем примере: IA = 1 3 5 7 9 11 JA = 1 4 2 5 3 4 1 4 2 5 AN = 3 5 4 1 8 2 5 6 7 9 Символический этап.
1. Заводим 5 целых списков по числу столбцов матрицы А. пронуме руем их от 1 до 6.
2. Обходим 1 строку и заносим 1 в те списки, номера которых ука заны в JA: 1: 1 2: 3: 4: 1 5: 3. Обходим вторую строку, вставляя аналогичным образом 2: 1: 1 2: 2 3: 4: 1 5: 2 В итоге получим: 1: 1 4 2: 2 5 3: 3 4: 1 3 4 5: 2 5
Массив ANT можно получить, вставляя численное значение каждого ННЭ, хранимое в AN, в вещественный список сразу после того, как соот ветствующее целое внесено в целый список. В итоге получим: IAT = 1 3 5 6 9 11 JAT = 1 4 2 5 3 1 3 4 2 5 ANT = 3 5 4 7 8 5 2 6 1 9 Произведение разреженной матрицы и заполненного вектора-столбца Алгоритм рассмотрим на нашем конкретном примере: IAT = 1 3 5 7 9 11 JAT = 1 4 2 5 3 1 3 4 2 5 ANT = 3 5 4 7 8 5 2 6 1 9 B = ( -5 4 7 2 6 ) Обработка 1 строки:
Просматриваем массив IA и обнаруживаем, что первая строка матрицы А соответствует первому и второму элементу массива JA: JA(1)=3, JA(2)=4, то есть ННЭ являются a 411 0 и a 414 0.
Просматриваем массив AN и устанавливаем, что a 411 0=3 и a 414 0=5. Обращаемся к вектору B, выбирая из него соответственно b 41 0=-5 и b 44 0=2. Вычисляем C 41 0= 3 * (-5) + 5 * 2 = -5. Абсолютно аналогично, вычисляя остальные строки, получим: C = (-5 58 56 1 58). Произведение двух разреженных матриц
Рассмотрим метод на конкретном примере. Предположим, что нам не обходимо перемножить две матрицы: IA = 1 3 5 7 9 11 JA = 1 4 2 5 3 4 1 4 2 5 AN = 3 5 4 1 8 2 5 6 7 9 IAT = 1 3 5 7 9 11 JAT = 1 4 2 5 3 1 3 4 2 5 ANT = 3 5 4 7 8 5 2 6 1 9
Суть метода состоит в том, что используя метод переменного перек лючателя и расширенный вещественный накопитель, изменяется порядок на копления сумм для формирования элементов матрицы С. Мы будем формиро вать элементы C 4ij 0 не сразу, а постепенно накапливая, обращаясь только к строкам матрицы B. Символический этап. Определяем мерность IX = 5 IX = 0 0 0 0 0 1-я строка матрицы JAT: 1 4 JA(1) = 1 4 JC(1) = 1 4 IX = 1 0 0 1 0 JA(4) = 1 4 IX(1) = 1 ? Да. JC(1) - без изменений IX(4) = 1 ? Да. JC(1) - без изменений 2-я строка матрицы JAT: 2 5 JA(2) = 2 5 JC(2) = 2 5 IX = 1 2 0 1 2 JA(5) = 2 5 -> JC(2) - без изменений ...... 4-я строка матрицы JAT: 1 3 4 JA(1) = 1 4 JC(4) = 1 4 IX = 4 2 2 4 2 JA(3) = 3 4 IX(3) = 4 ? Нет. JC(4) = 1 4 3 IX(4) = 4 ? Да. JC(4) - без изменений ...... в итоге получим: IC = 1 3 5 7 10 12 JC = 1 4 2 5 3 4 1 4 3 2 5 Численный этап. X = 0 0 0 0 0 1-я строка: JC(1) = 1 4 AN(1) = 3 5, AA = 3 ANT(1) = 3 5 AA * ANT = 9 15 X = 9 0 0 15 0 AA = 5 ANT(1) = 3 5 AA * ANT = 15 25 X = 24 0 0 40 0 CN(1) = 24 40 ......
Аналогично проделывая действия над другими строками получим: IC: 1 3 5 7 10 12 JC: 1 4 2 5 3 4 1 4 3 2 5 CN: 24 40 20 35 80 0 55 22 66 16 144
Все вышеприведенные операции были получены на компьютере и ре зультаты совпали для нашего контрольного примера. Описание программы на языке 2 C 0, занимающейся этими операциями не приводится, так как дан ная программа нами не разрабатывалась, однако в приложении присутству ет распечатка этой программы. 2ПРАКТИЧЕСКАЯ ЧАСТЬ. ОПИСАНИЯ ПРОГРАММ. 1. ЯВНЫЕ МЕТОДЫ.
Данная группа представлена 3 программами, реализующими методы Эй лера, Рунге-Кутта 2-го и 4-го порядков. В данной группе все программы индентичны, поэтому далее следует описание программы, реализующем ме тод Эйлера, с указанием различий для методов Рунге-Кутта 2-го и 4-го порядков. 1EILER. M Головной модуль. Входные и выходные данные: отсутствуют. Язык реализации: PC MathLab Операционная система: MS-DOS 3. 30 or higher Пояснения к тексту модуля:
Стандартный головной модуль. Происходит очистка экрана, задание начальных значений по времени и по Y. Затем происходит вызов подпрог раммы EIL. M (RG2. M или RG4. M для методов Рунге-Кутта 2 и 4 порядков) а после получения результатов строятся графики. 1EIL. M Вычислительный модуль. Входные данные:
FunFcn - имя подпрограммы, написанной пользователем, которая возвращает левые части уравнения для определенного момента времени. T0, Tfinal - начальные и конечные моменты времени. Y0 - начальные значения. Выходные данные: Tout - столбец времени Yout - столбцы решений по каждой координате Язык реализации: PC MathLab Операционная система: MS-DOS 3. 30 or higher Пояснения к тексту модуля:
Данный модуль и реализует собственно метод Эйлера (Рунге-Кутта 2 или 4-го порядков). В начале происходит инициализация внутренних пере менных, а также вычисление максимального шага, который затем использу ется для определения точности. Далее следует цикл While.... End внутри которого и выполняются все необходимые действия: по формуле (для каж дого метода своя! ) вычисляется значение Y на каждом шаге (при необхо димости вызывается подфункция FunFcn) а также формируются выходные значения Tout и Yout. Далее метод сам корректирует свой шаг, по форму ле описанной выше (в теоретическом разделе). Этот цикл выполняется до тех пор, пока значение Tfinal не будет достигнуто. Все выходные значе ния формируются внутри данного цикла и поэтому никаких дополнительных действий не требуется. В связи с этим с окончанием цикла заканчивается и подпрограмма EIL. M. Методы Рунге-Кутта реализованы аналогично, с учетом отличий в формулах вычислений (более сложные по сравнению с ме тодом Эйлера). 2. НЕЯВНЫЕ МЕТОДЫ.
Представлены группой из двух похожих между собой программ, реали зующих соответственно неявные методы Эйлера и Рунге-Кутта 2 порядка. Также как и в вышеприведенном случае, будет описан метод Эйлера, а от личия метода Рунге-Кутта будут отмечены в скобках. 1NME. M Головной модуль. Входные и выходные данные отсутствуют. Язык реализации: PC MathLab Операционная система: MS-DOS 3. 30 or higher Пояснения к тексту модуля:
Выполняет стандартные действия: очистка экрана, инициализация и ввод начальных значений, вызов подпрограмм обработки и вычислений, а также построение графиков. 1NMEF. M, NRG2. M Вычислительные модули. Входные данные: T0, Tfinal - начальные и конечные моменты времени X0 - вектор-столбец начальных значений. H - начальный шаг
A - матрица, на диагонали которой стоят собственные числа линеа ризованной системы ОДУ. Выходные данные: T - столбец времени X - столбец решений 7D 0X - столбец ошибки Пояснения к тексту модуля:
Стандартные действия: инициализация начальных значений , цикл While T 3. МЕТОДЫ РЕШЕНИЯ НЕЛИНЕЙНЫХ САУ
Представлены группой из 4-х методов: метод последовательных приб лижений, метод Ньютона, метод Ньютона дискретный, метод продолжения решения по параметру. Метод последовательных приближений. 1MMPP. M Головной модуль. Входные и выходные данные отсутствуют. Язык реализации: PC MathLab Операционная система: MS-DOS 3. 30 or higher Пояснения к тексту модуля:
Очистка экрана, инициализация начальных значений, вызов вычисли тельного модуля MPP. M, вывод результатов в виде графиков. 1MPP. M Вычислительный модуль. Входные данные: X0 - начальное приближение в виде вектора-строки Fun1 - функция, возвращающая вычисленные левые части
Fun2 - функция, возвращающая матрицу Якоби в определенной точке. E - допустимая ошибка. Выходные данные: Mout - номера итераций Xout - приближения на каждой итерации DXout - ошибка на каждой итерации Язык реализации: PC MathLab Операционная система: MS-DOS 3. 30 or higher Пояснения к тексту модуля:
Аналогичен вышеприведенным вычислительным модулям - инициализация начальных значений, вычисления по формулам, вывод промежуточных ре зультатов, формирование выходных значений. По мере необходимости вызы вает подпрограммы Fun1 и Fun2. Методы Ньютона и Ньютона дискретный 1MNEWT. M Головной модуль. Входные и выходные данные отсутствуют. Язык реализации: PC MathLab Операционная система: MS-DOS 3. 30 or higher Пояснения к тексту модуля:
Стандартный головной модуль - выполняет действия, аналогичные предыдущим головным модулям. Вызывает из себя NEWT. M и NEWTD. M - моду ли реализующие методы Ньютона и Ньютона дискретный, а также строит их графики на одной координатной сетке. 1NEWT. M, NEWTD. M Вычислительные модули. Входные данные: X0 - начальное приближение в виде вектора-строки Fun1 - функция, возвращающая левые части
Fun2 - функция, вычисляющая матрицу Якоби (только для метода Ньютона! ) E - допустимая ошибка Выходные данные: Mout - номера итераций Xout - приближения на каждой итерации DXout - ошибка на каждой итерации Язык реализации: PC MathLab Операционная система: MS-DOS 3. 30 or higher Пояснения к тексту модулей:
Стандартные вычислительные модули, производящие соответствующие им действия. Отличие их в том, что в первом случае для вычисления мат рицы Якоби вызывается подпрограмма, а во втором случае матрица Якоби вычисляется внутри модуля. Метод продолжения решения по параметру 1MMPRPP. M Головной модуль. Входные и выходные данные отсутствуют. Язык реализации: PC MathLab Операционная система: MS-DOS 3. 30 or higher Пояснения к тексту модуля:
Стандартный головной модуль со всеми вытекающими отсюда последс твиями. 1MPRPP. M Вычислительный модуль. Входные данные: Fun1 - имя подпрограммы, вычисляющей правые части Fun2 - имя подпрограммы, вычисляющем матрицу Якоби X0 - начальное приближение dT - начальный шаг Edop - допустимая ошибка Trace - вывод на экран Язык реализации: PC MathLab Операционная система: MS-DOS 3. 30 or higher Пояснения к тексту модуля:
Стандартный вычислительный модуль - инициализация, вычисление, вывод, формирование результата. Стоит отметить, что поскольку метод имеет глобальную сходимость, то объем вычислений тут значительный, а это значит, что при выполнении вычислений требуется значительное коли чество свободной оперативной памяти. 2ВЫВОДЫ
При выполнении данных лабораторных работ были изучены основные численные методы для решения ОДУ, САУ, а также технология разреженных матриц. Заодно были получены основные навыки в программировании мате матической системы PC MathLab. Каждый из представленных методов по своему хорош и применяется для систем определенного вида. 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. Затем путем последовательного деления
t 5k 0 на 2 находим такое t 5k 0, чтобы выполнялось неравенство: ¦ F(X 5k+1 0) ¦ Итерационный процесс заканчивается при выполнении усло вия: ¦ F(X 5k+1 0) ¦ где 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 с полученными вход ными данными. Формат вызова данного модуля описан далее (в описании самого модуля). После вычислений в головную программу возвращаются ре
зультаты вычислений на основе которых строятся графики а также выводятся оценки по затратам машинного времени и быстродейс твия. 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+1 0)¦ Формирование выходных значений-матриц происходит внутри
данного цикла и поэтому никаких дополнительных действий не требуется, то есть с окончанием данного цикла заканчивается и сама функция. 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]; метод последовательной ите рации обратной матрицы Якоби уже не справляется, в то время как дискретный метод Ньютона продолжает неплохо работать. Од нако метод Ньютона требует больших затрат машинного времени и поэтому при выборе метода необходимо исходить их конкретных условий задачи и если известно довольно точное приближение и требуется быстрота вычислений, то к таким условиям отлично подходит разработанный метод последовательной итерации обрат ной матрицы Якоби. 2Выводы В данной РГР был разработан и реализован метод последова
тельной итерации обратной матрицы Якоби, предназначенный для решения системы нелинейных уравнений. Программа, реализованная на языке PC MathLab хотя и не является оптимальной, однако вы полняет поставленную задачу и решает системы уравнений. Реали зованный метод не отличается повышенной сходимостью и требует довольно точного начального приближения, однако довольно быст ро сходится к точному решению, то есть его можно порекомендо вать для вычисления непростых систем нелинейных уравнений при наличии довольно точного начального приближения и наличия вре менных ограничений. 2Список литературы 1. О. М. Сарычева. "Численные методы в экономике. Конспект
лекций", Новосибирский государственный технический универси тет, Новосибирск 1995г. 2. Д. Мак-Кракен, У. Дорн. "Численные методы и программиро вание на Фортране", Издательство "Мир", М. 1977г. 3. Н. С. Бахвалов. "Численные методы", Издательство "Нау ка", М. 1975г.