%% -u"=12x.^2-6x megoldása 

% [0,1] szakaszon, homogén Dirichlet-perem, 
% szakaszonként lineáris elemek, hibaszámítás

close all, clear all
% jobb oldal
f=@(x)(12*x.^2-6*x);

% pontos megoldás

 upontos=@(x)-x.^4+x.^3;
 upontosder=@(x)-4*x.^3+3*x.^2;

% adatok
n=3;
h=1/(n+1);

% Merevségi mátrix

A=gallery('tridiag',n,-1,2,-1)/h;

% Terhelési vektor
feloszt=h:h:1-h;
b=h*f(feloszt'); %h/3*(f(feloszt-h/2)+f(feloszt)+f(feloszt+h/2))';

% Egyenletrendszer megoldása
y=A\b;

y=[0,y',0]';

%% Ábrázolás 
% Ha van pontos megoldás, akkor azt is kirajzoljuk. Lineáris elemeknél
% egyenes szakaszokat, kvadratikusnál parabolaíveket rajzolunk.

if length(y)== n+2
    dim=1;
elseif length(y)== 2*n+3
    dim=2;
end;

hold on
for i=1:n+1
    xkord=[(i-1)*h:h/dim:i*h]';
if dim==1
    ykord=y(i:i+1);
else
    ykord=y(2*i-1:2*i+1);
end   
egy=polyfit(xkord,ykord,dim);
xx=[(i-1)*h:h/30:i*h];
zz=polyval(egy,xx);
plot(xx,zz,'b-')
end
plot(0:h/dim:1,y','bo')
if exist('upontos','var')==1
    pg=ezplot(upontos,[0,1]);
    set(pg,'color','red')
end;
hold off

%% Hibaszámítás
% Ha van pontos megoldás, akkor kiszámítjuk az L2 és H1 hibákat.

if exist('upontos','var')==1 && exist('upontosder','var')==1

L2negyzet=0;
H1seminegyzet=0;

for i=1:n+1
xkord=[(i-1)*h:h/dim:i*h]';
if dim==1
    ykord=y(i:i+1);
else
    ykord=y(2*i-1:2*i+1);
end   
egy=polyfit(xkord,ykord,dim);
elteresnegyzet=@(x)(polyval(egy,x)-upontos(x)).^2;
L2negyzet=L2negyzet+quad(elteresnegyzet,(i-1)*h,i*h);
derivelteresnegyzet=@(x)(polyval(polyder(egy),x)-upontosder(x)).^2;
H1seminegyzet=H1seminegyzet+quad(derivelteresnegyzet,(i-1)*h,i*h);
end;
    
L2=sqrt(L2negyzet)
H1=sqrt(L2negyzet+H1seminegyzet)

end