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={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

Файл 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=-Jf;

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);