Исследование метода дифференцирования по параметру для решения нелинейных САУ (стр. 2 из 4)

3. Методы не требуют вычисления производных функций f_i (X,t), i=1,n, а только самой функции в нескольких точкахна шаге h_m .

Методом Рунге-Кутта 1-го порядка является явный метод Эйлера:

Х^/ =F(Х, t) ;

X^m+1 = X^m + h_m·F (X^m, t_m)

Ошибка аппроксимации ε^α ~ h²_т. Область абсолютной устойчивости – круг радиусом, равным 1 и центром в точке (0, -1) – см. рис. 1.3, кривая 1; область относительной устойчивости – вся правая полуплоскость.

Рассмотрим методы Рунге-Кутта 2-го и 4-го порядка, которые также используются довольно часто.

Алгоритм метода Рунге-Кутта 2-го порядка состоит в следующем:

X^m+1 = X^m + (h_m/2)·(K₁ + K₂), где

K₁= F(X^m, t_m), K₂= F(X^m+ h_m·K₁, t_m + h_m).

Ошибка аппроксимации ε^α=kh³_т. Область абсолютной устойчивости показана на рис. 1.3 (кривая 2). Область относительной устойчивости - вся правая полуплоскость.

Алгоритм метода Рунге-Кутта 4-го порядка:

X^m+1 = X^m + (h_m/6) · (K₁ + 2K₂ + 2K₃ + K₄), где

K₁= F(X^m, t_m), K₂= F(X^m + (h_m/2)·K₁, t_m + h_m/2);

K₃= F(X^m + (h_m/2) ·K₂, t_m + h_m/2);

K₄= F(X^m + h·K₃, t_m + h_m);

Ошибка аппроксимации ε^α=kh⁵_т. Область абсолютной устойчивости показана на рис. 1.3 (кривая 3). Область относительной устойчивости - вся правая полуплоскость.

1

2

3

Рис. 1.3

1.4 Метод Ньютона

Итерационная формула дискретного метода Ньютона имеет вид

X^m+1=X^m – J^-1(X^m) ·F(X^m) ,

где J (Х^m) = F^/_X/ _X=Xm– матрица Якоби.

Характеристики метода.

1. Сходимость.

Условие сходимости метода

ρ(G’_Х(Х))= ρ(I – (J^-1(Х) ·F(Х))^/_Х)<1.

Имеем ρ(I – J^-1(Х*) · F^/_Х(Х*))=0; это означает, что метод обладает локальной сходимостью, т.е. сходится только вблизи точного решения. Поэтому при реализации метода дифференцирования по параметру вначале нужно получить приближенное значение. Чем ближе к Х^*, тем быстрее сходится метод.

2. Выбор начального приближения Х⁰.

Поскольку метод сходится только вблизи точного решения, значит, начальное приближение также нужно выбирать вблизи Х*. Насколько близко, зависит от скорости изменения функции F(Х) вблизи Х*: чем выше скорость, тем меньше область, где соблюдается условие сходимости.

3. Скорость сходимости.

Она определяется соотношением

||E^m+1||=C_m||E^m||^{1+ Р},

0<р<1. При Х→Х* величина р→1. Это значит, что вблизи точного решения скорость сходимости близка к квадратичной.

4. Критерий окончания итераций.

Таким критерием может быть любое соотношение из описанных выше в обобщенном алгоритме решения нелинейных САУ.

1.5 Дискретный метод Ньютона

Трудности практического применения метода Ньютона состоят в следующем:

1.Необходимость определения матрицы J = F^/_X.

При этом существует два подхода:

- аналитический способ. Здесь метод Ньютона особенно эффективен. Однако точные формулы могут быть слишком громоздкими, что повышает вероятность ошибки. Кроме того функцииF(X)могут быть заданы таблично;

- конечно-разностная аппроксимация. При этом используется формула:

δf_i/δx_j = (f_i(x₁, …, x_j + ∆x_j, …, x_n) - f_i(x₁, …, x_j - ∆x_j, …, x_n)) / 2δx_j .

В этом случае имеем дискретный метод Ньютона, который уже не обладает квадратичной сходимостью. Скорость сходимости можно увеличить, уменьшая ∆x_j по мере приближения к X^*.

2. Вычисление матрицы J^-1 на каждом шаге требует значительных вычислительных затрат. Поэтому часто вместо этого решают систему линейных АУ, которая формируется следующим образом. Очевидно, что ∆X^m = X^m+1 - X^m. Тогда после алгебраических преобразований алгоритм X^m+1=X^m – J^-1(X^m) ·F(X^m) примет вид:

J(X^m) · ∆X^m = -F(X^m).

На каждом m-м шаге матрицыJ(X^m) и F(X^m) известны. Необходимо найти ∆X^m, как решение системы линейных АУ J(X^m) · ∆X^m = -F(X^m). Тогда

X^m+1 = X^m + ∆X^m.

Решение системы J(X^m) · ∆X^m = -F(X^m) – наиболее трудоемкий этап, который определяет вычислительную эффективность каждой итерации.

2. ОПИСАНИЕ ПРОГРАММНОГО ОБЕСПЕЧЕНИЯ

2.1 Общие сведения

Наименования файлов, из которых состоит пакет: main.m, rk1.m, rk2, rk4.m, funf.m, dmn.m, dif.m .

