Solving a problem with due natural 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) + y(1) = 3 * exp

The exact solution is:

y=x.*(exp(x))

For testing

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

The graphs should coincide, but they don't.

Can you help me?

My code is:

function [x,u]=nlCOLL(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 boundary conditions y(-1) = -1/exp(1) and y'(1) + y(1) = 3 * 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

% For the natural boundary condition, the Jacobian has an element equal to 1
% in the bottom right corner, and the epsilon vector has an element equal
% to U(1, delta) in the last position.

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 % Cycle of the Newton method
    J = zeros(m + 2); % Jacobian matrix of order m+2
    e = zeros(m + 2, 1); % Epsilon vector
      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 starts from 0, so we can consider the column l+1
            end
      J(j, m + 2) = fu(xcol(j), m + 1);
      % Additional entry for the natural boundary condition
      e(j) = -ddU(xcol(j), delta) + exp(-xcol(j)) * U(xcol(j), delta) * dU(xcol(j), delta) - ...
      exp(xcol(j)) * (xcol(j)^2 - 2);
      end
    % Add the new boundary conditions to the Jacobian and the epsilon vector
    J(m + 2, m + 1) = du(1, 0) + U(1, delta); % new boundary condition
    J(m + 2, m + 2) = 1; % natural condition
    e(m + 2) = U(1, delta) - 3 * exp(1); % new boundary condition
    deltanew = J(1:m+1, 1:m+1) \ e(1:m+1); % is the linear solution
    ndd = norm(delta - deltanew); % compute the norm of the difference
    delta = deltanew; % update the delta
    end
% Evaluate the solution computing u
h = 2 / nv; % |a + b| / nv
x = -1 : h : 1;
u = U(x, delta);
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=(l+2).*x.^(l+1)-(l).*x.^(l-1);
end

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

function y = U(x, delta)
% U is the approximated solution
% delta is the coefficient
u = zeros(1, length(x));
a = exp(1) - (exp(1) - 1 / exp(1)) / 2;
b = (exp(1) - 1 / exp(1)) / 2;
    m = length(delta) - 1; % m+1 components
    for l = 1:m
    u = u + delta(l + 1) .* fu(x, l);
    end
y =  u; % u(x) = omega(x) + sum(delta(l) * u(l)(x))

end

function y = dU(x, delta)
% dU/dx is the derivative of the approximated solution U
% delta is the coefficient
u = zeros(1, length(x));
a = exp(1) - (exp(1) - 1 / exp(1)) / 2;
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)
% ddU/ddx^2 is the second derivative of the approximated solution U
% 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) .* ddu(x, l);
    end
y = u;
end

[Shivansh]

Hi Francesca,

The code provided by you is not working. I am getting the following error on R2023a.

Can you provide the Matlab code file and the version of your MATLAB? It will be useful to reproduce the issue you are facing and getting a possible workaround.

[Francesca]

I have this problem too. I can't understand the errors.

[Torsten]

From function U, you call function dU in the line

y(end) = y(end) + dU(1, delta) - 3 * exp(1);

from function dU, you call back function U in the line

y(end) = y(end) + U(1, delta) - 3 * exp(1);

And this ping-pong game won't stop.

[Sam Chak]

Hi @Francesca

It appears that the nlCOLL() function is invoking other functions which utilize the input value 'u' to compute the output value 'u', thus creating a circular dependency that results in the Ouroboros anomaly. That's how I interpreted "infinite recursion within the program".

u = u + delta(l + 1).*fu(x, l);

[Francesca]

How can I eliminate this error so that the method works?

[Torsten]

After knowing the reason why the code fails, it should be easier for you as the code writer to repair the problem, shouldn't it ? Or aren't you the code writer ?

[Francesca]

This is a code given by the teacher to be corrected and finished at home

[Sam Chak]

@Francesca, So, in addition to solving the given boundary value problem using the custom algorithm nlCOLL(), are you also required to fix the bug in the nlCOLL() code?

[Francesca]

The task is to solve this problem with this code model, possibly to be corrected if necessary.

[Torsten]

Since it's very hard to understand what another person has coded, our teacher told us to better start anew.

It will usually take less time.

[Francesca]

I've modified the code but still have problems at the extremes of the range.