Why do I get NaN in Newton-Raphson

Hi,

I'm wiritng a script for solving a given equation with boundary conditions by using Finite Differences Method and Newton-Raphson Method. I'm quite confused as in the final while loop as a result I receive NaN. I don't know completly what I'm doing wrong. I would be greateful if sombedy checked my code.

The equation I'm solving: u''(x) + [1+(u(x))^2]*u'(x) = f(x), boundary conditions : u(0)=0 and u'(x) = 2*pi*cos(2*pi*x). I'm solving it for U vector.

N = input("set the number of nodes to use: ");
h = 1/N;
x_start = 0;
x_end = input("set the upper boundary in which to consider the function: ");

xj = linspace(x_start,x_end,N+2);

% N- is a number of nodes, for each node will be one equation
% x_start- first x out of xj
% x_end- last x out of xj
% xj- linear space in which the task is considered.
%     It is important to note that equations will be created for xj from
%     xj(2) to xj(N+1). Where xj(2) it is j1 and the boundary condition is
%     applied so this part will be hard coded here. Where we have x(N+1) it
%     is j_LAST and here also the upper boundary condition will be
%     applied, this time u'(x)=beta ==> 
%     => u_n = u_{n-1} + h*u'(1), where u_n is an artificial(extra) node!

%boundary conditions: u(0) = alfa = 0
%                     u'(N+1) = 2*pi*cos( 2*pi*x_{N+1} )
BETA = 2*pi*cos( 2*pi*xj(N+1) );

% right side function: f(x) = -4*pi^2 *sin( 2*pi*xj ) + 2*pi*(cos( 2*pi*xj ))^2

% bj = (cos( 2*pi*xj ))^2

% creating vectors: bj, fj. Pre-declared values for later modifications
bj = zeros(N,1);
fj = zeros(N,1);

for i=1:N
    bj(i) = (cos(2*pi*xj(i+1)))^2;
    fj(i) = -4*pi^2 *sin( 2*pi*xj(i+1) ) + 2*pi*(cos( 2*pi*xj(i+1) ))^2;
end

% creating matrix for function: u_{j-1}(1-0.5hbj)+u_{j+1}(1+0.5hbj)-2uj
x0 = zeros(N,1); %vector with values for initialization
F = zeros(N,1);
U = zeros(N,1); %vector U- values: here U(0)=0(ALWAYS)
%                  U(N+1)=U(N) + h*2*pi*cos( 2*pi*xj(N+1) )
U = update_U(U,x0,N); %setting U to its first value

%specyfing vector U according to the comment above
%U(N+1) = U(N) + h*2*pi*cos( 2*pi*xj(N+1) );

%setting values for function
%function's value in the first is different and cannot be looped
F(1) = U(2)*(1 + 0.5*h*bj(1)) - 2*U(1);

%function's value in between 2 and N-1
for i=2:N-1
    F(i) = U(i-1)*(1 - 0.5*h*bj(i)) + U(i+1)*(1 + 0.5*h*bj(i)) - 2*U(i);
end

%function's value in the last equation is different and cannot be looped
F(N) = U(N-1)*(1 - 0.5*h*bj(N)) + U(N)*(0.5*h*bj(N) - 1) + h*BETA*(0.5*h*bj(N) + 1);

%computing the Jacobian for function F:

%METHOD: [J_f(u)]_{i,j} = 1/epsilon *[f_i(u1,u2...u_{j+epsilon}..u_N) - f_i(u1,u2...u_j..u_N)]
%        where "epsilon" is smth around 0.00001(MUST be much smaller than "h")
%        scheme above is just a single cell of this matrix, in fact each
%        row is for function F(j) and each cell in this row is the same
%        formula but for another element from vector U

Jacobian = zeros(N,N); %initialization of Jacobian matrix for derivatives in certain vector U
epsilon = 0.001;
%the loop below fills matrix of size (i,j) where each "i" is a row and each
%"j" is a column. Each row represents function that is in vector F and
%coulmn is a variable to derive: dF(i)/dU(j). Where the function
%"calculate_exact_function()" takes the last argument as "true" an extra
%value is added to correlating U(i).


for i=1:N
    for j=1:N
        Jacobian(i,j) = (calculate_exact_function(U,F,bj,h,BETA,i,N,true) - calculate_exact_function(U,F,bj,h,BETA,i,N,false))/epsilon;
    end
end

%Newton-Raphson method:
stop_cryterium = 0.00001; %it will say whether the precision of results is enough
iteration_number = 0;

while iteration_number < 11 
    %all values F,U,fj,bj are initialized int the code above
    P = Jacobian\F;
    x0 = x0 - P;
    
    %updating vector U
    U = update_U(U,x0,N);
    %updating vector F
    %function's value in the first is different and cannot be looped
    F(1) = U(2)*(1 + 0.5*h*bj(1)) - 2*U(1);

    %function's value in between 2 and N-1
    for i=2:N-1
        F(i) = U(i-1)*(1 - 0.5*h*bj(i)) + U(i+1)*(1 + 0.5*h*bj(i)) - 2*U(i);
    end

    %function's value in the last equation is different and cannot be looped
    F(N) = U(N-1)*(1 - 0.5*h*bj(N)) + U(N)*(0.5*h*bj(N) - 1) + h*BETA*(0.5*h*bj(N) + 1);
    %updating Jacobian vector
    for i=1:N
        for j=1:N
            Jacobian(i,j) = (calculate_exact_function(U,F,bj,h,BETA,i,N,true) - calculate_exact_function(U,F,bj,h,BETA,i,N,false))/epsilon;
        end
    end
   
    iteration_number = iteration_number + 1;
