Главная страница

Численные методы. У которой известны входы и выходы, а то, что


Скачать 251.52 Kb.
НазваниеУ которой известны входы и выходы, а то, что
Дата28.10.2021
Размер251.52 Kb.
Формат файлаrtf
Имя файлаЧисленные методы.rtf
ТипЗакон
#257782



 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г.


написать администратору сайта