Реферат: "Комплект" заданий по численным методам
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 01+h* 7l 0 72 0< 7 2 0e 5hl 72 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 01/(1-h* 7l 0) 72 0 > 1.
С учетом ограничения на скоростьизменения приближенного ре-
шения относительно точного :
72 01/(1-h* 7l 0) 72 0 > 7 2 0e 5hl 72 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 78 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 14 2 5
│ 5 0 0 6 0 │ AN: 3 5 4 1 8 2 56 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
Все вышеприведенные операции были полученына компьютере и ре-
зультатысовпали для нашего контрольногопримера. Описание программы
на языке 2C 0, занимающейся этими операциямине приводится, так как дан-
ная программанами не разрабатывалась, однако в приложении присутству-
ет распечаткаэтой программы.
2ПРАКТИЧЕСКАЯ ЧАСТЬ.ОПИСАНИЯ ПРОГРАММ.
1. ЯВНЫЕ МЕТОДЫ.
Данная группа представлена 3 программами,реализующими методы Эй-
лера, Рунге-Кутта2-го и 4-го порядков. В данной группевсе программы
индентичны,поэтому далее следует описание программы, реализующем ме-
тод Эйлера, с указанием различий для методов Рунге-Кутта2-го и 4-го
порядков.
1EILER.M
Головной модуль.
Входные и выходные данные: отсутствуют.
Язык реализации: PC MathLab
Операционная система: MS-DOS 3.30 orhigher
Пояснения к тексту модуля:
Стандартный головной модуль. Происходит очистка экрана, задание
начальныхзначений по времени и по Y. Затемпроисходит вызов подпрог-
раммы EIL.M(RG2.M или RG4.M для методов Рунге-Кутта 2 и 4 порядков) а
после получениярезультатов строятся графики.
1EIL.M
Вычислительный модуль.
Входные данные:
FunFcn - имя подпрограммы, написанной пользователем, которая
возвращает левыечасти уравнения для определенного момента времени.
T0, Tfinal — начальные и конечные моментывремени.
Y0 — начальные значения.
Выходные данные:
Tout — столбец времени
Yout — столбцы решений по каждойкоординате
Язык реализации: PC MathLab
Операционная система: MS-DOS 3.30 orhigher
Пояснения к тексту модуля:
Данный модуль и реализует собственно метод Эйлера(Рунге-Кутта 2
или 4-гопорядков). В начале происходит инициализация внутренних пере-
менных, а такжевычисление максимального шага, который затем использу-
ется дляопределения точности. Далее следует циклWhile...End внутри
которого ивыполняются все необходимые действия: поформуле (для каж-
дого методасвоя!) вычисляется значение Y на каждом шаге (при необхо-
димостивызывается подфункция FunFcn) а также формируются выходные
значения Tout иYout. Далее метод сам корректирует свой шаг, по форму-
ле описанной выше(в теоретическом разделе). Этот циклвыполняется до
тех пор, показначение Tfinal не будет достигнуто. Все выходные значе-
нияформируются внутри данного цикла ипоэтому никаких дополнительных
действий нетребуется. В связи с этим с окончанием цикла заканчивается
и подпрограмма EIL.M. Методы Рунге-Кутта реализованыаналогично, с
учетом отличий вформулах вычислений (более сложные по сравнению с ме-
тодом Эйлера).
2. НЕЯВНЫЕ МЕТОДЫ.
Представлены группой из двух похожих междусобой программ, реали-
зующих соответственнонеявные методы Эйлера и Рунге-Кутта 2 порядка.
Также как и ввышеприведенном случае, будет описан метод Эйлера, а от-
личия методаРунге-Кутта будут отмечены в скобках.
1NME.M
Головной модуль.
Входные и выходные данные отсутствуют.
Язык реализации: PC MathLab
Операционная система: MS-DOS 3.30 orhigher
Пояснения к тексту модуля:
Выполняет стандартные действия: очистка экрана, инициализация и
ввод начальныхзначений, вызов подпрограмм обработки ивычислений, а
также построениеграфиков.
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 orhigher
Пояснения к тексту модуля:
Очистка экрана, инициализация начальных значений, вызоввычисли-
тельного модуляMPP.M, вывод результатов в виде графиков.
1MPP.M
Вычислительный модуль.
Входные данные:
X0 — начальное приближение в видевектора-строки
Fun1 — функция, возвращающая вычисленныелевые части
Fun2 — функция, возвращающая матрицу Якобив определенной точке.
E — допустимая ошибка.
Выходные данные:
Mout — номера итераций
Xout — приближения на каждой итерации
DXout — ошибка на каждой итерации
Язык реализации: PC MathLab
Операционная система: MS-DOS 3.30 orhigher
Пояснения к тексту модуля:
Аналогичен вышеприведенным вычислительныммодулям — инициализация
начальныхзначений, вычисления по формулам, вывод промежуточных ре-
зультатов,формирование выходных значений. По мере необходимости вызы-
вает подпрограммыFun1 и Fun2.
Методы Ньютона и Ньютонадискретный
1MNEWT.M
Головной модуль.
Входные и выходные данные отсутствуют.
Язык реализации: PC MathLab
Операционная система: MS-DOS 3.30 orhigher
Пояснения к тексту модуля:
Стандартный головной модуль - выполняет действия, аналогичные
предыдущимголовным модулям. Вызывает из себя NEWT.M и NEWTD.M — моду-
ли реализующиеметоды Ньютона и Ньютона дискретный, атакже строит их
графики на однойкоординатной сетке.
1NEWT.M, NEWTD.M
Вычислительные модули.
Входные данные:
X0 — начальное приближение в видевектора-строки
Fun1 — функция, возвращающая левые части
Fun2 — функция, вычисляющая матрицу Якоби (только для метода
Ньютона!)
E — допустимая ошибка
Выходные данные:
Mout — номера итераций
Xout — приближения на каждой итерации
DXout — ошибка на каждой итерации
Язык реализации: PC MathLab
Операционная система: MS-DOS 3.30 orhigher
Пояснения к тексту модулей:
Стандартные вычислительные модули, производящие соответствующие
им действия.Отличие их в том, что в первом случае для вычисления мат-
рицы Якобивызывается подпрограмма, а во второмслучае матрица Якоби
вычисляетсявнутри модуля.
Метод продолжения решения попараметру
1MMPRPP.M
Головной модуль.
Входные и выходные данные отсутствуют.
Язык реализации: PC MathLab
Операционная система: MS-DOS 3.30 orhigher
Пояснения к тексту модуля:
Стандартный головной модуль со всемивытекающими отсюда последс-
твиями.
1MPRPP.M
Вычислительный модуль.
Входные данные:
Fun1 — имя подпрограммы, вычисляющейправые части
Fun2 — имя подпрограммы, вычисляющемматрицу Якоби
X0 — начальное приближение
dT — начальны