Численные методы. У которой известны входы и выходы, а то, что
 Скачать 251.52 Kb.
|
|
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< 7 2 0e 5hl 7 2 0 (6) 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 < 2* 7t 0 . Для уравнений n-ого порядка : h 4i 0 <= 2* 7t 4i 0. А окончательно шаг выбирают по условиям устойчивости : h < 2* 7t 4min 0 , 7t 4min 0 = min 7t 4i Вначале задаётся допустимая ошибка аппроксимации , а в процессе ин- тегрирования шаг подбирается следующим образом : 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 <= E 4i 5aдоп 7/ 0[ 72 0f 4i 0(x 5m 0,t 4m 0) 72 0 + E 4i 5aдоп 7/ 0h 4max 0] (7) 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 <= 2* 7t 4min , 0 (2) где 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 < 0. Областью абсолютной устойчивости метода интегрирования системы (5) будет вся левая полуплоскость. Таким образом , шаг h должен на каждом интервале интегрирования подбираться таким образом, чтобы при этом не покидать эту область. Но в таком случае должно выполняться требование : h < 2* 7t 0, (6) где 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/ 7a 4max 0 или h < 2* 7t 4min 0. _Случай 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 <= t <= t 4m+1 0, i=[1;n]. 3. _Выбор шага интегрирования .. Должны соблюдаться условия абсолютной (6) или относительной (7) устойчивости в зависимости от характера интегрируемой системы. Для уравнения первого порядка шаг должен быть : h < 2* 7t 0 . Для уравнений n-ого порядка : h 4i 0 <= 2* 7t 4i 0. Однако область абсолютной устойчивости - вся левая полуплос- кость . Поэтому шаг с этой точки зрения может быть любым. Условия относительной устойчивости жестче, так как есть об- ласть , где они могут быть нарушены. Эти условия будут выполняться , если в процессе решения задачи контролировать ошибку аппроксимации , а шаг корректировать . Таким образом, шаг можно выбирать из условий точности, при этом условия устойчивости будут соблюдены автоматичес- ки. Сначала задается допустимая погрешность аппроксимации : E 4i 5aдоп 0 <= 0,001 72 0x 4i 72 4max 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 < E 4i 5aдоп 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)) < 1 должно выполняться для всех значений X 5m 0. Это условие осуществляется при достаточно жестких требованиях к F(x): функция должна быть непрерывна и дифференцируема по X, а также выпукла или вогнута вблизи X 5* 0. При этом выполняется лишь условие локальной сходимости. Причем можно показать, что чем ближе к X 5* 0, тем быстрее сходится метод. 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 < P < 1. При 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 < Tfinal, вычисление по формулам, вывод промежуточных резуль- татов, формирование выходных значений массивов. 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 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 с полученными вход- ными данными. Формат вызова данного модуля описан далее (в описании самого модуля). После вычислений в головную программу возвращаются ре- зультаты вычислений на основе которых строятся графики а также выводятся оценки по затратам машинного времени и быстродейс- твия. 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. Итерационный процесс заканчивается, если ¦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]; метод последовательной ите- рации обратной матрицы Якоби уже не справляется, в то время как дискретный метод Ньютона продолжает неплохо работать. Од- нако метод Ньютона требует больших затрат машинного времени и поэтому при выборе метода необходимо исходить их конкретных условий задачи и если известно довольно точное приближение и требуется быстрота вычислений, то к таким условиям отлично подходит разработанный метод последовательной итерации обрат- ной матрицы Якоби. 2Выводы В данной РГР был разработан и реализован метод последова- тельной итерации обратной матрицы Якоби, предназначенный для решения системы нелинейных уравнений. Программа, реализованная на языке PC MathLab хотя и не является оптимальной, однако вы- полняет поставленную задачу и решает системы уравнений. Реали- зованный метод не отличается повышенной сходимостью и требует довольно точного начального приближения, однако довольно быст- ро сходится к точному решению, то есть его можно порекомендо- вать для вычисления непростых систем нелинейных уравнений при наличии довольно точного начального приближения и наличия вре- менных ограничений. 2Список литературы 1. О.М.Сарычева. "Численные методы в экономике. Конспект лекций", Новосибирский государственный технический универси- тет, Новосибирск 1995г. 2. Д.Мак-Кракен, У.Дорн. "Численные методы и программиро- вание на Фортране", Издательство "Мир", М. 1977г. 3. Н.С.Бахвалов. "Численные методы", Издательство "Нау- ка", М. 1975г. |