МинистерствоОбразования РФ
НГТУ
Кафедраэкономической информатики
Курсовая работа по курсу«Численные методы в экономике»:
Исследование неявного метода Эйлерадля линейной системы
ОДУ с постоянным и переменным шагом
Введение
Целью данного проектаявляется исследование неявного метода Эйлера для решения линейных систем ОДУ спостоянным и переменным шагом. В ходе работы будет создана программа на языкематрично-ориентированной системы Mat LAB и приведенаматематическая интерпретация метода. Также будет представлено влияние величинышага интегрирования и начальных значений на качество и точность вычислений.Будут проанализированы результаты и сделаны выводы.
1. Постановказадачи. Математическое описание метода
Для линейных систем
( Х = АХ+ВU(t) )
эта проблема решаетсяочень просто. Нелинейные системы приходится линеаризовать в точке Xm, tm, затем уже решать приведенным нижеспособом:
/>
Все значения в этойформуле, кроме Xm+1, которое нужно найти, известны (I – единичная матрица). Это получаетсялинейная система, которая решается стандартными методами.
Неявный метод Эйлера: Группа неявных методовРунге-Кутта используется для интегрирования «жестких» систем. Неявныйметод Эйлера (Рунге-Кутта 1-го порядка) описывается с помощью следующей формулы:
Xm+1= Xm+hmF(Xm+1, tm+1)
Какуже говорилось ранее, чтобы определить Xm+1 надо решить даннуюсистему. При известных значениях величин Xm, hm, tm+1 — это системанелинейных уравнений относительно Xm+1. Ее необходиморешать на каждом шаге по времени m.
Для линейных систем
( Х = АХ+ВU(t) )
эта проблема решаетсяочень просто. Нелинейные системы приходится линеаризовать в точке Xm, tm, затем уже решать приведенным нижеспособом:
/>
Все значения в этойформуле, кроме Xm+1, которое нужно найти, известны (I – единичная матрица). Это получаетсялинейная система, которая решается стандартными методами.
Рассмотримхарактеристики метода.
1.Точность.Ошибка аппроксимации по величине равнаошибке аппроксимации явного метода Эйлера, но она противоположна по знаку:
iam=-0.5h2mX..(t-)
гдеhm
2. Устойчивость метода
Сделавлинеаризацию функции F(X,t)вточке Xm, hm, tm+1 получимуравнение относительно ym+1
Ym+1 = ym+hym+1
Характеристическогоуравнение r-hlr-1=0 «дает» корень r=1/(1-hl).
1. Условие абсолютной устойчивости (Re(hl)
|1-hl|>1. Последнее неравенство можно преобразовать квиду:
[1-Re(h)]2+ Im((h)]2>1
/>
Сучетом того, что мы рассматриваем ситуацию, когда Re(hl)[1-Re(h)]2+ Im((h)]2>1 — вся левая полуплоскость.
2.Условия относительной устойчивости
(Re(hl)>0): |1/(1-hl)|>1.
Сучетом ограничения на скорость изменения приближенного решения относительноготочного:
|1/(1-hl)|ehl|
Изэтого соотношения следует, что при |hl|®1 левая часть стремится к бесконечности.Это значит, что в правой полуплоскости имеется некоторая область, где неравенство
|1/(1-hl)|ehl| не выполняется.
3. Выбор шага
4. Условия выборашага диктуются требованиями абсолютной или относительной устойчивости. Однакообласть абсолютной устойчивости – вся левая полуплоскость. Поэтому шаг с этойточки зрения может быть любым.
Условияустойчивости жестче, так как есть область, где они могут быть нарушены. Можнопоказать, что эти условия будут выполнены, если в процессе решения задачиконтролировать eami, а шагкорректировать. Таким образом, шаг можно выбирать только из условий точности,при этом условия устойчивости будут соблюдены автоматически. Сначала задаетсядопустимая погрешность аппроксимации:
eagoni i |max, i=1,n
Процедуравыбора шага в процессе численного интегрирования состоит в следующем:
1. Решая систему
Xm+1= Xm+hmF(Xm+1, tm+1
относительноXm+1 с шагом hm,
получаемочередную точку решения системы
Xo=F(X,t).
2. Оцениваемвеличину ошибки аппроксимации по формуле:
eam i=|hm/(hm+hm+1)[(Xim+1 – Xim)– hm/hm-1(Xmi – Xm-1i]|
3. Действительное значение ошибки аппроксимации сравнивается с допустимым
eam ieagoni, i=1,n..
4. Если хотя бы для одного i неравенство не соблюдается, то шаг hmуменьшается, впротивном случае – увеличивается по сравнению с ранее принятым. Вопрос состоитв том как это сделать. Самый простой вариант – уменьшить в два раза. Вычисленияповторяются с п.1.
5. Если упомянутые выше неравенства выполняются для всех i, торассчитывается следующий шаг:
him+1 = ( eagoni/ |eam i|)*hm, i=1,n
6.Шаг выбирается одинаковым для всех элементов вектора Х:
hm+1=minhm+1i.
6. Вычисляется новый момент времени
7. tm+2=tm+1+hm+1и алгоритм повторяетсяс п.1.
Интегрирование«жесткой» системы ОДУ:«Жесткими» называются системы, число обусловленности которых:
/>
В решении системыприсутствуют экспоненты, сильно отличающиеся друг от друга по скоростизатухания. После затухания быстрой экспоненты наступает медленная частьрешения, казалось бы шаг интегрирования можно увеличить. Однако требованияустойчивости должны выполняться на всем участке решения. Это требует большихзатрат машинного времени. Таким образом, ограниченность области абсолютнойустойчивости делает явный метод Рунге-Кутта 4-го порядка непригодным дляинтегрирования «жестких» линейных систем ОДУ.
4. Описание программного обеспечения
Общие сведения итребования к ПО: Программасостоит из 3-х файлов, написанных на языке Mat LAB — 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
if nargin
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
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
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
if nargin
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);
while (t
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/61/3;
1/3-1/3];
Жесткаяматрица А:
function A=a();
A=[-5050;
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 – предыдущие значения функций;
yp1 – следующие значения функций;
yt –значения функций, высчитанное споловинчатым шагом (для определения ошибки по методу Рунге);
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, X X0=[0.1;0.1] X0=[0.5;0.5] X0=[1;1]
График2:
Var, E X0=[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– начальное значение X.
Рисунки с 7 -12аналогичны рисункам 1-6 (соответственно) по начальным данным и отличаются лишьжесткостью решаемой системы.
Выводы
Видим, что при увеличенииначального значения происходит лишь параллельный сдвиг решения, однако графикошибки перемещается в противоположную сторону (при увеличении начальногозначения она уменьшается и наоборот), причем непропорционально. Следовательно,при увеличении начального значения X ошибка уменьшается, однако с увеличением t это изменение становится все менее и менее заметным. Этохарактерно как для метода с переменным шагом, так и для метода с постояннымшагом. Кроме того, для метода с постоянным шагом увеличение шага приводит куменьшению точности решения. Что касается жесткой системы, то для неехарактерно прямолинейное конечное решение. С уменьшением шага ошибка, посчитаннаяпо методу Рунге, уменьшается. При увеличении начального X происходит параллельный сдвигграфика конечного решения. Если в случае нежесткой системы графики решенийразных составляющих X находятся далекодруг от друга, то при решении жесткой системы получаются очень близкие графики.
Заключение
В ходе выполнения работыбыл изучен неявный метод Эйлера для решения линейных систем ОДУ. В ходереализации проекта и проведения тестирования была проверена справедливостьтеоретических выкладок. Получены сведения о зависимости точности интегрированияот величины и способа выбора шага, а также от начальных значений переменных.
Используемая литература
1. Сарычева О.М. Численные методы. –Новосибирск, 1995г. – 65с.
2. Бахвалов Н.С. Численные методы.Ч1.- М: Наука, 1975г. – 632с., илл.
3. Копченова Н.В., Марон И.А.Вычислительная математика в примерах и задачах. – М: Наука, 1972г. — 368с