"Комплект" заданий по численным методам
Формат:
Дата создания: 10.10.2016
Размер: 10.9 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 │ └ ┘ Данный способ представления называется полным (так как представ- лена вся матрица А) и упорядоченным (так как элементы каждой строки хранятся в соответствии с возрастанием столбцовых индексов). Обознача- ется 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. Каждый из представленных методов по своему хорош и применяется для систем определенного вида.