Головная программа – файл main.m, остальные – внешние функции данной программы.

Программное обеспечение: ОС класса Windows.

Язык, на котором написана программа: MatLab (Version 5.1).

2.2 Функциональное назначение

Программа разработана для решения нелинейных систем алгебраических уравнений методом дифференцирования по параметру. Решение проводится с использованием трех различных методов Рунге – Кутта первого, второго и четвертого порядка. Решение уточняется с помощью дискретного метода Ньютона.

2.3 Описание логической структуры

Файл описания системы - funf.m. С помощью файла dif.m формируется система ОДУ, затем система интегрируется с помощью явных методов Рунге-Кутта (rk1.m, rk2, rk4.m), и полученное решение уточняется дискретным методом Ньютона (dmn.m). В файле main.m можно изменять шаг интегрирования, начальное решение системы, допустимую ошибку при уточнении, начальный и конечный моменты времени…

Головная программа

В головной программе задаются начальные значения необходимых параметров, а также производится вызов основных функций, необходимых для реализации метода дифференцирования по параметру. В течение одного цикла работы программы вызываются последовательно функции rk1.m, rk2.m, rk4.m, вычисляющие приближенные значения по методам Рунге – Кутта 1-го, 2-го и 4-го порядков. После каждого из них производится вызов функции dmn.m, вычисляющей уточненное значение по дискретному методу Ньютона, строятся графики и выводятся значения полученных решений системы и ошибок интегрирования; значения уточненных решений системы и ошибок интегрирования, полученных в процессе итераций. В конце работы каждого цикла выводится число шагов, за которое было получено приближенное решение и число итераций, за которое было получено уточненное решение.

Текст программы (файл Main.m):

t0=0;

tfinal=1;

y0=[2 0];

h=0.1;

trace=1;

disp('Метод Рунге - Кутта 1го порядка');pause;

[tout,yout,x]=rk1('dif',t0,tfinal,y0,h,trace);

subplot(2,1,1);

plot(tout,yout(:,1));

grid;

title(sprintf('МетодРунге - Кутта 1го порядка'));

ylabel('y(1)'); xlabel('t');

subplot(2,1,2);

plot(tout,yout(:,2));

grid;

ylabel('y(2)'); xlabel('t');

pause;

x0=x';

dx=[0.01;0.01];

Ed=0.000001;

[x,dx,m] = dmn('funf',x0,dx,Ed);

subplot(111);

plot(m,x(:,1),m,x(:,2));

grid;

title(sprintf(‘Уточнение решения системы’));

ylabel('x(m)'); xlabel('m');

pause;

plot(m,dx(:,1),m,dx(:,2));

grid;

ylabel('dx(m)'); xlabel('m');

pause;

out=[m,x,dx]

disp('Метод Рунге - Кутта 2го порядка');pause;

[tout,yout,x]=rk2('dif',t0,tfinal,y0,h,trace);

subplot(2,1,1);

plot(tout,yout(:,1));

grid;

title(sprintf('МетодРунге-Кутта 2гопорядка'));

ylabel('y(1)'); xlabel('t');

subplot(2,1,2);

plot(tout,yout(:,2));

grid;

ylabel('y(2)'); xlabel('t');

pause;

x0=x';

dx=[0.01;0.01];

[x,dx,m] = dmn('funf',x0,dx,Ed);

subplot(111);

plot(m,x(:,1),m,x(:,2));

grid;

title(sprintf(‘Уточнение решения системы’));

ylabel('x(m)'); xlabel('m');

pause;

plot(m,dx(:,1),m,dx(:,2));

grid;

ylabel('dx(m)'); xlabel('m');

pause;

out=[m,x,dx]

disp('Метод Рунге - Кутта 4го порядка');pause;

[tout,yout,x]=rk4('dif',t0,tfinal,y0,h,trace);

subplot(2,1,1);

plot(tout,yout(:,1));

grid;

title(sprintf('МетодРунге-Кутта 4гопорядка'));

ylabel('y(1)'); xlabel('t');

subplot(2,1,2);

plot(tout,yout(:,2));

grid;

ylabel('y(2)'); xlabel('t');

pause;

x0=x';

dx=[0.01;0.01];

[x,dx,m] = dmn('funf',x0,dx,Ed);

subplot(111);

plot(m,x(:,1),m,x(:,2));

grid;

title(sprintf(‘Уточнение решения системы’));

ylabel('x(m)'); xlabel('m');

pause;

plot(m,dx(:,1),m,dx(:,2));

grid;

ylabel('dx(m)'); xlabel('m');

pause;

out=[m,x,dx]

Функции rk1, rk2, rk4

Данные функции являются программной реализацией методов Рунге-Кутта 1-го, 2-го и 4-го порядка. Здесь производится расчет приближенных значений при интегрировании системы с t = [0; 1]. Входные параметры функции – правые части системы, начальный и конечный моменты времени, начальное значение вектора Y, шаг и признак трассировки. Выходные параметры – вектор моментов времени tout, матрица yout={nx 2}, где n – число шагов метода( по строкам матрицы располагаются вектора Y^m, m = 0,n), вектор х, содержащий значения, полученные на последнем шаге работы функции и предназначенный для того, чтобы затем можно было передать значения данного вектора в подпрограмму, реализующую дискретный метод Ньютона.