zadanie2

tspan=[0 5];
y0=[1/sqrt(2) 0];
N=50;
[t,Y]=RK3another(@ftest, tspan, y0, N);
y1=Y(:,1);%Приближенное решение
y2=Y(:,2);%Приближенное решение
u1=cos(t)./(1+ exp(2.*t)).^(1/2);
u2=sin(t)./(1+ exp(2.*t)).^(1/2);

% Сравнение приближенного решения с точным%%
fig1 = figure;
plot(t,u1)
hold on
plot(t,y1,'o');
xlabel('t')
ylabel('y_1')
legend('точное','прибл.')
title('N=50')

fig2 = figure;
plot(t,u2)
hold on
plot(t,y2,'o');
xlabel('t')
ylabel('y_2')
legend('точное','прибл.')
title('N=50')

% tspan = [0 5];
% N = [30 50 80 100 150 200 300 400 500 1000];
% y0=[1/sqrt(2) 0];
% e = zeros(1,numel(N));
% hs = zeros(1,numel(N));
% for i=1:numel(N)
%     [t,Y,h]=RK3another(@ftest, tspan, y0, N(i));
%     hs(i)=h;
%     y1=Y(:,1);
%     y2=Y(:,2);
%     u1=cos(t)./(1+ exp(2.*t)).^(1/2);
%     u2=sin(t)./(1+ exp(2.*t)).^(1/2);
%     e1=norm(y1-u1,inf);
%     e2=norm(y2-u2,inf);
%     e(i)=max(e1,e2);
% end

% alpha=zeros(1,numel(N)-1);
% for i=1:numel(N)-1
%     alpha(i)=log10(e(i+1)/e(i))/log10(hs(i+1)/hs(i));
% end

% f2 = figure;
% loglog(hs,e)
% xlabel('h')
% ylabel('e')


% e_alpha = zeros(1, numel(alpha));
% h_alpha = zeros(1, numel(alpha));
% for i=1:numel(alpha)
%     h_alpha(i) = hs(i);
%     e_alpha(i) = e(i)/h_alpha(i)^(alpha(i));
% end
% f3 = figure;
% loglog(h_alpha, e_alpha)
% xlabel('h')
% ylabel('e/h^α')

% C_h = zeros(1, numel(hs));
% for i=1:numel(hs)
%     C_h(i)=e(i)/hs(i)^3;
% end
%
% f4 = figure;
% loglog(hs,C_h);
% xlabel('h')
% ylabel('C(h)')


%%%%%
function [T,Y,h]=RK3another(f,tspan,y0,N)
    n=numel(y0);
    T=zeros(N,1);
    Y=zeros(N,n);
    h=(tspan(2)-tspan(1))/(N-1);
    Y(1,:)=(y0(:))';
    T(1)=tspan(1);
    for i=1:N-1
        K1=(f(T(i),Y(i,:)))';
        K2=(f(T(i)+h/3,Y(i,:)+(h/3)*K1))';
        K3=(f(T(i)+2*h/3,Y(i,:)+(2*h/3)*K2 ))';
        Y(i+1,:)=Y(i,:)+(h/4)*(K1+3*K3);
        T(i+1)=T(i)+h;
    end
end



%%%%%%
function dydt=ftest(t,y)
    dydt=zeros(2,1);
    dydt(1)=-y(2)+y(1)*(y(1)^2+y(2)^2-1);
    dydt(2)=y(1)+y(2)*(y(1)^2+y(2)^2-1);