code éléments finis P1

funcprot(0);
lines(0);
clear;

function y=a(x)
    if (x < 1/2)  then y= x^{1-alpha};
    else y=1;
    end
endfunction

function y=f(x)
    c1= -4*(1/2)^{alpha} -2*alpha;
    if (x < 1/2)  then y=0;
    else y=-2*c1;
    end
endfunction

function I=integrale(a,b,r)
    q=sqrt(3/5);
    [H1,H2,H3]=(5/9,8/9,5/9);
    [e1, e2, e3]=(a+(1-q)*(b-a)/2, (a+b)/2, a+(1+q)*(b-a)/2);
    
    I=((b-a)/2)*(H1*r(e1)+H2*r(e2)+H3*r(e3));
    
endfunction

function b=CalculSMD(N, f)
    h=1/(N+1);
    for i=1:N 
        b(i)=f(i*h)*h;
    end
endfunction

function A=CalculMatD(N,a)
    A=zeros(N,N);
    h=1/(N+1);
    for i=1:N-1
        A(i,i)=integrale((i-1)*h,i*h,a)+integrale(i*h,(i+1)*h,a);
        A(i,i+1)=-integrale(i*h,(i+1)*h,a);
        A(i+1,i)=-integrale(i*h,(i+1)*h,a);
    end;
    i=N;
    A(N,N)=integrale((N-1)*h,N*h,a)+integrale(N*h,(N+1)*h,a);
    A=A/h^2;
endfunction

function y=uapproche(N,a,f)
    y=[0;CalculMatD(N,a)\CalculSMD(N, f);0];
endfunction

function y=uderive(N,a,f)
    h=1/(N+1);
    b=uapproche(N,a,f);
    y(1:N+1)=(b((1:N+1)+1)-b(1:N+1))/h;
    y(N+2)=y(N+1);
endfunction

function y=u(x)
    [c1,c2,c3]=(4*(1/2)^{alpha} -2*alpha, 4*(1/2)^{alpha}+ 3*alpha, alpha);
    if x<1/2 then y=x^{alpha};
    else y=c1*x^2 + c2*x + c3;
    end
endfunction

function y=ud(x)
    [c1,c2]= (4*(1/2)^{alpha} -2*alpha, 4*(1/2)^{alpha}+ 3*alpha);
    if x<1/2 then y=alpha*x^(alpha-1);
    else y = 2*c1*x + c2;
    end
endfunction

function y=u(x,i,N,d)
    h=1/(N+1);
    y=d(i+1) + (x-i*h)*((d(i+2)-d(i+1))/h);
endfunction

function y=Erreur(N,a,f)
    d=uapproche(N,a,f);
    h=1/(N+1);
    [H1,H2,H3]=(5/9, 8/9, 5/8);
    q=sqrt(3/5);
    y=0;
    for i=0:N
        [a,b]=(i*h, (i+1)*h);
        [e1, e2, e3]=(a+(1-q)*(b-a)/2, (a+b)/2, a+(1+q)*(b-a)/2);
        y = y + H1*(u(e1,i,N,d)-u(e1))^2;
        y = y + H2*(u(e2,i,N,d)-u(e1))^2;
        y = y + H3*(u(e3,i,N,d)-u(e3))^2;
        //disp(y);
    end
    y=((b-a)/2)*y;
endfunction



//Pour la questioin 5 de l'exercice 6 
function AfficheConditionnement()
    alpha=0.7;
    for N=1:100
        h(N)=1/(N+1);
        A=CalculMatD(N,a);
        condit(N)=cond(A);
    end
        plot2d(h,condit,13,logflag="ll")
        
    alpha=1;
    for N=1:100
        A=CalculMatD(N,a);
        condit(N)=cond(A);
    end
        plot2d(h,condit,5,logflag="ll")
    
        plot2d(h,h^{-1},26,logflag="ll")
        plot2d(h,h^{-2},2,logflag="ll")
        
    xtitle('Conditionnement de la matrice A_h')
    legends(['1/h^2','alpha=0.7','alpha=1','1/h'],[2,13,5,26],1)
    xs2pdf (0,'ex6_5.pdf');
endfunction

// Pour la question 1 de l'exercice 7
function AfficheAppExacte()
    N=2^6-1;
    h=1/(N+1);
    r(1:N+2)=(0:N+1)*h;
    alpha=0.7;
        plot2d(r,uapproche(N,a,f),1);
    alpha=1;
        plot2d(r,uapproche(N,a,f),22);
        
    xtitle('Solution approché (sa) en vis à vis de la solution exacte (se)')
    legends(['(sa) alpha=0.7','(se) alpha=0.7','(sa) alpha=1','(se) alpha=1'],[1,2,22,5],1)
    xs2pdf (0,'ex7_1.pdf');
endfunction

// Pour la question 2 de l'exercice 7
function AfficheDerAppExacte()
    N=2^6-1;
    h=1/(N+1);
    r(1:N+2)=(0:N+1)*h;
    alpha=0.7;
        plot2d(r,uderive(N,a,f),1);
    alpha=1;
        plot2d(r,uderive(N,a,f),22);
        
    xtitle('Dérivée de la solution approché (sa) en vis à vis de la solution exacte (se)')
    legends(['(sa) alpha=0.7','(se) alpha=0.7','(sa) alpha=1','(se) alpha=1'],[1,2,22,5],1)
    xs2pdf (0,'ex7_2.pdf');
endfunction

//Pour la question 3 de l'exercice 7
function AfficheErreur()
    clf();
    
    alpha=1;
    x(1:8)= 2^(1:8)-1;
    for j=1:8
        y(j)=Erreur(2^j-1,a,f);
        H(j)= 1/(x(j)+1);
    end
        plot2d(H,y,5,logflag="ll");
        
    alpha=0.7;
    for j=1:8
        y(j)=Erreur(2^j-1,a,f);
    end
        plot2d(H,y,2,logflag="ll");
        plot2d(H,H^2,26,logflag="ll");
        
    xtitle('Erreur en norme L2');
    legends(['alpha=0.7','alpha=1','h^2'],[2,5,26],1);
    xs2pdf (0,'ex7_3.pdf');
endfunction