Исследование неявного метода Эйлера для линейной системы ОДУ с постоянным и переменным шагом (стр. 1 из 2)

Министерство Образования РФ

НГТУ

Кафедра экономической информатики

Курсовая работа по курсу

«Численные методы в экономике»:

Исследование неявного метода Эйлера для линейной системы

ОДУ с постоянным и переменным шагом

Введение

Целью данного проекта является исследование неявного метода Эйлера для решения линейных систем ОДУ с постоянным и переменным шагом. В ходе работы будет создана программа на языке матрично-ориентированной системы MatLAB и приведена математическая интерпретация метода. Также будет представлено влияние величины шага интегрирования и начальных значений на качество и точность вычислений. Будут проанализированы результаты и сделаны выводы.

1. Постановка задачи. Математическое описание метода

Для линейных систем

( Х = АХ+ВU(t) )

эта проблема решается очень просто. Нелинейные системы приходится линеаризовать в точке X^m, t_m, затем уже решать приведенным ниже способом:

Все значения в этой формуле, кроме X^m+1, которое нужно найти, известны (I – единичная матрица). Это получается линейная система, которая решается стандартными методами.

Неявный метод Эйлера:Группа неявных методов Рунге-Кутта используется для интегрирования "жестких" систем. Неявный метод Эйлера (Рунге-Кутта 1-го порядка) описывается с помощью следующей формулы:

X^m+1 = X^m+h_mF(X^m+1, t_m+1)

Как уже говорилось ранее, чтобы определить X^m+1 надо решить данную систему. При известных значениях величин X^m, h_m, t_m+1- это система нелинейных уравнений относительно X^m+1. Ее необходимо решать на каждом шаге по времени m.

Для линейных систем

( Х = АХ+ВU(t) )

эта проблема решается очень просто. Нелинейные системы приходится линеаризовать в точке X^m, t_m, затем уже решать приведенным ниже способом:

Все значения в этой формуле, кроме X^m+1, которое нужно найти, известны (I – единичная матрица). Это получается линейная система, которая решается стандартными методами.

Рассмотрим характеристики метода.

1.Точность.Ошибка аппроксимации по величине равна ошибке аппроксимации явного метода Эйлера, но она противоположна по знаку:

_i^am = -0.5h²_mX^..(t^-)

где h_m<= t^-<= t_m+1

2. Устойчивость метода

Сделав линеаризацию функции F(X,t)в точке X^m, h_m, t_m+1 получим уравнение относительно y^m+1

Y^m+1 = y^m+hy^m+1

Характеристического уравнение r-hlr-1=0 "дает" корень r=1/(1-hl).

1. Условие абсолютной устойчивости (Re(hl)<0): |1/(1-hl)|<1 или по другому

|1-hl|>1. Последнее неравенство можно преобразовать к виду:

[1-Re(h)]² + Im((h)]²>1

С учетом того, что мы рассматриваем ситуацию, когда Re(hl)<0, область абсолютной устойчивости, как следует из неравенства [1-Re(h)]² + Im((h)]²>1 - вся левая полуплоскость.

2. Условия относительной устойчивости

(Re(hl)>0): |1/(1-hl)|>1.

С учетом ограничения на скорость изменения приближенного решения относительного точного:

|1/(1-hl)|<|e^hl|

Из этого соотношения следует, что при |hl|®1 левая часть стремится к бесконечности. Это значит, что в правой полуплоскости имеется некоторая область, где неравенство

|1/(1-hl)|<|e^hl| не выполняется.

3. Выбор шага

4. Условия выбора шага диктуются требованиями абсолютной или относительной устойчивости. Однако область абсолютной устойчивости – вся левая полуплоскость. Поэтому шаг с этой точки зрения может быть любым.

Условия устойчивости жестче, так как есть область, где они могут быть нарушены. Можно показать, что эти условия будут выполнены, если в процессе решения задачи контролировать e^am_i , а шаг корректировать. Таким образом, шаг можно выбирать только из условий точности, при этом условия устойчивости будут соблюдены автоматически. Сначала задается допустимая погрешность аппроксимации:

e^{a gon} _i <=0.001| X_i | _max, i=1,n

Процедура выбора шага в процессе численного интегрирования состоит в следующем:

1. Решая систему

X^m+1 = X^m+h_mF(X^m+1, t_m+1

относительно X^m+1 с шагом h_m,

получаем очередную точку решения системы

X^o=F(X,t).

2. Оцениваем величину ошибки аппроксимации по формуле:

e^am _i = |h_m/(h_m+h_m+1)[(X_i^m+1 – X_i^m) – h_m/h_m-1(X^m_i – X^m-1_i]|

3. Действительное значение ошибки аппроксимации сравнивается с допустимым

e^am _i <e^{a gon} _i , i=1,n..

4. Если хотя бы для одного i неравенство не соблюдается, то шаг h_mуменьшается, в противном случае – увеличивается по сравнению с ранее принятым. Вопрос состоит в том как это сделать. Самый простой вариант – уменьшить в два раза. Вычисления повторяются с п.1.

5. Если упомянутые выше неравенства выполняются для всех i, то рассчитывается следующий шаг:

h_i^m+1 = ( e^{a gon} _i / |e^am _i |)*h^m, i=1,n

6. Шаг выбирается одинаковым для всех элементов вектора Х:

h_m+1=minh^m+1_i.

6. Вычисляется новый момент времени

7. t_m+2=t_m+1+h_m+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)&bsol;(y+h*B*U);

eam=abs(h*((yp1-y)-h*(y-ym1)/h1)/(h+h1));

if eam<=ep

yt=(I-A*h/2)&bsol;((I-A*h/2)&bsol;(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);

while (t<tfinal)

U=feval(funU, t);

if (t+h)>tfinal, h=tfinal-t; end

yt=(I-A*h/2)&bsol;((I-A*h/2)&bsol;(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);