%% Merev rendszerek
% 1. Oldjuk meg az alábbi reakciókinetikai kezdetiértékfeladatot a [0,50] időintervallumon! A kezdeti érékek rendre: 1,1 és 0. Próbáljuk ki a megoldást a Matlab ode45, ode23s függvényeivel ill. implicit Euler-módszerrel!
% c_1^\prime=-0.013c_1-1000c_1c_3\\ c_2^\prime=-2500c_2c_3\\ c_3^\prime=-0.013c_1-1000c_1c_3-2500c_2c_3

f=@(x,y)[-0.013*y(1)-1000*y(1)*y(3);-2500*y(2)*y(3);-0.013*y(1)-1000*y(1)*y(3)-2500*y(2)*y(3)];

% ode45(f,[0,50],[1,1,0])

% ode23s(f,[0,50],[1,1,0])

N=50;
h=50/N;

x=zeros(1,N+1); 
x(1)=0;
y=zeros(3,N+1);
y(:,1)=[1;1;0];

J=@(y)[-0.013-1000*y(3)-1/h,0,-1000*y(1);0,-2500*y(3)-1/h,-2500*y(2);-0.013-1000*y(3),-2500*y(3),-1000*y(1)-2500*y(2)-1/h];

for i=1:N
    x(i+1)=x(i)+h;
    z=y(:,i);
    for k=1:10
    z=z-J(z)\(y(:,i)+h*f(x(i+1),z)-z);
    end
    y(:,i+1)=z;
end

plot(x,y,'bo-')

%
%% Többlépéses módszerek
% 2. Írjunk programot, ami meghatározza a többlépéses módszerek esetén a konzisztenciavizsgálathoz szükséges d_i együtthatókat!

format rat
lepesszam=2;
 a=[1,-4/3,1/3], b=[2/3,0,0]
% a=[1,-1,0], b=[0,3/2,-1/2]

i=lepesszam:-1:0;
d0=sum(a); disp('d0:'); disp(d0*(abs(d0)>10^-10))
d1=sum(i.*a-b); disp('d1:'); disp(d1*(abs(d1)>10^-10))
d2=sum((i.^2).*a/2-i.*b); disp('d2:'); disp(d2*(abs(d2)>10^-10))
d3=sum((i.^3).*a/6-(i.^2).*b/2); disp('d3:'); disp(d3*(abs(d3)>10^-10))
d4=sum((i.^4).*a/24-(i.^3).*b/6); disp('d4:'); disp(d4*(abs(d4)>10^-10)) 
d5=sum((i.^5).*a/120-(i.^4).*b/24); disp('d5:'); disp(d5*(abs(d5)>10^-10)) 
d6=sum((i.^6).*a/720-(i.^5).*b/120); disp('d6:'); disp(d6*(abs(d6)>10^-10))


format short

%
% 3. Készítsük el a kétlépéses BDF formulát (BDF2)! Vizsgáljuk meg stabilitás, konzisztencia és konvergencia szempontjából! Ellenőrizzük le a konvergenciarendet az y^\prime=x^2(1+y), y(0)=3 kezdetiértékfeladat megoldásával a [0,1] intervallumon! Az első lépést számítsuk a RK2 módszerrel!

clear all
ypontos=@(x)4*exp(x.^3/3)-1;
f=@(x,y)x.^2.*(1+y);
x0=0, y0=3,
xmax=1, N=80, xx=linspace(x0,xmax,500);

h=(xmax-x0)/N;

x=zeros(1,N+1); y=x;
x(1)=x0;
y(1)=y0;
x(2)=x(1)+h;
y(2)=y(1)+h*f(x(1)+h/2,y(1)+h/2*f(x(1),y(1))); % RK2

for i=2:N
    x(i+1)=x(i)+h;
    z=y(i);
    hiba=1;
    while hiba>10^-6      
    zuj=2*h/3*f(x(i+1),z)+4/3*y(i)-1/3*y(i-1);
    hiba=abs(z-zuj);
    z=zuj;
    end
    y(i+1)=z;
end

hiba=y(end)-ypontos(xmax)

plot(x,y,'bo-',xx,ypontos(xx),'r-')


%
% 4. Készítsük el a kétlépéses Adams-Bashforth-formulát, és végezzünk hasonló vizsgálatot, mint az előző feladatban!

clear all
ypontos=@(x)4*exp(x.^3/3)-1;
f=@(x,y)x.^2.*(1+y);
x0=0, y0=3,
xmax=1, N=160, xx=linspace(x0,xmax,500);

h=(xmax-x0)/N;

x=zeros(1,N+1); y=x;
x(1)=x0;
y(1)=y0;
x(2)=x(1)+h;
y(2)=y(1)+h*f(x(1)+h/2,y(1)+h/2*f(x(1),y(1))); % RK2

for i=2:N
    x(i+1)=x(i)+h;
    y(i+1)=y(i)+h*(3/2*f(x(i),y(i))-1/2*f(x(i-1),y(i-1)));
end

hiba=y(end)-ypontos(xmax)

plot(x,y,'bo-',xx,ypontos(xx),'r-')

%
%% Peremértékfeladatok
% 5. Oldjuk meg a belövéses módszer segítségével az y^{\prime\prime}=y+4e^x, y(0)=1, y(1/2)=2e^{1/2} peremértékfeladatot! A kezdetiérték-feladat megoldásához használjuk a Heun-módszert ill. az intervallumfelezési módszert! (Pontos megoldás: y(x)=(2x+1)e^x)

clear all
ypontos=@(x)(2*x+1).*exp(x);
a=0; A=1; b=1/2; xx=linspace(a,b,500);
N=100; h=(b-a)/N;
D=1; % kezdő derivált
B=2*exp(1/2); % a megcélzott jobb oldali fvérték

x(1)=a;

f=@(x,y)[y(2);y(1)+4*exp(x)];


y=[A;D];    

for i=1:N
    x(i+1)=x(i)+h;
    y(:,i+1)=y(:,i)+h*(1/2*f(x(i),y(:,i))+1/2*f(x(i)+h,y(:,i)+h*f(x(i),y(:,i))));
end

plot(x,y(1,:),'bo-',xx,ypontos(xx),'r-')

%
% 6. Oldjuk meg az előző feladatban szereplő peremértékfeladatot a véges differencia módszerrel!

clear all
ypontos=@(x)(2*x+1).*exp(x);
a=0; A=1; b=1/2; xx=linspace(a,b,500);
N=4; h=(b-a)/(N+1);
x=[a+h:h:b-h]';

B=2*exp(1/2); % a megcélzott jobb oldali fvérték

y=(gallery('tridiag',N,1,-2,1)/h^2-speye(N))\(4*exp(x)-[A/h^2,zeros(1,N-2),B/h^2]');
plot([a,x',b],[A,y',B],'bo-',xx,ypontos(xx),'r-')

hiba=(norm(ypontos(x)-y,inf))

%
% 7. Oldjuk meg az alábbi lehajló rúd egyenletet az előző programok segítségével!
% y^{\prime\prime}=4x^2\cos y, y(0)=0, y(1)=0