Министерство Образования РФ
НГТУ
Кафедра экономической информатики
Курсовая работа по курсу
«Численные методы в экономике»:
Исследование неявного метода Эйлера для линейной системы
ОДУ с постоянным и переменным шагом
Введение
Целью данного проекта является исследование неявного метода Эйлера для решения линейных систем ОДУ с постоянным и переменным шагом. В ходе работы будет создана программа на языке матрично-ориентированной системы MatLAB и приведена математическая интерпретация метода. Также будет представлено влияние величины шага интегрирования и начальных значений на качество и точность вычислений. Будут проанализированы результаты и сделаны выводы.
1.
Постановка задачи. Математическое описание метода
Для линейных систем
( Х = АХ+В
U
(
t
)
)
эта проблема решается очень просто. Нелинейные системы приходится линеаризовать в точке Xm
, tm
, затем уже решать приведенным ниже способом:
Все значения в этой формуле, кроме Xm
+1
, которое нужно найти, известны (I – единичная матрица). Это получается линейная система, которая решается стандартными методами.
Неявный метод Эйлера:
Группа неявных методов Рунге-Кутта используется для интегрирования "жестких" систем. Неявный метод Эйлера (Рунге-Кутта 1-го порядка) описывается с помощью следующей формулы:
Xm
+1
=
Xm
+
hm
F
(
Xm
+1
,
tm
+1
)
Как уже говорилось ранее, чтобы определить Xm
+1
надо решить данную систему. При известных значениях величин Xm
, hm
, tm
+1
- это система нелинейных уравнений относительно Xm
+1
.
Ее необходимо решать на каждом шаге по времени m.
Для линейных систем
( Х = АХ+В
U
(
t
)
)
эта проблема решается очень просто. Нелинейные системы приходится линеаризовать в точке Xm
, tm
, затем уже решать приведенным ниже способом:
Все значения в этой формуле, кроме Xm
+1
, которое нужно найти, известны (I – единичная матрица). Это получается линейная система, которая решается стандартными методами.
Рассмотрим характеристики метода.
1.Точность.
Ошибка аппроксимации по величине равна ошибке аппроксимации явного метода Эйлера, но она противоположна по знаку:
i
am
= -0.5h2
m
X..
(t-
)
где hm
<= t-
<= tm
+1
2.
Устойчивость метода
Сделав линеаризацию функции F(X,t)в точке Xm
, hm
, tm
+
1 получим уравнение относительно ym
+1
Ym+1
= ym
+h
ym+1
Характеристического уравнение r-hlr-1=0 "дает" корень r=1/(1-hl).
1. Условие абсолютной устойчивости (Re(hl)<0): |1/(1-hl)|<1 или по другому
|1-hl|>1. Последнее неравенство можно преобразовать к виду:
[1-Re(h
)]2
+ Im((h
)]2
>1
С учетом того, что мы рассматриваем ситуацию, когда Re(hl)<0, область абсолютной устойчивости, как следует из неравенства [1-
Re
(
h
)]2
+
Im
((
h
)]2
>1
- вся левая полуплоскость.
2. Условия относительной устойчивости
(Re(hl)>0): |1/(1-hl)|>1.
С учетом ограничения на скорость изменения приближенного решения относительного точного:
|1/(1-
h
l
)|<|
eh
l
|
Из этого соотношения следует, что при |hl|®1 левая часть стремится к бесконечности. Это значит, что в правой полуплоскости имеется некоторая область, где неравенство
|1/(1-
h
l
)|<|
eh
l
|
не выполняется.
3.
Выбор шага
4. Условия выбора шага диктуются требованиями абсолютной или относительной устойчивости. Однако область абсолютной устойчивости – вся левая полуплоскость. Поэтому шаг с этой точки зрения может быть любым.
Условия устойчивости жестче, так как есть область, где они могут быть нарушены. Можно показать, что эти условия будут выполнены, если в процессе решения задачи контролировать ea
m
i
, а шаг корректировать. Таким образом, шаг можно выбирать только из условий точности, при этом условия устойчивости будут соблюдены автоматически. Сначала задается допустимая погрешность аппроксимации:
e
a
gon
i
<=0.001| Xi
| max,
i=1,n
Процедура выбора шага в процессе численного интегрирования состоит в следующем:
1. Решая систему
Xm
+1
=
Xm
+
hm
F
(
Xm
+1
,
tm
+1
относительно Xm
+1
с шагом hm
,
получаем очередную точку решения системы
Xo
=F(X,t).
2.
Оцениваем величину ошибки аппроксимации по формуле:
e
a
m
i
= |hm
/(hm
+hm+1
)[(Xi
m+1
– Xi
m
) – hm
/hm-1
(Xm
i
– Xm-1
i
]|
3. Действительное значение ошибки аппроксимации сравнивается с допустимым
e
a
m
i
<
e
a
gon
i
, i=1,n..
4. Если хотя бы для одного i неравенство не соблюдается, то шаг hm
уменьшается, в противном случае – увеличивается по сравнению с ранее принятым. Вопрос состоит в том как это сделать. Самый простой вариант – уменьшить в два раза. Вычисления повторяются с п.1.
5. Если упомянутые выше неравенства выполняются для всех i, то рассчитывается следующий шаг:
hi
m+1
=
(
e
a
gon
i
/ |
e
a
m
i
|)*hm
, i=1,n
6. Шаг выбирается одинаковым для всех элементов вектора Х:
hm
+1
=
min
hm
+1
i
.
6. Вычисляется новый момент времени
7. tm
+2
=
tm
+1
+
hm
+1
и алгоритм повторяется с п.1.
Интегрирование «жесткой» системы ОДУ:
«Жесткими» называются системы, число обусловленности которых:
В решении системы присутствуют экспоненты, сильно отличающиеся друг от друга по скорости затухания. После затухания быстрой экспоненты наступает медленная часть решения, казалось бы шаг интегрирования можно увеличить. Однако требования устойчивости должны выполняться на всем участке решения. Это требует больших затрат машинного времени. Таким образом, ограниченность области абсолютной устойчивости делает явный метод Рунге-Кутта 4-го порядка непригодным для интегрирования «жестких» линейных систем ОДУ.
4. Описание программного обеспечения
Общие сведения и требования к ПО
:
Программа состоит из 3-х файлов, написанных на языке MatLAB - rkpost1.m (для метода с постоянным шагом), rkper1.m(для метода с переменным шагом) и тестовая программа для них (test.m).
Функциональное назначение:
Программа предназначена для решения линейных систем ОДУ в среде MatLABметодом численного интегрирования Эйлера (Рунге-Кутта 1-го порядка).
Программа:
- rkper1.m
function [tout, yout, eout]=rkper1(funA, funB, funU, t0, tfinal, y0, ep, trace)
if nargin<7, ep=0.001; end
if nargin<8, trace=1; end
t=t0; y=y0;
tout=t; yout=y.';
h1=(tfinal-t)/20000;
h=h1*200;
if trace
clc, t, h, y
end
A=feval(funA);
B=feval(funB);
n=ones(max(size(y0)),1);
I=diag(n,0);
ym1=y0; yp1=y0;
while (t<tfinal)
U=feval(funU, t+h);
if (t+h)>tfinal, h=tfinal-t; end
yp1=(I-A*h)(y+h*B*U);
eam=abs(h*((yp1-y)-h*(y-ym1)/h1)/(h+h1));
if eam<=ep
yt=(I-A*h/2)((I-A*h/2)(y+h/2*B*U)+h/2*B*U);
h1=h;
ym1=y;
y=yp1;
eout=[eout;abs(y-yt).'];
tout=[tout;t]; yout=[yout;y.'];
h=min(sqrt((n*ep)./eam)*h);
t=t+h;
else h=h/2;
end
if trace
home, t, h, y
end
end
- rkpost1.m
function [tout, yout, eout]=rkpost1(funA, funB, funU, t0, tfinal, y0, h, trace)
if nargin<7, h=(tfinal-t0)/5; end
if nargin<8, trace=1; end
t=t0; y=y0;
tout=t; yout=y.';
if trace
clc, t, h, y
end
A=feval(funA);
B=feval(funB);
I=diag(ones(1,max(size(y0))),0);
whi
U=feval(funU, t);
if (t+h)>tfinal, h=tfinal-t; end
yt=(I-A*h/2)((I-A*h/2)(y+h/2*B*U)+h/2*B*U);
t=t+h;
eout=[eout;abs(y-yt).'];
tout=[tout;t]; yout=[yout;y.'];
if trace
home, t, h, y
end
end
- test.m
disp('Решаем нежесткую систему:')
pause
disp('Решаем переменный шаг:')
pause
% Переменный шаг
[t1,y1,e1]=rkper1 ('a','b','u',0,3.5,[0.1;0.1],0.01);
[t2,y2,e2]= rkper1 ('a','b','u',0,3.5,[0.5;0.5],0.01);
[t3,y3,e3]= rkper1 ('a','b','u',0,3.5,[1;1],0.01);
plot(t1,y1,t2,y2,t3,y3)
pause
t1e=t1(1:max(size(t1))-1);
t2e=t2(1:max(size(t2))-1);
t3e=t3(1:max(size(t3))-1);
plot(t1e,e1,t2e,e2,t3e,e3)
pause
disp('Решаем постоянный шаг:')
pause;
% Постоянный шаг
[tc1,yc1,ec1]=nrk1('a','b','u',0,3.5,[0.2;0.2],0.1);
[tc2,yc2,ec2]=nrk1('a','b','u',0,3.5,[0.2;0.2],0.01);
[tc3,yc3,ec3]=nrk1('a','b','u',0,3.5,[0.2;0.2],0.005);
plot(tc1,yc1,tc2,yc2,tc3,yc3)
pause
t1ec=tc1(1:max(size(tc1))-1);
t2ec=tc2(1:max(size(tc2))-1);
t3ec=tc3(1:max(size(tc3))-1);
plot(t1ec,ec1,t2ec,ec2,t3ec,ec3)
pause
[tc1,yc1,ec1]= rkpost1 ('a','b','u',0,3.5,[0.1;0.1],0.1);
[tc2,yc2,ec2]= rkpost1 ('a','b','u',0,3.5,[0.5;0.5],0.1);
[tc3,yc3,ec3]= rkpost1 ('a','b','u',0,3.5,[1;1],0.1);
plot(tc1,yc1,tc2,yc2,tc3,yc3)
pause
t1ec=tc1(1:max(size(tc1))-1);
t2ec=tc2(1:max(size(tc2))-1);
t3ec=tc3(1:max(size(tc3))-1);
plot(t1ec,ec1,t2ec,ec2,t3ec,ec3)
pause
disp('Решаем жесткую систему:')
pause
disp('Решаем переменный шаг:')
pause
% Переменный шаг
[t1,y1,e1]=nrk1var('a1','b','u',0,3.5,[0.1;0.1],0.01);
[t2,y2,e2]=nrk1var('a1','b','u',0,3.5,[0.5;0.5],0.01);
[t3,y3,e3]=nrk1var('a1','b','u',0,3.5,[1;1],0.01);
plot(t1,y1,t2,y2,t3,y3)
pause
t1e=t1(1:max(size(t1))-1);
t2e=t2(1:max(size(t2))-1);
t3e=t3(1:max(size(t3))-1);
plot(t1e,e1,t2e,e2,t3e,e3)
pause
disp('Решаем постоянный шаг:')
pause;
% Постоянный шаг
[tc1,yc1,ec1]=nrk1('a1','b','u',0,3.5,[0.2;0.2],0.1);
[tc2,yc2,ec2]=nrk1('a1','b','u',0,3.5,[0.2;0.2],0.01);
[tc3,yc3,ec3]=nrk1('a1','b','u',0,3.5,[0.2;0.2],0.005);
plot(tc1,yc1,tc2,yc2,tc3,yc3)
pause
t1ec=tc1(1:max(size(tc1))-1);
t2ec=tc2(1:max(size(tc2))-1);
t3ec=tc3(1:max(size(tc3))-1);
plot(t1ec,ec1,t2ec,ec2,t3ec,ec3)
pause
[tc1,yc1,ec1]=nrk1('a1','b','u',0,3.5,[0.1;0.1],0.1);
[tc2,yc2,ec2]=nrk1('a1','b','u',0,3.5,[0.5;0.5],0.1);
[tc3,yc3,ec3]=nrk1('a1','b','u',0,3.5,[1;1],0.1);
plot(tc1,yc1,tc2,yc2,tc3,yc3)
pause
t1ec=tc1(1:max(size(tc1))-1);
t2ec=tc2(1:max(size(tc2))-1);
t3ec=tc3(1:max(size(tc3))-1);
plot(t1ec,ec1,t2ec,ec2,t3ec,ec3)
pause
- Нежесткая матрица А:
function A=a();
A=[-5/6 1/3;
1/3 -1/3];
Жесткая матрица А:
function A=a();
A=[-50 50;
50 -50.1];
Нежесткая матрица А:
function B=b();
B=[5/2; 0];
Матрица U:
function U=u(t);
U=[1];
Переменные
- Выходные переменные:
tout
– выходной вектор времени;
yout
– выходное значение функций (соответствует вектору времени);
eout
– выходное значение ошибок определения функций (высчитано по методу Рунге).
- Внутренние переменные:
A
,
B
,
U
– соответствующие матрицы (вектора) из формулы Х = АХ+В
U
(
t
)
.
t
– текущее время;
y
– текущие значения функций;
ym1
– предыдущие значения функций;
y
p
1
– следующие значения функций;
y
t
–значения функций, высчитанное с половинчатым шагом (для определения ошибки по методу Рунге);
h
– текущий шаг;
h1
– предыдущий шаг;
n
– вектор с единицами с размерностями как y;
I
– единичная матрица;
eam
– текущая ошибка, вычисляемая для определения допустимого шага;
- Входы:
funA, funB, funU
– внешние функции вычисления A, B и U;
t0, tfinal
– начальное и конечное значение времени;
y0
– начальное значение функций (с которых начинаем решение системы ОДУ);
ep
– допустимая ошибка;
trace
– выводить или нет на экран промежуточные значения;
h
– шаг в методе с постоянным шагом.
Описание алгоритма:
Сначала производятся все начальные инициализации и установки. Затем запускается основной цикл расчета значений функций. В нем для переменного шага высчитывается критерий. Затем по этому критерию определяется сам шаг. Цикл повторяется. Для постоянного шага производится непосредственно пересчет значения функции. Также, для определения ошибки по формуле Рунге, производится вычисление функции через два половинных шага. На выходе получается значения времени, функций и ошибок. "Длина" массива ошибок на 1 меньше, чем "длины" других массивов (в силу специфики его вычисления).
5. Описание тестовых задач
В программе тестирования сначала находится решение нежесткой системы с разными начальными значениями y, шагами и способами: с переменным и постоянным шагами. Строятся совокупные графики функций и ошибок для каждого случая. Затем решается нежесткая система с разными допустимыми ошибками и разными начальными значениями y и разными шагами. Строятся совокупные графики функций и ошибок для каждого случая. Для способа с переменным шагом изменяем только начальные значения, для способа с постоянным шагом – начальные значения, и величины шага.
График1:
График2.
График 3.
График 4.
График5.
График 6.
График 7.
График 8.
График 9.
График 10.
эйлер линейный программа интегрирование
График 11.
График 12.
6. Анализ результатов. Выводы
Ниже показаны полученные графики функций и ошибок:
График1:
Var, XX0=[0.1;0.1] X0=[0.5;0.5] X0=[1;1]
График2
:
Var, EX0=[0.1;0.1] X0=[0.5;0.5] X0=[1;1]
График
3:
Const, X h=0.1, X0=[0.2;0.2] h=0.01, X0=[0.2;0.2] h=0.005, X0=[0.2;0.2]
График4
:
Const, E h=0.1, X0=[0.2;0.2] h=0.01, X0=[0.2;0.2] h=0.005, X0=[0.2;0.2]
График
5
:
Const, X h=0.1, X0=[0.1;0.1] h=0.1, X0=[0.5;0.5] h=0.1, X0=[1;1]
График6:
Const, E h=0.1, X0=[0.1;0.1] h=0.1, X0=[0.5;0.5] h=0.1, X0=[1;1]
Пояснения
:
Var
– метод с переменным шагом;
Const
– метод с постоянным шагом.
X0=[0.1;0.1], X0=[0.5;0.5], X0=[1;1] отражают цветовую гамму на графике.
X
– графики функций;
E
– графики ошибок.
h
– шаг, X
0
– начальное значение X.
Рисунки с 7 -12 аналогичны рисункам 1-6 (соответственно) по начальным данным и отличаются лишь жесткостью решаемой системы.
Выводы
Видим, что при увеличении начального значения происходит лишь параллельный сдвиг решения, однако график ошибки перемещается в противоположную сторону (при увеличении начального значения она уменьшается и наоборот), причем непропорционально. Следовательно, при увеличении начального значения X ошибка уменьшается, однако с увеличением t это изменение становится все менее и менее заметным. Это характерно как для метода с переменным шагом, так и для метода с постоянным шагом. Кроме того, для метода с постоянным шагом увеличение шага приводит к уменьшению точности решения. Что касается жесткой системы, то для нее характерно прямолинейное конечное решение. С уменьшением шага ошибка, посчитанная по методу Рунге, уменьшается. При увеличении начального X происходит параллельный сдвиг графика конечного решения. Если в случае нежесткой системы графики решений разных составляющих X находятся далеко друг от друга, то при решении жесткой системы получаются очень близкие графики.
Заключение
В ходе выполнения работы был изучен неявный метод Эйлера для решения линейных систем ОДУ. В ходе реализации проекта и проведения тестирования была проверена справедливость теоретических выкладок. Получены сведения о зависимости точности интегрирования от величины и способа выбора шага, а также от начальных значений переменных.
Используемая литература
1. Сарычева О.М. Численные методы. – Новосибирск, 1995г. – 65с.
2. Бахвалов Н.С. Численные методы. Ч1.- М: Наука, 1975г. – 632с., илл.
3. Копченова Н.В., Марон И.А. Вычислительная математика в примерах и задачах . – М: Наука, 1972г. - 368с