"Комплект" заданий по численным методам


Сдавался/использовался	1996г.
Примечание	4 лабораторных + 2 РГР
	Загрузить архив:
	Файл: chisl.zip (32kb [zip], Скачиваний: 351) скачать

2ВВЕДЕНИЕ

В экономике очень часто используется модель, называемая "черный

ящик", то есть система у которой известны входы и выходы, а то, что

происходит внутри - неизвестно. Законы по которым происходят изменения

выходных сигналов в зависимости от входных могут быть различными, в

том числе это могут быть и дифференциальные законы. Тогда встает проб-

лема решения систем дифференциальных уравнений, которые в зависимости

от своей структуры могут быть решены различными методами. Точное реше-

ние получить очень часто не удается, поэтому мы рассмотрим численные

методы решения таких систем. Далее мы представим две группы методов:

явные и неявные. Для решения систем ОДУ неявными методами придется ре-

шать системы нелинейных уравнений, поэтому придется ввести в рассмот-

рение группу методов решения систем нелинейных уравнений, которые в

свою очередь будут представлены двумя методами. Далее следуют теорети-

ческие аспекты описанных методов, а затем будут представлены описания

программ. Сами программы, а также их графики приведены в приложении.

Также стоит отметить, что в принципе все численные методы так или

иначе сводятся к матричной алгебре, а в экономических задачах очень

часто матрицы имеют слабую заполненность и большие размеры и поэтому

неэффективно работать с полными матрицами. Одна из технологий, позво-

ляющая разрешить данную проблему - технология разреженных матриц.

В связи с этим, мы рассмотрим данную технологию и операции умножения и

транспонирования над такими матрицами.

Таким образом мы рассмотрим весь спектр лабораторных работ. Опи-

сания всех программ приводятся после теоретической части. Все тексты

программ и распечатки графиков приведены в приложении.

2ТЕОРЕТИЧЕСКАЯ ЧАСТЬ

1. ОПИСАНИЕ МЕТОДОВ ИНТЕГРИРОВАНИЯ СИСТЕМ ОДУ

И ИХ ХАРАКТЕРИСТИК

ЯВНЫЙ МЕТОД ЭЙЛЕРА И ЕГО ХАРАКТЕРИСТИКИ

Алгоритм этого метода определяется формулой:

x5m+10 = x5m0 + h4m0*F(x5m0, t4m0)4,0 (1)

которая получается путём аппроксимации ряда Тейлора до членов перво-

го порядка производной x'(t4m0),т.к. порядок точности равен 1 (s=1).

Для аналитического исследования свойств метода Эйлера линеари-

зуется исходная система ОДУ x' = F(x, t) в точке (x5m0,t4m0):

x' = A*x, (2)

где матрица А зависит от точки линеаризации (x5m0,t4m0).

Входной сигнал при линеаризации является неизвестной функцией

времени и при фиксированном t4m0 на шаге h4m0 может считаться констан-

той. Ввиду того ,что для линейной системы свойство устойчивости за-

висит лишь от А, то входной сигнал в системе (2) не показан. Элемен-

ты матрицы А меняются с изменением точки линеаризации,т.е. с измене-

нием m.

Характеристики метода:

1. _Численная устойчивость..

Приведя матрицу А к диагональной форме : А = Р*7l0*Р5-10 и перейдя

к новым переменным y = P5-10*x , система (2) примет вид :

y' = 7l0*y (3)

Тогда метод Эйлера для уравнения (3) будет иметь вид :

y5m+10 = y5m0 + h*7l0 * y5m0, (4)

где величина h*7l0 изменяется от шага к шагу.

В этом случае характеристическое уравнение для дискретной сис-

темы (4) будет иметь вид :

1 + h*7l0 - r = 0.

А корень характеристического уравнения равен:

r = 1+ h*7l0,

где 7l0 =7 a0 _+.7 b0 .

_Случай 1.. Исходная система (3) является асимптотически устой-

чивой , т.е. нулевое состояние равновесия системы асимптотически ус-

тойчиво и 7 a0 < 0.

Областью абсолютной устойчивости метода интегрирования системы

(4) будет круг радиусом , равным 1 , и с центром в точке (0, -1).

Таким образом , шаг h должен на каждом интервале интегрирования под-

бираться таким образом, чтобы при этом не покидать область А . Но

