Исследование метода дифференцирования по параметру для решения нелинейных САУ (стр. 3 из 4)
Тексты программ (файл 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
Файл funf.m содержит исходную систему нелинейных САУ (описание ее правых частей). Пользователь может ввести в файл любую нелинейную систему ОДУ.
Текст файла funf.m:
function [F] = funf(x)
F=[x(1)*x(1)+x(2)*x(2)-4; x(1)*x(2)-1];
Файл dif.m
Файл формируется для того, чтобы создать систему дифференциальных уравнений, интегрируемых впоследствии по методам Рунге – Кутта. Он содержит матрицы начальных значений системы и значения матрицы Якоби. Функция вычисляет произведение обратной матрицы Якоби и начального решения системы.
Текст программы(файл dif.m):
function yp= dif(t,y)
J=[2*y(1) 2*y(2); y(2) y(1)];
J=-inv(J);
a=[0; -1];
b=J*a;
yp=b;
Файл Dmn.m
Файл содержит подпрограмму, вычисляющую уточненное решение системы по дискретному методу. Он выводит количество пройденных итераций на экран. Выходными данными является вектор mout и матрицы xout, dxout.
Текст программы(файл dif.m):
function[xout,dxout,mout]=dmn(funf,x,dx,edop);
xout=x';
dxout=dx';
x1=x;
m=0; it=0;
mout=m;
nv=[1;1];
n=size(x);
while(max(nv)>edop)
f=feval(funf,x);
nf=norm(f);
for j=1:n,
x1(j)=x(j)+dx(j);
f1=feval(funf,x1);
x1(j)=x(j)-dx(j);
f2=feval(funf,x1);
J(:,j)=(f1-f2)/(2*dx(j));
x1(j)=x(j);
end;
dx=-J\f;
ndx=norm(dx);
nv=[nf;ndx];
x=x+dx ;
m=m+1; it=it+1;
xout=[xout;x1'];
dxout=[dxout;dx'];
mout=[mout;m];
end;
disp('Количество итераций равно'); disp(it);
2.4 Используемые технические средства
ЭВМ, совместимые с IBM PC. Процессор не ниже Pentium1. ОП не меньше 32 Мб, мышь, стандартная клавиатура, видеокарта.
2.5 Вызов и загрузка
Для запуска программы из среды MatLab необходимо вызвать головную программу, набрав в командной строке ее имя main. При этом необходимо перейти в тот каталог, в котором находится данный пакет программ. Затем производятся все предусмотренные программой операции. После вывода на экран каждого графика или сообщения система ожидает нажатия любой клавиши.
Программа со всеми сопутствующими файлами, занимает 4,3 Кб.
2.6 Входные данные
Правые части системы в файле funf.m, матрица Якоби в файле dif.m, система дифференциальных уравнений, составленная по исходной системе нелинейных САУ (в файле dif.m), начальное приближение, значение шага, интервал времени для интегрирования, допустимая ошибка.
2.7 Выходные данные
На выходе программы на экран выводятся приближенные решения, полученные по методам Рунге-Кутта, значения шага и времени на каждом шаге интегрирования, количество шагов интегрирования, графики приближенных решений по методам Рунге-Кутта, количество итераций, необходимых для нахождения уточненного решения дискретным методом Ньютона, уточненное решение при каждой итерации, графики уточненных значений по методу Ньютона и ошибки для каждой составляющей решения.
3. ОПИСАНИЕ ТЕСТОВЫХ ЗАДАЧ
Решение системы нелинейных САУ.
Для интегрирования возьмем систему:
x2+y2-4=0xy – 1=0
Тогда при запуске программы на экране появляются следующие сообщения:
Метод Рунге - Кутта 1го порядка
t =
0
h =
0.1000
y =
2 0
…
t =
1
h =
1.1102e-016
y =
1.9398 0.5139
Количество шагов =
11
График решения системы ДУ:
 |
Количество итераций равно
3
Графики значений системы и ошибок при каждой итерации:
 |
 |
out =
0 1.9398 0.5139 0.0100 0.0100
1.0000 1.9398 0.5139 -0.0079 0.0037
2.0000 1.9319 0.5176 0.0000 0.0000
3.0000 1.9319 0.5176 0.0000 0.0000
Метод Рунге - Кутта 2го порядка
t =
0
h =
0.1000
y =
2 0
…
t =
1
h =
1.1102e-016
y =
1.9319 0.5176
Kоличество шагов =
11
График решения системы ДУ:
 |
Количество итераций равно
2
Графики значений системы и ошибок при каждой итерации:
 |
 |
out =
0 1.9319 0.5176 0.0100 0.0100
1.0000 1.9319 0.5176 0.0000 0.0000
2.0000 1.9319 0.5176 0.0000 0.0000
Метод Рунге - Кутта 4го порядка
t =
0
h =
0.1000
y =
2 0
…
t =
1
h =
1.1102e-016
y =
1.9291 0.5190
Kоличество шагов =
11
График решения системы ДУ:
 |
Количество итераций равно
3
Графики значений системы и ошибок при каждой итерации:
 |
 |
out =
0 1.9291 0.5190 0.0100 0.0100
1.0000 1.9291 0.5190 0.0027 -0.0014
2.0000 1.9319 0.5176 0.0000 0.0000
3.0000 1.9319 0.5176 0.0000 0.0000
Проверим теперь влияние задаваемого шага интегрирования на точность получаемого решения: зададим h = 0.5 вместо 0.1. Тогда получим:
Метод Рунге - Кутта 1го порядка
t =
0
h =
0.5000
y =
2 0
…
t =
1
h =
0.5000
y =
1.9683 0.5040
Количество шагов =
2
Количество итераций равно
4
out =
0 1.9683 0.5040 0.0100 0.0100
1.0000 1.9683 0.5040 -0.0359 0.0133
2.0000 1.9323 0.5173 -0.0005 0.0004
3.0000 1.9319 0.5176 0.0000 0.0000
4.0000 1.9319 0.5176 0.0000 0.0000
Метод Рунге - Кутта 2го порядка
t =
0
h =
0.5000
y =
2 0
t =
1
h =
0.5000
y =
1.9321 0.5175
Kоличество шагов =
2
Количество итераций равно
2
out =
0 1.9321 0.5175 0.0100 0.0100
1.0000 1.9321 0.5175 -0.0002 0.0002
2.0000 1.9319 0.5176 0.0000 0.0000
Метод Рунге - Кутта 4го порядка
t =
0
h =
0.5000
y =
2 0
t =
1
h =
0.5000
y =
1.9183 0.5243
Kоличество шагов =
2
Количество итераций равно
3
out =
0 1.9183 0.5243 0.0100 0.0100
1.0000 1.9183 0.5243 0.0137 -0.0068
2.0000 1.9319 0.5176 -0.0001 0.0001
3.0000 1.9319 0.5176 0.0000 0.0000
Видим, что при увеличении h снизилась точность получаемого приближенного решения, уменьшилось количество шагов по методу Рунге – Кутта (их стало не 11, а 3), и, вследствие этого, увеличилось количество итераций по дискретному методу Ньютона.
Проверим влияние задаваемой допустимой ошибки для дискретного метода Ньютона: зададим edop = 0.001 вместо edop = 0.00001. Получаем:
Метод Рунге - Кутта 1го порядка
t =
0
h =
0.1000
y =
2 0
…
t =
1
h =
1.1102e-016
y =
1.9398 0.5139
Количество шагов =
11
Количество итераций равно
2
out =
0 1.9398 0.5139 0.0100 0.0100
1.0000 1.9398 0.5139 -0.0079 0.0037
2.0000 1.9319 0.5176 0.0000 0.0000
Метод Рунге - Кутта 2го порядка
t =
0
h =
0.1000
y =
2 0
…
t =
1
h =
1.1102e-016
y =
1.9319 0.5176
Kоличество шагов =
11
Количество итераций равно
1
out =
0 1.9319 0.5176 0.0100 0.0100
1.0000 1.9319 0.5176 0.0000 0.0000
Метод Рунге - Кутта 4го порядка
t =
0
h =
0.1000
y =
2 0
…
t =
1
h =
1.1102e-016
y =
1.9291 0.5190
Kоличество шагов =
11
Количество итераций равно
2
out =
0 1.9291 0.5190 0.0100 0.0100
1.0000 1.9291 0.5190 0.0027 -0.0014
2.0000 1.9319 0.5176 0.0000 0.0000
Видим, что при увеличении допустимой ошибки для дискретного метода Ньютона уменьшается число итераций, так как уже при второй, и даже первой, итерации достигается заданная точность решения.
Решим эту же систему при другом начальном приближении Х0 = (3 0).
Метод Рунге - Кутта 1го порядка
t =
0
h =
0.1000
y =
3 0
…
t =
1
h =
1.1102e-016
y =
3.7401 0.2728
Количество шагов =
11
Количество итераций равно
6
out =
0 3.7401 0.2728 0.0100 0.0100
1.0000 3.7401 0.2728 -1.3520 0.0932
2.0000 2.3880 0.3660 -0.3996 0.0984
3.0000 1.9884 0.4644 -0.0539 0.0484
4.0000 1.9345 0.5128 -0.0026 0.0047
5.0000 1.9319 0.5176 0.0000 0.0001
6.0000 1.9319 0.5176 0.0000 0.0000
Метод Рунге - Кутта 2го порядка
t =
0
h =
0.1000
y =
3 0
…
t =
1
h =
1.1102e-016
y =
3.7321 0.2680
Kоличество шагов =
11
Количество итераций равно