"Комплект" заданий по численным методам
ВВЕДЕНИЕ
В экономике очень часто используется модель, называемая "черный
ящик", то есть система у которой известны входы и выходы, то, что
происходит внутри - неизвестно. Законы по которым происходят изменения
выходных сигналов в зависимости от входных могут быть различными, в
том числе это могут быть и дифференциальные законы. Тогда встает проб-
лема решения систем дифференциальных равнений, которые в зависимости
от своей структуры могут быть решены различными методами. Точное реше-
ние получить очень часто не дается, поэтому мы рассмотрим численные
методы решения таких систем. Далее мы представим две группы методов:
явные и неявные. Для решения систем ОДУ неявными методами придется ре-
шать системы нелинейных равнений, поэтому придется ввести в рассмот-
рение группу методов решения систем нелинейных равнений, которые в
свою очередь будут представлены двумя методами. Далее следуют теорети-
ческие аспекты описанных методов, затем будут представлены описания
программ. Сами программы, также их графики приведены в приложении.
Также стоит отметить, что в принципе все численные методы так или
иначе сводятся к матричной алгебре, в экономических задачах очень
часто матрицы имеют слабую заполненность и большие размеры и поэтому
неэффективно работать с полными матрицами. Одна из технологий, позво-
ляющая разрешить данную проблему - технология разреженных матриц.
В связи с этим, мы рассмотрим данную технологию и операции множения и
транспонирования над такими матрицами.
Таким образом мы рассмотрим весь спектр лабораторных работ. Опи-
сания всех программ приводятся после теоретической части. Все тексты
программ и распечатки графиков приведены в приложении.
ТЕОРЕТИЧЕСКАЯ ЧАСТЬ
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 Е 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, занимающейся этими операциями не приводится, так как дан-
ная программа нами не разрабатывалась, однако в приложении присутству-
ет распечатка этой программы.
ПРАКТИЧЕСКАЯ ЧАСТЬ. ОПИСАНИЯ ПРОГРАММ.
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
Пояснения к тексту модуля:
Стандартный вычислительный модуль - инициализация, вычисление,
вывод, формирование результата. Стоит отметить, что поскольку метод
имеет глобальную сходимость, то объем вычислений тут значительный, а
это значит, что при выполнении вычислений требуется значительное коли-
чество свободной оперативной памяти.
ВЫВОДЫ
При выполнении данных лабораторных работ были изучены основные
численные методы для решения ОДУ, САУ, также технология разреженных
матриц. Заодно были получены основные навыки в программировании мате-
матической системы PC MathLab. Каждый из представленных методов по
своему хорош и применяется для систем определенного вида.