в таком случае должно выполняться требование :

h < 2*7t0, (5)

где 7t0=1/72a20 - постоянная времени системы (3) . Она определяет ско-

рость затухания переходных процессов в ней. А время переходного про-

цесса T4пп0 = 3*7t0 , где 7t0 = 72a25-10.

Если иметь ввиду, что система (3) n-го порядка, то

T4пп0 > 3*7t4max0,

где 7t4max0 = 72a4min725-17 0; 7a4min 0= min 7a4i0 . Из всех 7a4i0 (i=[1;n]) выбирает-

ся максимальное значение : max72a4i720=7a4max0 и тогда можно вычислить :

7t4min 0= 1/7a4max0,

а шаг должен выбираться следующим образом :

h < 2/7a4max0 или h < 2*7t4min0.

_Случай 2.. Нулевое состояние равновесия системы (2) неустойчи-

во, т.е. 7a0 > 0 . Т.е. система тоже неустойчива , т.е. 720r720>1. Поэтому

областью относительной устойчивости явного метода Эйлера является

вся правая полуплоскость . Выполняется требование :

720 1+h*7l0 72 0<7 2 0e5hl7 20 (6)

2. _Точность метода..

Так как формула интегрирования (1) аппроксимирует ряд Тейлора

для функции x(t4m0+1) до линейного по h члена включительно. Существует

такое значение t в интервале [t4m0 , t4m+10], при котором

Е4i5am0 =1/2! * h4m520*x4i0''(t),

где i=[1;n].

3. _Выбор шага интегрирования..

Должны соблюдаться условия абсолютной (5) или относительной

(6) устойчивости в зависимости от характера интегрируемой системы.

Для уравнения первого порядка шаг должен быть :

h < 2*7t0 .

Для уравнений n-ого порядка :

h4i0 <= 2*7t4i 0.

А окончательно шаг выбирают по условиям устойчивости :

h < 2*7t4min0 , 7t4min0 = min 7t4i

Вначале задаётся допустимая ошибка аппроксимации , а в процессе ин-

тегрирования шаг подбирается следующим образом :

1) по формуле (1) определяется очередное значение x5m+10;

2) определяется dx4i5m0 = x4i5m+10 - x4i5m0 ;

3) условие соблюдения точности имеет вид :

h4i5m0 <= E4i5aдоп7/0[720f4i0(x5m0,t4m0)720 + E4i5aдоп7/0h4max0] (7)

4) окончательно на m-м интервале времени выбирается в виде:

h4m0 = min h4i5m0.

ЯВНЫЕ МЕТОДЫ РУНГЕ-КУТТА

Метод Эйлера является методом Рунге-Кутта 1-го порядка .

Методы Рунге-Кутта 2-го и 4-го порядка являются одношаговыми ,

согласуются с рядом Тейлора до порядка точности s , который равен

порядку метода . Эти методы не требуют вычисления производных

функций , а только самой функции в нескольких точках на шаге h4m0.

Алгоритм метода Рунге-Кутта 2-го порядка состоит в следующем :

x5m+10 = x5m0 + h4m0/2 (k410+k420),

где k410=f(x5m0,t4m0) ; k420=f(x5m0+h4m0*k410,t4m0+h4m0).

Ошибка аппроксимации Е5a0 = k*h4m530 .

Алгоритм метода Рунге-Кутта 4-го порядка

x5m+10=x5m0+h4m0/6(k410+2k420+2k430+k440),

где k410=f(x5m0,t4m0); k420=f(x5m0+h4m0/2*k410,t4m0+h4m0/2); k430=f(x5m0+h4m0/2*k420,t4m0+h4m0/2);

k440=f(x5m0+h4m0*k430,t4m0+h4m0).

Ошибка аппроксимации Е5a0 = k*h4m550.

НЕЯВНЫЙ МЕТОД ЭЙЛЕРА И ЕГО ХАРАКТЕРИСТИКИ

Неявный метод Эйлера используется для интегрирования " жест-

ких " систем. "Жесткие" системы это такие системы, в которых7 0 (7l4max0)

и (7l4min0) сильно отключаются друг от друга , то в решениях системы

x' = A*x (1)

будут присутствовать экспоненты, сильно отличаются друг от друга по

скорости затухания . Шаг интегрирования для таких систем должен вы-