end

disp(x0);

function f = calculate_exact_function(U,F,bj,h,BETA,i,N,eps)
%U- vector of U values
%F- function vector, exact F(i) will be modified as certain U(i) will get
%its value modified like => U(i) = U(i) + eps
%bj- coefficient calculated above
%h- coefficient calculated above
%BETA- coefficient calculated above
%i- informs which iteration is now; necessary to make this function knows
%to which element from vector U add "eps" correction
%eps- for approximation. In this function eps is a volountary argument, if
%not passed then this function returns simply F(i) unmodified!!!
if eps == true
    epsilon = 0.001;
    U(i) = U(i) + epsilon;
    
    if i == 1
        F(1) = U(2)*(1 + 0.5*h*bj(1)) - 2*U(1);
        f = F(1);
    elseif i == N
        F(N) = U(N-1)*(1 - 0.5*h*bj(N)) + U(N)*(0.5*h*bj(N) - 1) + h*BETA*(0.5*h*bj(N) + 1);
        f = F(N);
    elseif 2 <= i && i <=(N-1)
        F(i) = U(i-1)*(1 - 0.5*h*bj(i)) + U(i+1)*(1 + 0.5*h*bj(i)) - 2*U(i);
        f = F(i);
    end
    
elseif eps == false
    if i == 1
        F(1) = U(2)*(1 + 0.5*h*bj(1)) - 2*U(1);
        f = F(1);
    elseif i == N
        F(N) = U(N-1)*(1 - 0.5*h*bj(N)) + U(N)*(0.5*h*bj(N) - 1) + h*BETA*(0.5*h*bj(N) + 1);
        f = F(N);
    elseif 2 <= i && i <=(N-1)
        F(i) = U(i-1)*(1 - 0.5*h*bj(i)) + U(i+1)*(1 + 0.5*h*bj(i)) - 2*U(i);
        f = F(i);
    end
end

end
function U = update_U(U,x_temp,N)
%U- vector of U values(u1,u2....u_N)
%x_temp- temporary values taken from Newton-Raphson method
%U = U + x_temp;
for j=1:N
    U(j) = x_temp(j);
end
end

[J. Alex Lee]

Before going down the rabbit hole, I think you should clarify if you are doing this to solve a problem in real life, or as a learning exercise (because this method is probably not the best/most accurate/most efficient way to solve).

The operations that result in NaN are things like 0/0, operations on Infs (which can be generated by division by zero of anything other than zero).

Do you also get a warning that says your matrix is singular to working precisions?

[Szymon Bogus]

This is a project for on of my university course. This equation as well as the methodology was given to me my professor. All I'm doing is implementing what my professor told me to implement.

"Do you also get a warning that says your matrix is singular to working precisions?": Yes I do.

I know that NaN appears when there is 0/0 but it seems to me I didn't make any mistake coding this algorithm. Function vector F and Jacobian update every iteration according to new values of U. Maybe the problem lays in incorrect initializing value of vectro x0 which is zeros(N,1).

[Walter Roberson]

When people get NaN doing this kind of exercise, most of the time it is triggered by the \ operator, which in turn has been asked to work with large values . In Jacobian\F the F probably has values close to +/- 1E308 at the time the problem first occurs. A bit of multiplication and addition leads to overflows to +/- infinity, then that leads within a step or two to trying to subtract infinities, giving NaN.

Why the values explode to about +/- 1e308 is a different question. Sometimes the system is just unstable. Sometimes poor initial conditions were chosen. Sometimes there is implementation error, such as the jacobian turning out to have been programmed incorrectly. Sometimes the system is discontinuous -- these kinds of NR approaches are only valid on systems that are continuously differentiable anywhere near the system is likely to go.

        F(N) = U(N-1)*(1 - 0.5*h*bj(N)) + U(N)*(0.5*h*bj(N) - 1) + h*BETA*(0.5*h*bj(N) + 1);
        f = F(N);

Why do you bother to assign to F(N) instead of directly to f ? F is an input vector, but you are not returning the updated F as an output.

[Szymon Bogus]

@Walter Roberson yeah you're right , that might be unnecessary here to assign it in this way. I'll check what happens when I fix it to f(N).

[Szymon Bogus]

@Walter Roberson I applied those changes so it is right now like below, however this time erro like this one occured.

f(N) = U(N-1)*(1 - 0.5*h*bj(N)) + U(N)*(0.5*h*bj(N) - 1) + h*BETA*(0.5*h*bj(N) + 1);

[J. Alex Lee]

Along the lines of what Walter explained, I can suggest a few diagnostic tools:

• monitor the NR iterations
• are you getting the errors on the very first evaluation? maybe your initial guess is bad
• are you getting the errors only after several iterations? maybe your jacobian is misprogrammed to take your solution far away, or maybe your initial guess is in a bad "neighborhood"
• calculate a numerical jacobian to check your jacobian against

TO elaborate on the efficacy of NR and validity, in a lot of real world cases you can find solutions to tweaked versions of the problem of interest, so in principle you can always start from a known solution and continuously change the problem using the solutions as initial guesses toward the problem that you care about, and NR should work well as long as you take small enough "steps" so that your old solutions (as initial guesses) remain close to the region of convergence. Of course many problems of real world interest will also fail with this methodology when you have more interesting behavior like critical points. Long story short, NR is a useful workhorse for these types of problems, I suggest to be favored over more "robust" methods that are geared toward more generalized applications.