Blatt4_13

% B4A13 geschrieben für Octave, Adams Moulton, Zeitvariable x
%Exakte Lösung ist y(x)= e^(x-x^2) (trivial), damit erhalten wir y(1,5)=0.47237


function Blatt3_13 %B4A13
    Moulton()
    exakt()
    BDF()
endfunction

%--------------------------Moulton-----------------------------------
function Moulton

f = inline('(1-2*x)*y','x','y');
h=0.1;
x=  1 : h : 1.5  ;
y = x;
y(1)=1;
% Startwerte suchen (Runge-Kutta)

    k1 = f(x(1), y(1));
    k2 = f(x(1) + h/2, y(1) + h/2*k1);

    y(2) = y(1) + h*k2;

% Iterieren:
for i = 3:length(y)
    % Adams-Moulton
  % Löse y(i) - (y(i-1) + h/12*(5*f(x(i),y(i)) + 8*f(x(i-1),y(i-1))- f(x(i-2),y(i-2))))=0;
  % Löse y(i) + h/12*5*f(x(i),y(i)) - (y(i-1)  + h/12 (8*f(x(i-1),y(i-1))- f(x(i-2),y(i-2))))=0;
  %Löse nach y(i) mit Sekantenverfahren, Startwerte aus Runge Kutta
    %Vereinfachung des Terms
  c= h/12*(8*f(x(i-1),y(i-1))- f(x(i-2),y(i-2)))+y(i-1); %Konstanter Wert

  y(i)=approx(y(i-2),y(i-1),@(t) t - h/12 * 5* f(x(i),t)-c);

endfor
y
disp('Damit ist der Wert an der Stelle 1.5:')
y(length(y))

figure('Name','Moulton-Plot','NumberTitle','off')
plot(x,y)
title ("Ergebnis mit Moulton");


endfunction

%------------------------------Exakte Lösung----------------------
function exakt
h=0.1;
x=1:h:1.5;
y=arrayfun(@(t) e^(t-t^2), x);
figure('Name','exakter Plot','NumberTitle','off')
disp('Exakte Lösung ist y(x)= e^(x-x^2) (trivial), damit erhalten wir y(1,5)=0.47237')
plot(x,y)
title("Exakte Lösung");

endfunction

%-------------------------BDF-Formel---------------------------------------
function BDF


f = inline('(1-2*x)*y','x','y');
h=0.1;
x=  1 : h : 1.5  ;
y = x;
y(1) = 1;
% Startwerte suchen (Runge-Kutta)

    k1 = f(x(1), y(1));
    k2 = f(x(1) + h/2, y(1) + h/2*k1);

    y(2) = y(1) + h*k2;

% Löse auf: 0= 2/3*h*f(x(i),y(i)) - y(i) -(-4/3 y(i-1)+1/3*y(i-2)). Setze:
for i=3:length(y)
    c= (-4/3*y(i-1)+1/3*y(i-2))
    y(i)=approx(y(i-2),y(i-1),@(t) 2/3*h*f(x(i),t)-t-c)




endfor
y
disp('Damit ist der Wert an der Stelle 1.5:')
y(length(y))

figure('Name','BDF-Plot','NumberTitle','off')
plot(x,y)
title('BDF Lösung');

endfunction

%-----------------------------------
function r=approx(in1,in2,g)
  %test=1
  while abs(g(in2))>0.00001
    in1=in2-(in2-in1)/(g(in2)-g(in1)) * g(in2);
    in2=in1-(in1-in2)/(g(in1)-g(in2)) * g(in1);
    %test=test+1
  endwhile
  r=in2;
endfunction

%erstellt mit Tobias Post/Philipp Heering