Навигатор > Математика и Геометрия > Рефераты > Программа для решения системы нелинейных уравнений методом последовательной итерации...

Программа для решения системы нелинейных уравнений методом последовательной итерации обратной матрицы Якоби

	Скачать реферат: Программа для решения системы нелинейных уравнений методом последовательной итерации обратной матрицы Якоби

В экономике очень часто используется модель, называемая "черныйящик", то есть система у которой известны входы и выходы, а то, что происходит внутри - неизвестно. Законы по которым происходят изменения выходных сигналов в зависимости от входных могут быть различными, в том числе это могут быть и дифференциальные законы. Тогда встает проблема решения систем дифференциальных уравнений, которые в зависимости от своей структуры могут быть решены различными методами. Точное решение получить очень часто не удается, поэтому мы рассмотрим численные методы решения таких систем. Далее мы представим две группы методов: явные и неявные. Для решения систем ОДУ неявными методами придется решать системы нелинейных уравнений, поэтому придется ввести в рассмотрение группу методов решения систем нелинейных уравнений, которые в свою очередь будут представлены двумя методами. Далее следуют теоретические аспекты описанных методов, а затем будут представлены описания программ. Сами программы, а также их графики приведены в приложении.

Также стоит отметить, что в принципе все численные методы так или иначе сводятся к матричной алгебре, а в экономических задачах очень часто матрицы имеют слабую заполненность и большие размеры и поэтому неэффективно работать с полными матрицами. Одна из технологий, позволяющая разрешить данную проблему - технология разреженных матриц.

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

Таким образом мы рассмотрим весь спектр лабораторных работ. Описания всех программ приводятся после теоретической части. Все тексты программ и распечатки графиков приведены в приложении.

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-го порядка, то

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-го порядка состоит в следующем:

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.

Методы решения нелинейных САУ метод ньютона

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

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.

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

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

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

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

не обладает квадратичной сходимостью. Скорость сходимости можно увеличить, уменьшая 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 утверждения:

1. Хранить в памяти только ненулевые элементы.

2. В процессе решения принимать меры, уменьшающие возможность появления новых ненулевых элементов.

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

Разреженный строчный формат

Это одна из широко используемых схем для хранения разреженных матриц, которая предъявляет минимальные требования к памяти и очень удобна для выполнения основных операций с матрицами.

Значения ненулевых элементов и соответствующие столбцовые индексы хранятся по строкам в двух массивах: NA и JA. Для разметки строк в этих массивах используется массив указателей IA, отмечающих позиции массивов AN и JA, с которых начинается описание очередной строки. Последняя цифра в массиве IA содержит указатель первой свободной позиции в JA и AN. Рассмотрим конкретный пример, который будет также и далее применятся для демонстрации других операций и который был использован в качестве контрольного примера для программы, выполняющей основные операции с разреженными матрицами.

Данный способ представления называется полным (так как представлена вся матрица А) и упорядоченным (так как элементы каждой строки хранятся в соответствии с возрастанием столбцовых индексов). Обозначается 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

Символический этап.

Массив ANT можно получить, вставляя численное значение каждого

ННЭ, хранимое в AN, в вещественный список сразу после того, как соответствующее целое внесено в целый список. В итоге получим:

IAT = 1 3 5 6 9 11

JAT = 1 4 2 5 3 1 3 4 2 5

ANT = 3 5 4 7 8 5 2 6 1 9

Произведение разреженной матрицы и заполненного вектора-столбца

Алгоритм рассмотрим на нашем конкретном примере:

IAT = 1 3 5 7 9 11

JAT = 1 4 2 5 3 1 3 4 2 5

ANT = 3 5 4 7 8 5 2 6 1 9

B = ( -5 4 7 2 6 )

Обработка 1 строки:

Просматриваем массив IA и обнаруживаем, что первая строка матрицы

А соответствует первому и второму элементу массива JA: JA(1)=3, JA(2)=4, то есть ННЭ являются a 411 0 и a 414 0.

Просматриваем массив AN и устанавливаем, что a 411 0=3 и a 414 0=5.

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

b 44 0=2.

Вычисляем C 41 0= 3 * (-5) + 5 * 2 = -5.

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

C = (-5 58 56 1 58).

Произведение двух разреженных матриц

Рассмотрим метод на конкретном примере. Предположим, что нам необходимо перемножить две матрицы:

IA = 1 3 5 7 9 11

JA = 1 4 2 5 3 4 1 4 2 5

AN = 3 5 4 1 8 2 5 6 7 9

IAT = 1 3 5 7 9 11

JAT = 1 4 2 5 3 1 3 4 2 5

ANT = 3 5 4 7 8 5 2 6 1 9

Суть метода состоит в том, что используя метод переменного переключателя и расширенный вещественный накопитель, изменяется порядок накопления сумм для формирования элементов матрицы С. Мы будем формировать элементы C 4ij 0 не сразу, а постепенно накапливая, обращаясь только к строкам матрицы B.

Символический этап.

Определяем мерность IX = 5

IX = 0 0 0 0 0

1-я строка матрицы JAT: 1 4

JA(1) = 1 4 JC(1) = 1 4

IX = 1 0 0 1 0

JA(4) = 1 4

IX(1) = 1 ? Да. JC(1) - без изменений

IX(4) = 1 ? Да. JC(1) - без изменений

2-я строка матрицы JAT: 2 5

JA(2) = 2 5 JC(2) = 2 5

