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

 

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

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 ) – наиболее трудоемкий этап, который определяет вычислительную эффективность каждой итерации.

 

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

 

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

 

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

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

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

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

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

 

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

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

 

Файл описания системы - 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={n x 2}, где n – число шагов метода( по строкам матрицы располагаются вектора Y^m, m = 0,n), вектор х, содержащий значения, полученные на последнем шаге работы функции и предназначенный для того, чтобы затем можно было передать значения данного вектора в подпрограмму, реализующую дискретный метод Ньютона.

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

 

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

t=t0;y=y0;

tout=t;

yout=y;

b=0;

if trace

clc, t, h, y

end

while (t<tfinal)

if t+h>tfinal

h=tfinal-t;

end

if trace

clc, t, h, y

end

k1=feval(dif, t, y);

q=h*k1;

y=y+q';

t=t+h;

b=b+1;

tout=[tout; t];

yout=[yout; y];

if trace

home,t, h, y

end;

end; x=y;

disp('Количество шагов = '); disp(b);

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

t=t0;y=y0;

tout=t;

yout=y;

b=0;

if trace

clc, t, h, y

end

while (t<tfinal)

if t+h>tfinal

h=tfinal-t;

end

if trace

clc, t, h, y

end

k1=feval(dif, t, y);

z=h*k1;

k2 = feval(dif, t+h, y+z');

q=h/2*(k1+k2);

y=y+q';

t=t+h;

b=b+1;

tout=[tout; t];

yout=[yout; y];

if trace

home,t, h, y

end;

end; x=y;

disp('Kоличество шагов = '); disp(b);

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

t=t0;y=y0;

tout=t;

yout=y;

b=0;

if trace

clc, t, h, y

end

while (t<tfinal)

if t+h>tfinal

h=tfinal-t;

end

if trace

clc, t, h, y

end

k1=feval(dif, t, y);

z=h*k1;

k2 = feval(dif, t+h, y+z');

z=h*k2/2;

k3 = feval(dif, t+h/2, y+z');

z=h*k3;

k4 = feval(dif, t+h, y+z');

q=h*(k1+2*k2+2*k3+k4)/6;

y=y+q';

t=t+h;

b=b+1;

tout=[tout; t];

yout=[yout; y];

if trace

home,t, h, y

end;

end; x=y;

disp('Kоличество шагов = '); disp(b);

 

Файл funf.m