Solving a problem with due essential boundery conditions with the Non Linear Collocation Method

I tried to implement the nonlinear collocation method to solve a differential equation, but the exact solution recommended by the professor does not coincide with the numerical solution.

Differential equation:

-y'' + exp(-x) * y * y' = exp(x) * (x^2 - 2)

Boundary conditions:

y(-1) = -1/exp

y(1) = exp

The exact solution is:

y=x.*(exp(x))

The graphic should coincides, but they don't.

Can you help me?

Thank you in advance.

%For testing 
[x,u]=NonLinearCollocationEssential(5,25,0,0.01);

function [x, u] = NonLinearCollocationEssential(m, nv, d0, tol)
% Non-linear collocation method to approximate the solution of the differential equation
% -y'' + exp(-x) * y * y' = exp(x) * (x^2 - 2)
% with essential boundary conditions y(-1) = -1/exp(1), y(1) = exp(1)

% Input
% d0: initial value for delta0 for the iteration of the method
%     (Attention to the choice of d0)
% nv: number of visualization points
% m: number of elements in the basis
% tol: tolerance for the stopping criteria

%Equation:
% -y'' + exp(-x) * y * y' = exp(x) * (x^2 - 2)

%BC:
%y(-1) = -1/exp(1) 
%y(1) = exp(1)

%Exact solution: y=x.*(exp(x));

%For testing 
%[x,u]=NonLinearCollocationEssential(5,25,0,0.01);
%y=x.*(exp(x));
%plot(x,u,'r*',x,y,'k-');

hc = 2 / (m + 2); % Collocation step (size of interval / m+2)
xcol = -1 + hc : hc : 1 - hc; % Collocation points (a + hc, b - hc)
delta = d0 * ones(m + 1, 1); % Initial values of delta, we have m+1 elements
ndd = 1 + tol; % Norm of the delta difference, to enter the cycle


while ndd > tol  % Ciclo del metodo di Newton
    J = zeros(m + 1); % Matrice Jacobiana di dimensione m+1
    epsilon = zeros(m + 1, 1); % Vettore epsilon

    for j = 1:m + 1
        for l = 0:m
            J(j, l + 1) = -ddu(xcol(j), l) + exp(-xcol(j)) * (fu(xcol(j), l) * dU(xcol(j), delta) + ...
                du(xcol(j), l) * U(xcol(j), delta));
            % l inizia da 0, quindi consideriamo la colonna l+1
        end
        epsilon(j) = -ddU(xcol(j), delta) + exp(-xcol(j)) * U(xcol(j), delta) * dU(xcol(j), delta) - ...
            exp(xcol(j)) * (xcol(j)^2 - 2); 
    end

    deltanew = J \ epsilon; % Calcolo della soluzione lineare
    ndd = norm(delta - deltanew); % Calcolo della norma della differenza
    delta = deltanew; % Aggiornamento di delta
end

% Evaluate the solution computing u
h = 2 / nv; % |a + b| / nv
x = -1 : h : 1;

u = U(x, delta);

% BC
u(1,1) = -1/exp(1); % y(-1) = -1/exp(1)
u(end,end) = exp(1); % y(1) = exp(1)

% Plot the results for testing
y = x .* exp(x);
plot(x, u, 'r*', x, y, 'k-');
xlabel('x');
ylabel('y');
legend('Approximate Solution', 'Exact Solution');
title('Non-Linear Collocation Method');
end


function y = fu(x, l)
    % Function for the elements
    y = x .^ l .* (x - 1) .* (x + 1);
end

function y = du(x, l)
    % First derivative of the base
    y = (x .^ (l-1)) .* (l.* (x.^2 -1) +2*x.^2);
end

function y = ddu(x, l)
    % Second derivative of the base
   y = (x .^ (l-2)) .* (l^2 .* (x.^2 -1) + 3*l*x.^2 + l +2*x.^2);
end


function y = U(x, delta)
    % U is the approximated solution
    % delta is the coefficient
    u = zeros(1, length(x));
    m = length(delta) - 1; % m+1 components

    for l = 1:m
        u = u + delta(l + 1) .* fu(x, l);
 
    end
    y=u;

end

function y = dU(x, delta)
    u = zeros(1, length(x));
    m = length(delta) - 1; % m+1 components

    for l = 1:m
        u = u + delta(l + 1) .* du(x, l);
    end

    y = u;
end

function y = ddU(x, delta)
    u = zeros(1, length(x));
    m = length(delta) - 1; % m+1 components
    
    for l = 1:m
        u = u + delta(l + 1) .* ddu(x, l);
    end

   y = u;
end

[John D'Errico]

Without even a glance at your code, I would suggest the plot indicates your problem is a failure of the boundary conditions. That is, you made a mistake at the end points in some form. Interior to the curve, the two solutions nearly overlap. So look closely at your BCs and how you implemented them.

[Torsten]

Your code seems to impose the boundary conditions y(-1) = y(1) = 0. But y(-1) = -1/exp(1), y(1) = exp(1) is required.

[Francesca]

The conditions at the edge have been fixed. But the points that approach the extremes are far from the exact solution.

[Sam Chak]

Hi @Francesca

The approximate solution depends on U(x, delta). If you check the local function U, you will see that it initializes 'u' as follows:

u = zeros(1, length(x));

Additionally, there's a loop that depends on 'fu(x, l)' and ultimately, 'y' is assigned the value of 'u':

for l = 1:m
    u = u + delta(l + 1).*fu(x, l);
end
y = u;

Upon further examination of 'fu(x, l),' you will find this formula:

y = x.^l.*(x - 1).*(x + 1);

where it shows when because .

Consequently, the initial value of 'U(x, delta)' can be simplified as follows:

0 + delta(l + 1).*(0) = 0.

[Sam Chak]

Also notice that these two computed values are unused in the local function U(x, delta):

a = (+ 1 + exp(2))/(2*exp(1));
b = (- 1 + exp(2))/(2*exp(1));

On which lines did you implement the boundary conditions?

Try pressing "Control+F" (or "Command+F" on a Mac) to find exp(1) in the code.

[Francesca]

I tried inserting the boundary conditions inside the U function (even outside the loop), but the code crashes.

[Sam Chak]

Hi @Francesca, have you figured out how to address the issue of the boundary conditions?