IX = 1 2 0 1 2

JA(5) = 2 5 -> JC(2) - без изменений

....

4-я строка матрицы JAT: 1 3 4

JA(1) = 1 4 JC(4) = 1 4

IX = 4 2 2 4 2

JA(3) = 3 4

IX(3) = 4 ? Нет. JC(4) = 1 4 3

IX(4) = 4 ? Да. JC(4) - без изменений

....

в итоге получим:

IC = 1 3 5 7 10 12

JC = 1 4 2 5 3 4 1 4 3 2 5

Численный этап.

X = 0 0 0 0 0

1-я строка: JC(1) = 1 4

AN(1) = 3 5, AA = 3

ANT(1) = 3 5

AA * ANT = 9 15

X = 9 0 0 15 0

AA = 5

ANT(1) = 3 5

AA * ANT = 15 25

X = 24 0 0 40 0

CN(1) = 24 40

....

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

IC: 1 3 5 7 10 12

JC: 1 4 2 5 3 4 1 4 3 2 5

CN: 24 40 20 35 80 0 55 22 66 16 144

Все вышеприведенные операции были получены на компьютере и результаты совпали для нашего контрольного примера. Описание программы

на языке 2 C 0, занимающейся этими операциями не приводится, так как данная программа нами не разрабатывалась, однако в приложении присутствует распечатка этой программы.

Практическая часть. описания программ.

1. Явные методы.

Данная группа представлена 3 программами, реализующими методы Эйлера,Рунге-Кутта 2-го и 4-го порядков. В данной группе все программы индентичны, поэтому далее следует описание программы, реализующем метод Эйлера, с указанием различий для методов Рунге-Кутта 2-го и 4-го порядков.

EILER.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 порядка.

Также как и в вышеприведенном случае, будет описан метод Эйлера, а отличия метода Рунге-Кутта будут отмечены в скобках.

NME.M

Головной модуль.

Входные и выходные данные отсутствуют.

Язык реализации: PC MathLab

Операционная система: MS-DOS 3.30 or higher

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

Выполняет стандартные действия: очистка экрана, инициализация и ввод начальных значений, вызов подпрограмм обработки и вычислений, а также построение графиков.

NMEF.M, NRG2.M

Вычислительные модули.

Входные данные:

T0, Tfinal - начальные и конечные моменты времени

X0 - вектор-столбец начальных значений.

H - начальный шаг

A - матрица, на диагонали которой стоят собственные числа линеаризованной системы ОДУ.

Выходные данные:

T - столбец времени

X - столбец решений

7D 0X - столбец ошибки

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

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

While T < Tfinal, вычисление по формулам, вывод промежуточных результатов, формирование выходных значений массивов.

3. Методы решения нелинейных САУ

Представлены группой из 4-х методов: метод последовательных приближений, метод Ньютона, метод Ньютона дискретный, метод продолжения решения по параметру.

Метод последовательных приближений.

MMPP.M

Головной модуль.

Входные и выходные данные отсутствуют.

Язык реализации: PC MathLab

Операционная система: MS-DOS 3.30 or higher

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

Очистка экрана, инициализация начальных значений, вызов вычислительного модуля MPP.M, вывод результатов в виде графиков.

MPP.M

Вычислительный модуль.

Входные данные:

X0 - начальное приближение в виде вектора-строки

Fun1 - функция, возвращающая вычисленные левые части

Fun2 - функция, возвращающая матрицу Якоби в определенной точке.

E - допустимая ошибка.

Выходные данные:

Mout - номера итераций

Xout - приближения на каждой итерации

DXout - ошибка на каждой итерации

Язык реализации: PC MathLab

Операционная система: MS-DOS 3.30 or higher

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

Аналогичен вышеприведенным вычислительным модулям - инициализация начальных значений, вычисления по формулам, вывод промежуточных результатов, формирование выходных значений. По мере необходимости вызывает подпрограммы Fun1 и Fun2.

Методы Ньютона и Ньютона дискретный

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

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

Стандартные вычислительные модули, производящие соответствующие им действия. Отличие их в том, что в первом случае для вычисления матрицы Якоби вызывается подпрограмма, а во втором случае матрица Якоби вычисляется внутри модуля.

Метод продолжения решения по параметру

MMPRPP.M

Головной модуль.

Входные и выходные данные отсутствуют.

Язык реализации: PC MathLab

Операционная система: MS-DOS 3.30 or higher

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

Стандартный головной модуль со всеми вытекающими отсюда последствиями.

MPRPP.M

Вычислительный модуль.

Входные данные:

Fun1 - имя подпрограммы, вычисляющей правые части

Fun2 - имя подпрограммы, вычисляющем матрицу Якоби

X0 - начальное приближение

dT - начальный шаг

Edop - допустимая ошибка

Trace - вывод на экран

Язык реализации: PC MathLab

Операционная система: MS-DOS 3.30 or higher

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

Стандартный вычислительный модуль - инициализация, вычисление, вывод, формирование результата. Стоит отметить, что поскольку метод имеет глобальную сходимость, то объем вычислений тут значительный, а это значит, что при выполнении вычислений требуется значительное количество свободной оперативной памяти.

Выводы

При выполнении данных лабораторных работ были изучены основные численные методы для решения ОДУ, САУ, а также технология разреженных матриц. Заодно были получены основные навыки в программировании математической системы PC MathLab. Каждый из представленных методов по своему хорош и применяется для систем определенного вида.

Поиск