Численные методы в экономике
2ВВЕДЕНИЕ
В экономике очень часто используется модель, называемая "черный
ящик", то есть система у которой известны входы и выходы, а то, что
происходит внутри - неизвестно. Законы по которым происходят изменения
выходных сигналов в зависимости от входных могут быть различными, в
том числе это могут быть и дифференциальные законы. Тогда встает проб-
лема решения систем дифференциальных уравнений, которые в зависимости
от своей структуры могут быть решены различными методами. Точное реше-
ние получить очень часто не удается, поэтому мы рассмотрим численные
методы решения таких систем. Далее мы представим две группы методов:
явные и неявные. Для решения систем ОДУ неявными методами придется ре-
шать системы нелинейных уравнений, поэтому придется ввести в рассмот-
рение группу методов решения систем нелинейных уравнений, которые в
свою очередь будут представлены двумя методами. Далее следуют теорети-
ческие аспекты описанных методов, а затем будут представлены описания
программ. Сами программы, а также их графики приведены в приложении.
Также стоит отметить, что в принципе все численные методы так или
иначе сводятся к матричной алгебре, а в экономических задачах очень
часто матрицы имеют слабую заполненность и большие размеры и поэтому
неэффективно работать с полными матрицами. Одна из технологий, позво-
ляющая разрешить данную проблему - технология разреженных матриц.
В связи с этим, мы рассмотрим данную технологию и операции умножения и
транспонирования над такими матрицами.
Таким образом мы рассмотрим весь спектр лабораторных работ. Опи-
сания всех программ приводятся после теоретической части. Все тексты
программ и распечатки
графиков приведены в приложении.
2ТЕОРЕТИЧЕСКАЯ ЧАСТЬ
1. ОПИСАНИЕ МЕТОДОВ ИНТЕГРИРОВАНИЯ СИСТЕМ ОДУ
И ИХ ХАРАКТЕРИСТИК
ЯВНЫЙ МЕТОД ЭЙЛЕРА И ЕГО ХАРАКТЕРИСТИКИ
Алгоритм этого метода определяется формулой:
x 5m+1 0 = x 5m 0 + h 4m 0*F(x 5m 0, t 4m 0) 4, 0
которая получается путём аппроксимации ряда Тейлора до членов перво-
го порядка производной x'(t 4m 0),т.к. порядок точности равен 1 (s=1).
Для аналитического исследования свойств метода Эйлера линеари-
зуется исходная система ОДУ 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 может считаться констан-
той. Ввиду того ,что для линейной системы свойство устойчивости за-
висит лишь от А, то входной сигнал в системе (2) не показан. Элемен-
ты матрицы А меняются с изменением точки линеаризации,т.е. с измене-
нием m.
Характеристики метода:
1. _Численная устойчивость ..
Приведя матрицу А к диагональной форме : А = Р* 7l 0*Р 5-1 0 и перейдя
к новым переменным y = P 5-1 0*x , система (2) примет вид :
y' = 7l 0*y
Тогда метод Эйлера для уравнения (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,
где 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
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,
Ошибка аппроксимации Е 5a 0 = k*h 4m 53 0 .
Алгоритм метода Рунге-Кутта 4-го порядка
x 5m+1 0=x 5m 0+h 4m 0/6(k 41
где k 41 0=f(x 5m 0,t 4m 0);
k 42 0=f(x 5m 0+h 4m 0/2*k 41
k 44 0=f(x 5m 0+h 4m 0*k 43 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,
где 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
3. Действительное значение аппроксимации сравнивается с до-
пустимым. Если Е 4i 5am 0 < E 4i 5aдоп 0, то выполняется следующий пункт, в про-
тивном случае шаг уменьшается в два раза , и вычисления повторяются
с п.1.
4. Рассчитываем следующий шаг:
h 4i 5m+1 0 = SQR(E 4i 5aдоп 7/2 0Е 4i 5am
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
Ошибка аппроксимации Е 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
где
k 41 0=f(x 5m+1 0,t 4m+1 0);
k 43 0=f(x 5m+1 0-h 4m 0/2*k
Ошибка аппроксимации Е 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
--- = ------------------------------
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)
Нетрудно проверить, что для этих вектор-функций выполняются усло-
- Численные методы, использованные в методах оценки недвижимости
- Численные методы решения систем линейных алгебраических уравнений
- Численные методы решения систем линейных уравнений
- Численость и структура населения
- Числовая последовательность
- Числовая последовательность
- Числові ряди Збіжність і розбіжність Сума ряду Дії над збіжними рядами Необхідна ознака збіж
- Численность и структура населения
- Численность и структура населения
- Численность и структура населения, тенденции их изменения
- Численность населения, особенности его географического распределения
- Численность населения России
- Численность населения России
- Численные методы
