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