Исследование неявного метода Эйлера для линейной системы ОДУ с постоянным и переменным шагом
Министерство Образования РФ
НГТУ
Кафедра экономической информатики
Курсовая работа по курсу
ВлЧисленные методы в экономикеВ»:
Исследование неявного метода Эйлера для линейной системы
ОДУ с постоянным и переменным шагом
Введение
Целью данного проекта является исследование неявного метода Эйлера для решения линейных систем ОДУ с постоянным и переменным шагом. В ходе работы будет создана программа на языке матрично-ориентированной системы 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<= t-<= tm+1
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)<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-hl)|<|ehl
Из этого соотношения следует, что при |hl|Во1 левая часть стремится к бесконечности. Это значит, что в правой полуплоскости имеется некоторая область, где неравенство
|1/(1-hl)|<|ehlне выполняется.
3. Выбор шага
4. Условия выбора шага диктуются требованиями абсолютной или относительной устойчивости. Однако область абсолютной устойчивости тАУ вся левая полуплоскость. Поэтому шаг с этой точки зрения может быть любым.
Условия устойчивости жестче, так как есть область, где они могут быть нарушены. Можно показать, что эти условия будут выполнены, если в процессе решения задачи контролировать eami , а шаг корректировать. Таким образом, шаг можно выбирать только из условий точности, при этом условия устойчивости будут соблюдены автоматически. Сначала задается допустимая погрешность аппроксимации:
ea gon i <=0.001| Xi | 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 i <ea gon i , i=1,n.
4. Если хотя бы для одного i неравенство не соблюдается, то шаг hmуменьшается, в противном случае тАУ увеличивается по сравнению с ранее принятым. Вопрос состоит в том как это сделать. Самый простой вариант тАУ уменьшить в два раза. Вычисления повторяются с п.1.
5. Если упомянутые выше неравенства выполняются для всех i, то рассчитывается следующий шаг:
him+1 = ( ea gon i / |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<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 (tU=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);
while (tU=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 тАУ предыдущие значения функций;
y1 тАУ следующие значения функций;
yt тАУзначения функций, высчитанное с половинчатым шагом (для определения ошибки по методу Рунге);
h тАУ текущий шаг;
h1 тАУ предыдущий шаг;
тАУ вектор с единицами с размерностями как 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 тАУ шаг, X0 тАУ начальное значение X.
Рисунки с 7 -12 аналогичны рисункам 1-6 (соответственно) по начальным данным и отличаются лишь жесткостью решаемой системы.
Выводы
Видим, что при увеличении начального значения происходит лишь параллельный сдвиг решения, однако график ошибки перемещается в противоположную сторону (при увеличении начального значения она уменьшается и наоборот), причем непропорционально. Следовательно, при увеличении начального значения X ошибка уменьшается, однако с увеличением t это изменение становится все менее и менее заметным. Это характерно как для метода с переменным шагом, так и для метода с постоянным шагом. Кроме того, для метода с постоянным шагом увеличение шага приводит к уменьшению точности решения. Что касается жесткой системы, то для нее характерно прямолинейное конечное решение. С уменьшением шага ошибка, посчитанная по методу Рунге, уменьшается. При увеличении начального X происходит параллельный сдвиг графика конечного решения. Если в случае нежесткой системы графики решений разных составляющих X находятся далеко друг от друга, то при решении жесткой системы получаются очень близкие графики.
Заключение
В ходе выполнения работы был изучен неявный метод Эйлера для решения линейных систем ОДУ. В ходе реализации проекта и проведения тестирования была проверена справедливость теоретических выкладок. Получены сведения о зависимости точности интегрирования от величины и способа выбора шага, а также от начальных значений переменных.
Используемая литература
1. Сарычева О.М. Численные методы. тАУ Новосибирск, 1995г. тАУ 65с.
2. Бахвалов Н.С. Численные методы. Ч1.- М: Наука, 1975г. тАУ 632с., илл.
3. Копченова Н.В., Марон И.А. Вычислительная математика в примерах и задачах . тАУ М: Наука, 1972г. - 368с
Вместе с этим смотрят:
РЖнварiантнi пiдпростори. Власнi вектори i власнi значення лiнiйного оператора
РЖнтерполювання функцiй
Автокорреляционная функция. Примеры расчётов
Актуальные проблемы квантовой механики
Алгебра и алгебраические системы