Solving a problem with the Energy method with Hermite cubic function

I tried to solve the problem written below with the energy method using the cubic Hermite function, but only the final value of the approximation is incorrect.

How can I fix this error?

Equation:

-((1+x)y')'+(2+x)/(1+x)y=1

BC:

y(0)=0 ESSENTIAL

y(2)+3y'(2)=1

Exact solution (to be able to compare results): y=x./(1+x)

For testing

>>[x,u]= EnergyHermiteCubicFunctionEssential(58,63);

>> y=x./(1+x);

>>plot(x,u,'r*',x,y,'k-')

My code is:

[x,u]= EnergyHermiteCubicFunctionEssential(58,63);
y=x./(1+x);
plot(x,u,'r*',x,y,'k-')

function [x,u]= EnergyHermiteCubicFunctionEssential(m,nv)
%Energy method with Hermite cubic function

%Input:
%m: number of element of basis
%nv: number of visualization point

%Output:
% x: Vector containing the points at which the solution is approximated
% u: Approximated solution

%Equation:
%-((1+x)y')'+(2+x)/(1+x)y=1
%BC:
%y(0)=0 ESSENTIAL
%y(2)+3y'(2)=1 
%Exact solution: y=x./(1+x);


h=2/m;
xc=0:h:2;
%we have 2m+1 elements on the basis
A=zeros(2*m+1); %matrix squared
b=zeros(2*m+1,1); %coloumn vector

for l=1:2*m+1
    for j=1:2*m+1
        if abs(l-j)<=3 %indices have distance minore o uguale a 3
            [al,bl]=support(l,xc);
            [aj,bj]=support(j,xc);
        lb=max(al,aj); %massimo estremo sinistro
        rb=min(bl,bj); %minimo estremo destro
        
        if abs(lb-rb)>10^(-10) %tol
% 
% i termini di bordo li hO già inseriti (vedi linea 35)
                A(l,j)=integral (@(x) fA(x,l,j,xc),lb,rb);

        end
        
      end
    end
    
    [lb,rb]=support(l,xc);
    if abs(rb-lb)>10^(-10) %is the toll
    b(l)=integral(@(x) fB(x,l,xc),lb,rb);
    end
end

A(2*m,2*m)=A(2*m,2*m)+1; %1 is the boundery therm

delta=A\b; %delta of approximated solution

h=2/nv; %insert the visualization points
x=0:h:2;

u=zeros(1,nv+1); %row vector %nv+1

%for l=0:m-1 %there are m elements
% devi usare gli indici formali
    for l=1:2*m+1 %there are m elements
        u=u+delta(l)*fu(x,l,xc);
    end
u=u+(1/5.*x); %u(x)=z(x)+l(x); l(x)=1/5 * x;
end

function y=fH(x,l,xc)
m=length(xc)-1;
y=zeros(size(x));
    for j=1:length(x)
        if l==0 %first special case
            if xc(1)<=x(j) && x(j)< xc(2) %mi trovo in [x0,x1) left
            y(j)=3.*((x(j)-xc(1))./(xc(2)-xc(1))).^2-2.*((x(j)-xc(1))./(xc(2)-xc(1))).^3;
            end
        elseif l==m %second special case
              if xc(m)<x(j) && x(j)<= xc(m+1) %xc(m) is xm-1, xc(m+1) is xm
                  %mi trovo in (xm-1,xm] right
            y(j)=3.*((xc(m+1)-x(j))./(xc(m+1)-xc(m))).^2-2.*((xc(m+1)-x(j))./(xc(m+1)-xc(m))).^3;
              end
        else
             if xc(l)<x(j) && x(j)<= xc(l+1) %right
            y(j)=3.*((xc(l+1)-x(j))./(xc(l+1)-xc(l))).^2-2.*((xc(l+1)-x(j))./(xc(l+1)-xc(l))).^3;
        
             elseif xc(l+1)<=x(j) && x(j)< xc(l+2) %left
              y(j)=3.*((x(j)-xc(l+1))./(xc(l+2)-xc(l+1))).^2-2.*((x(j)-xc(l+1))./(xc(l+2)-xc(l+1))).^3;
             end
        end 
    end
end

function y=fS(x,l,xc)
m=length(xc)-1;
y=zeros(size(x));
    for j=1:length(x)
        if l==0 %first special case
            if xc(1)<=x(j) && x(j)< xc(2) %xc(1) is x0, xc(2) is x1
            y(j)=-(x(j)-xc(1)).^2./(xc(2)-xc(1))+(x(j)-xc(1)).^3./(xc(2)-xc(1)).^2;
            end
        elseif l==m %second special case
              if xc(m)<x(j) && x(j)<= xc(m+1) %xc(m) is xm-1, xc(m+1) is xm
                  %ho cambiato l'uguale perchè il valore della funzione era uguale
                  %a 0
            y(j)=((xc(m+1)-x(j)).^2./(xc(m+1)-xc(m)))-((xc(m+1)-x(j)).^3./(xc(m+1)-xc(m)).^2);
              end
        else
             if xc(l)<x(j) && x(j)<= xc(l+1) 
            y(j)=((xc(l+1)-x(j)).^2./(xc(l+1)-xc(l)))-((xc(l+1)-x(j)).^3./(xc(l+1)-xc(l)).^2);
        
             elseif xc(l+1)<=x(j) && x(j)< xc(l+2) 
                 y(j)=-(x(j)-xc(l+1)).^2./(xc(l+2)-xc(l+1))+((x(j)-xc(l+1)).^3./(xc(l+2)-xc(l+1)).^2);
             end
        end 
    end
end

function y=dH(x,l,xc)
m=length(xc)-1;
y=zeros(size(x));
    for j=1:length(x)
        if l==0 %first special case
            if xc(1)<=x(j) && x(j)< xc(2) %xc(1) is x0, xc(2) is x1
            y(j)= 6 .* (x(j)-xc(1)) ./ (xc(2)-xc(1)).^2 + ...
                                6 .*  (x(j)-xc(1)).^2 ./ (xc(1)-xc(2)) .^3;
            end
        elseif l==m %second special case
              if xc(m)<x(j) && x(j)<= xc(m+1) %xc(m) is xm-1, xc(m+1) is xm
        y(j)= - 6 .* (xc(m+1)- x(j)) ./ (xc(m+1)-xc(m)).^2 + ...
                                6 .*  (xc(m+1)-x(j)).^2 ./ (xc(m+1)-xc(m)) .^3;
              end
        else
             if xc(l)<x(j) && x(j)<= xc(l+1) 
             y(j)= - 6.*(xc(l+1)-x(j))./(xc(l+1)-xc(l)).^2+ ...
                 6.*(xc(l+1)-x(j)).^2./(xc(l+1)-xc(l)).^3;
        
             elseif xc(l+1)<=x(j) && x(j)< xc(l+2) 
         y(j)= 6 .* (x(j)-xc(l+1)) ./ (xc(l+2)-xc(l+1)).^2 + ...
                                6 .* ( x(j)- xc(l+1)).^2 ./ (xc(l+1)-xc(l+2)) .^3;
             end
        end 
    end
end

function y=dS(x,l,xc)
m=length(xc)-1;
y=zeros(size(x));
    for j=1:length(x)
        if l==0 %first special case
            if xc(1)<=x(j) && x(j)< xc(2) %xc(1) is x0, xc(2) is x1
            y(j)= 2 .*  (x(j)-xc(1))  ./ (xc(1)-xc(2)) + ...
                                 3 .* ( x(j)-xc(1) ).^2 ./ (xc(1)-xc(2)) .^2;
            end
        elseif l==m %second special case
              if xc(m)<x(j) && x(j)<= xc(m+1) %xc(m) is xm-1, xc(m+1) is xm
           y(j)= 2 .* (xc(m+1)-x(j)) ./(xc(m)-xc(m+1)) +...
                                3 .*  (xc(m+1)-x(j)).^2 ./ (xc(m+1)-xc(m)).^2;
              end
        else
             if xc(l)<x(j) && x(j)<= xc(l+1) 
          y(j)= 2 .* (xc(l+1)-x(j)) ./(xc(l)-xc(l+1)) +...
                                3 .* ( xc(l+1)-x(j)).^2 ./ (xc(l+1)-xc(l)).^2;
                        
             elseif xc(l+1)<=x(j) && x(j)< xc(l+2)  
            y(j)= 2 .* (x(j)-xc(l+1)) ./ (xc(l+1)-xc(l+2))  + ...
                                 3 .* (x(j)-xc(l+1)).^2 ./  (xc(l+1)-xc(l+2)).^2;
            end
        end 
    end
end

function y=fu(x,l,xc)
%selection function
y=zeros(size(x));
    if mod(l,2)==1 %dispari, verifico il resto della divisione per 2
        y=fS(x,(l-1)/2, xc);
    else
        y=fH(x,l/2, xc);
    end
    end
    
    function y=du(x,l,xc)
    y=zeros(size(x));
    if mod(l,2)==1 %dispari
        y=dS(x,(l-1)/2, xc);
    else
        y=dH(x,l/2, xc);
    end
end

function [lb,rb]=support(l,xc)
m=length(xc)-1;
    if l==1 %case s0
        lb=xc(1);
        rb=xc(2);
    elseif l>=2*m
        lb=xc(m);
        rb=xc(m+1);
        
    elseif mod(l,2)==1 %caso dispari
        lb=xc((l-1)/2); %(l-1)/2 because of the matlab translation
        rb=xc((l-1)/2+2);
    else
         %caso pari
         lb=xc(l/2);
         rb=xc(l/2+2);
    end
end

function y=fA(x,l,j,xc)
%we may write in therms of fu and du, not with fS and fH
y=(1+x).*du(x,l,xc).*du(x,j,xc)+((2+x)./(1+x)).*fu(x,l,xc).*fu(x,j,xc);
end

function y=fB(x,l,xc)
y=fu(x,l,xc).*(6/5-1/5.*x.*(2+x)./(1+x));
end

[Aurora]

Ho un problema con il supporto nel caso in cui si hanno due condizioni essenziali, per caso potresti aiutarmi?

Answer 1

[Ishu]

Hi Francesca,

As in your code, you have created a 'for loop' for 'l = 1:2*m+1' for the calculation of 'u'. And for every iteration you are calling 'fu()' function. So for 'l = 2*m' the function 'fH()' within the 'fu()' will be called and for 'l = 2*m+1' the function 'fS()' within the 'fu()' will be called. To correct the 'u' value you have to make changes in these two functions.

In these functions you have a special case with 'l == m' :

When 'j = lenght(x)', in this case "x(j) == xc(m+1)" and hence 'y(j)' becomes 0, because of which your last approximation is wrong. Therefore, in the approximation instead of 'xc(m+1)' you can use 'xc(m)'.

You can change your code as:

For fH()

elseif l==m %second special case
    % make a little change in if condition
              if xc(m)<x(j) && x(j)< xc(m+1) %xc(m) is xm-1, xc(m+1) is xm
                  %mi trovo in (xm-1,xm] right
            y(j)=3.*((xc(m+1)-x(j))./(xc(m+1)-xc(m))).^2-2.*((xc(m+1)-x(j))./(xc(m+1)-xc(m))).^3;
            % you can add the following code
              elseif x(j) == xc(m+1)
                  y(j)=3.*((xc(m)-x(j))./(xc(m+1)-xc(m))).^2+2.*((xc(m)-x(j))./(xc(m+1)-xc(m))).^3;
              end

For fS()

elseif l==m %second special case
    % make a little change in if condition
              if xc(m)<x(j) && x(j)< xc(m+1) %xc(m) is xm-1, xc(m+1) is xm
                  %ho cambiato l'uguale perchè il valore della funzione era uguale
                  %a 0
            y(j)=((xc(m+1)-x(j)).^2./(xc(m+1)-xc(m)))-((xc(m+1)-x(j)).^3./(xc(m+1)-xc(m)).^2);
            % add the following code
              elseif x(j) == xc(m+1)
                  y(j)=((xc(m)-x(j)).^2./(xc(m+1)-xc(m)))-((xc(m)-x(j)).^3./(xc(m+1)-xc(m)).^2);
              end              

Hope it helps!

[Francesca]

Thank you. You were very helpful