бираться по условиям устойчивости из неравенства

h <= 2*7t4min ,0 (2)

где 7t0=1/72a20 - постоянная времени системы y' = 7l0*y . Она определяет

скорость затухания переходных процессов в ней . Неравенство (2)

должно выполняться на всем участке решения , что соответственно тре-

бует значительных затрат машинного времени.

Алгоритм этого метода определяется формулой:

x5m+10 = x5m0 + h4m0*F(x5m+10, t4m+10) 4 0(3)

Если h4m0 мал, то x5m0 и х5m+10 близки друг к другу. В качестве на-

чального приближения берется точка x5m0 , а следовательно , между x5m0 и

x5m+10 будет существовать итерационный процесс.

Для аналитического исследования свойств метода Эйлера линеа-

ризуется исходная система ОДУ x' = F(x, t) в точке (x5m0,t4m0):

x' = A*x,

где матрица А зависит от точки линеаризации (x5m0,t4m0).

Входной сигнал при линеаризации является неизвестной функцией

времени и при фиксированном t4m0 на шаге h4m0 может считаться констан-

той. Ввиду того ,что для линейной системы свойство устойчивости за-

висит лишь от А, то входной сигнал в системе (1) не показан. Элемен-

ты матрицы А меняются с изменением точки линеаризации,т.е. с измене-

нием m.

Характеристики метода:

1. _Численная устойчивость..

Приведя матрицу А к диагональной форме : А = Р*7l0*Р5-10 и перейдя

к новым переменным y = P5-10*x , система (3) примет вид :

y' = 7l0*y (4)

Тогда метод Эйлера для уравнения (4) будет иметь вид :

y5m+10 = y5m0 + h*7l0 * y5m+10, (5)

где величина h*7l0 изменяется от шага к шагу.

В этом случае характеристическое уравнение для дискретной сис-

темы (5) будет иметь вид :

1 - h*7l0*r - 1 = 0.

А корень характеристического уравнения равен:

r = 1/(1-h*7l0) ,

где 7l0 =7 a0 _+.7 b0 .

_Случай 1.. Исходная система (4) является асимптотически устой-

чивой , т.е. нулевое состояние равновесия системы асимптотически ус-

тойчиво и 7 a0 < 0.

Областью абсолютной устойчивости метода интегрирования системы

(5) будет вся левая полуплоскость. Таким образом , шаг h должен на

каждом интервале интегрирования подбираться таким образом, чтобы при

этом не покидать эту область. Но в таком случае должно выполняться

требование :

h < 2*7t0, (6)

где 7t0=1/72a20 - постоянная времени системы (4) . Она определяет ско-

рость затухания переходных процессов в ней. А время переходного про-

цесса T4пп0 = 3*7t0 , где 7t0 = 72a25-10.

Если иметь ввиду, что система (4) n-го порядка, то

T4пп0 > 3*7t4max0,

где 7t4max0 = 72a4min725-17 0; 7a4min 0= min 7a4i0 . Из всех 7a4i0 (i=[1;n]) выбирает-

ся максимальное значение : max72a4i720=7a4max0 и тогда можно вычислить :

7t4min 0= 1/7a4max0,

а шаг должен выбираться следующим образом :

h < 2/7a4max0 или h < 2*7t4min0.

_Случай 2.. Нулевое состояние равновесия системы (4) неустойчи-

во, т.е. 7a0 > 0 . Т.е. система тоже неустойчива , т.е. 720r720>1,

а следовательно :

720 1/(1-h*7l0) 720 > 1.

С учетом ограничения на скорость изменения приближенного ре-

шения относительно точного :

720 1/(1-h*7l0) 720 >7 2 0e5hl7 20. (7)

Из этого соотношения следует , что при 720h*7l20 -> 1 левая часть

стремится к бесконечности . Это означает , что в правой полуплоскос-

ти есть некоторый круг радиусом , равным 1 , и с центром в точке

(0, 1), где неравенство (7) не выполняется.

2. _Точность метода..

Ошибка аппроксимации по величине равна ошибке аппроксимации

явного метода Эйлера , но она противоположна по знаку :

Е4i5am0 =-1/2! * h4m520*x4i0''(t),

где t4m0 <= t <= t4m+10,

i=[1;n].

3. _Выбор шага интегрирования..

Должны соблюдаться условия абсолютной (6) или относительной

(7) устойчивости в зависимости от характера интегрируемой системы.

Для уравнения первого порядка шаг должен быть :

h < 2*7t0 .

Для уравнений n-ого порядка :

h4i0 <= 2*7t4i 0.

Однако область абсолютной устойчивости - вся левая полуплос-

кость . Поэтому шаг с этой точки зрения может быть любым.

Условия относительной устойчивости жестче, так как есть об-

ласть , где они могут быть нарушены. Эти условия будут выполняться ,

если в процессе решения задачи контролировать ошибку аппроксимации ,

а шаг корректировать . Таким образом, шаг можно выбирать из условий

точности, при этом условия устойчивости будут соблюдены автоматичес-

ки. Сначала задается допустимая погрешность аппроксимации :

E4i5aдоп0 <= 0,001 720x4i724max0,

где i=[1;n].

Шаг выбирается в процессе интегрирования следующим образом:

1. Решая систему (3) относительно x5m+10 с шагом h4m0, получаем

очередную точку решения системы x = F(x,t) ;

2. Оцениваем величину ошибки аппроксимации по формуле:

Е4i5am0 = 720h4m7/0(h4m0+h4m-10)*[(x4i5m+10 - x4i5m0) - h4m7/0h4m-10*(x4i5m0 -x4i5m-10)]72

3. Действительное значение аппроксимации сравнивается с до-

пустимым. Если Е4i5am0 < E4i5aдоп0, то выполняется следующий пункт, в про-

тивном случае шаг уменьшается в два раза , и вычисления повторяются

с п.1.

4. Рассчитываем следующий шаг:

h4i5m+10 = SQR(E4i5aдоп7/20Е4i5am720) * h4m0.

5. Шаг выбирается одинаковым для всех элементов вектора X :

h4m+10 = min h4i5m+10.

6. Вычисляется новый момент времени и алгоритм повторяется .

НЕЯВНЫЕ МЕТОДЫ РУНГЕ-КУТТА

Метод Эйлера является методом Рунге-Кутта 1-го порядка .

Методы Рунге-Кутта 2-го и 4-го порядка являются одношаговыми ,

согласуются с рядом Тейлора до порядка точности s , который равен

порядку метода . Эти методы не требуют вычисления производных

функций , а только самой функции в нескольких точках на шаге h4m0.

Алгоритм метода Рунге-Кутта 2-го порядка состоит в следующем:

x5m+10 = x5m0 + h4m0/2 (k415m+10+k425m+10),

где k415m+10=f(x5m+10,t4m+10); k425m+10=f(x5m+10-h4m0*k415m+10,t4m+10).

Ошибка аппроксимации Е4m5a0 = k*h4m530 .

Алгоритм метода Рунге-Кутта 4-го порядка

x5m+10 = x5m0 + h4m0/6 (k415m+10+2k425m+10+2k435m+10+k445m+10),

где k410=f(x5m+10,t4m+10); k420=f(x5m+10-h4m0/2*k415m+10,t4m+10-h4m0/2);

k430=f(x5m+10-h4m0/2*k425m+10,t4m+10-h4m0/2); k440=f(x5m+10-h4m0*k435m+10,t4m0).

Ошибка аппроксимации Е5a0 = k*h4m550.

2. МЕТОДЫ РЕШЕНИЯ НЕЛИНЕЙНЫХ САУ

МЕТОД НЬЮТОНА

Итерационная формула метода Ньютона имеет вид

X5m+10=X5m 0-5 0J5-1 0*5 0(X5m0)5 0*5 0F(X5m0),

где J(X)=F4x5|0/4x=xm

Характеристики метода:

1. Сходимость. Покажем, что в точке P(G4x5|0(X5*0))=0

Здесь G(x)=X - J5-10(x) * F(x) - изображение итерационного процес-

са. Условие сходимости метода: P(G4x5|0(X)) < 1 должно выполняться для

всех значений X5m0. Это условие осуществляется при достаточно жестких

требованиях к F(x): функция должна быть непрерывна и дифференцируема

по X, а также выпукла или вогнута вблизи X5*0. При этом выполняется лишь

условие локальной сходимости. Причем можно показать, что чем ближе к

X5*0, тем быстрее сходится метод.

2. Выбор начального приближения X500 - выбирается достаточно близко

к точному решению. Как именно близко - зависит от скорости изменения

функции F(x) вблизи X5*0: чем выше скорость, тем меньше область, где

соблюдается условие сходимости и следовательно необходимо ближе выби-

рать X500 к X5*0.

3. Скорость сходимости определяется соотношением

║E5m+10║ = C4m0 * ║E5m0║51+p0, 0 < P < 1.

При X -> X5*0 величина P -> 1. Это значит, что вблизи точного реше-

ния скорость сходимости близка к квадратичной. Известно, что метода

Ньютона сходится на 6-7 итерации.

4. Критерий окончания итераций - здесь могут использоваться кри-

терии точности, то есть максимальная из норм изменений X и F(x).

ДИСКРЕТНЫЙ МЕТОД НЬЮТОНА

Трудность использования метода Ньютона состоит в

1. Необходимости вычисления на каждом этапе J=F4x5|0.

Возможно несколько путей решения:

- аналитический способ. Наиболее эффективен, однако точные форму-

лы могут быть слишком большими, а также функции могут быть заданы таб-

лично.

- конечно-разностная аппроксимация:

dF4i0 F4i0(x410,...,x4j0+dx4j0,...,x4n0) - F4i0(x410,...,x4j0-dx4j0,...x4n0)

─── = ──────────────────────────────────────────────────

dX4j0 2*dX4j

В этом случае мы имеем уже дискретный метод Ньютона, который уже

не обладает квадратичной сходимостью. Скорость сходимости можно увели-

чить, уменьшая dX4j0 по мере приближения к X5*0.

2. Вычисление матрицы J5-10 на каждом шаге требует значительных вы-

числительных затрат, поэтому часто решают такую систему:

dX5 0=5 0J5-10(X5m0)5 0*5 0F(X5m0)

или умножая правую и левую часть на J, получим:

J(X5m0)5 0*5 0dX5m 0=5 0F(X5m0)

На каждом M-ом шаге матрицы J и F известны. Необходимо найти dX5m0,

как решение вышеприведенной системы линейных АУ, тогда

X5m+10=X5m0+dX5m

Основной недостаток метода Ньютона - его локальная сходимость,

поэтому рассмотрим еще один метод.

МЕТОД ПРОДОЛЖЕНИЯ РЕШЕНИЯ ПО ПАРАМЕТРУ

Пусть t - параметр, меняющийся от 0 до 1. Введем в рассмотрение

некоторую систему

H(X,t)=0,

такую, что:

1. при t=0 система имеет решение X50

2. при t=1 система имеет решение X5*

3. вектор-функция H(X,t) непрерывна по t.

Вектор функция может быть выбрана несколькими способами, например

H(X,t) = F(X) + (t-1)

или

H(X,t) = t * F(X)

Нетрудно проверить, что для этих вектор-функций выполняются усло-

вия 1-3.

Идея метода состоит в следующем:

Полагаем t410=7D0t и решаем систему H(X,t410)=0 при выбранном X500. Полу-

чаем вектор X5t10. Далее берем его в качестве начального приближения и

решаем при новом t420=t410+7D0t систему H(X,t420)=0, получаем X5t20 и так далее

до тех пор, пока не будет достигнута заданная точность.

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

4: 1

3. Обходим вторую строку, вставляя аналогичным образом 2:

1: 1

2: 2

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, то есть ННЭ являются a4110 и a4140.

Просматриваем массив AN и устанавливаем, что a4110=3 и a4140=5.

Обращаемся к вектору B, выбирая из него соответственно b410=-5 и

b440=2.

Вычисляем C410= 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

Суть метода состоит в том, что используя метод переменного перек-

лючателя и расширенный вещественный накопитель, изменяется порядок на-

копления сумм для формирования элементов матрицы С. Мы будем формиро-

вать элементы C4ij0 не сразу, а постепенно накапливая, обращаясь только

к строкам матрицы 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 C0, занимающейся этими операциями не приводится, так как дан-

ная программа нами не разрабатывалась, однако в приложении присутству-

ет распечатка этой программы.

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 - столбец решений

7D0X - столбец ошибки

Пояснения к тексту модуля:

Стандартные действия: инициализация начальных значений , цикл

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. Каждый из представленных методов по

своему хорош и применяется для систем определенного вида.