Newton's solver not converging for 1D nonlinear diffusion equation.

Hello all, I am having trouble with a newton solver that I have been working on over the last week. It is in the context of modelling 2D random with proliferation walks via column averaging but that's beside the point. The PDE that describes this interaction is where D is the diffusion (migration) terma and lambda is the nonlinear (proliferation).

I have checked the equations used for the Jacobian and the f vector a dozen times to the notes in class so I'm 99% sure that's not the issue.

An implicit Euler method is used for those interested.

Newton's method fails to converge when proliferation is 'turned on' (pp > 0). The code works fine for pp = 0 (compared to the random walk data). Whenever I try pp = 0.01 I get the error('Newton failed to converge\n') message.

Any help would be greatly appreciated :D, code below

%% 1D Non-Linear Diffusion Equaltion Solver using Newton's Method
clear all
clc

% Set up variables
deltax = 0.1;
tau    = 0.1;  % Time step
xlims  = 20;   % X boundaries
xIC    = 10;   % Lattice Sites with the IC

pp     = 0.01; % Probablity to proliferate
pm     = 0.5;  % Probability for a migration


Newtol = 0.01; % newton tolerance
xarr = linspace(-xlims,xlims,2*xlims/deltax+1) + xlims; % x array

I = length(xarr); % number of lattice point
f = zeros(I,1);   % the f vector used in newton iterations
J = zeros(I);     % Jacobian for newton iterations
N = zeros(I,1);   % Vector with averaged particle position
deltaN = N;       % deltax vector for newton iterations
N(round(xlims/deltax+1)-xIC:round(xlims/deltax+1)+xIC) = 40; % Set up casona IC [-inf,0] = 1; (0,Inf] = 0


Nold = N; % used for newton iterations

curtol      = Newtol + 1; % current tol, inital value must be larger than newtol
numsteps    = 0;          % newton steps
maxnewsteps = 10;         % max newton steps
cond        = 1;          % convergence condition

D = pm * deltax^2/(4*tau); % diffusion constant
lam = pp / tau; % Nonlinear constant


% Set Tmax and the timestep array
Tmax = 100;

% for all the timesteps
for i = 1:(Tmax / tau)
    cond = 1;     % reset convergence condition
    numsteps = 0; % reset newton step count

    % while newton has not converged
    while cond

        % Calculate the Jacobian and the f vector

        % Left Boundary
        J(1,1) = -1;
        J(1,2) = 1;
        f(1) = N(2) - N(1);

        % Right Boundary
        J(I,I) = 1;
        J(I,I-1) = -1;
        f(I) = N(I) - N(I-1);


        for n = 2:I-1 % for each internal lattice site

            J(n,[n-1, n+1]) =  D / (deltax^2);
            J(n,n) = - 1 / tau - 2 * D / (deltax^2)...
                + lam - lam * 2 * N(n);

            f(n) = - 1 / tau * (N(n) - Nold(n)) + D /(deltax^2) ...
                * (N(n-1) - 2*N(n) + N(n+1)) + lam * N(n) - lam * N(n) ^ 2;
        end

        % Solve J dN = -f using the tridiagonal matrix algorithm
        deltaN = J\(-f);

        % set N_old to the previous N state except for the first
        % iteration as it is the same as the IC
        if i > 1
            Nold = N;
        end
        % generate new N
        N = deltaN + N;
        fprintf(2,'Timestep %g\n',i)
        % Calculate error metric the 2 norm of the f vector
        curtol = norm(f)

        % if newton converged end the loop
        if (Newtol > curtol)
            cond = 0;
        end

        % Stop running the code if max newton iterations reached
        if numsteps >= maxnewsteps
            fprintf(2,'Timestep %g\n',i)
            error('Newton failed to converge\n')
        end

        % Newton  step counter
        numsteps = numsteps + 1;

    end

end
% Save final solution data


% Plotting
clf
xlabel('\fontsize{12}x')
ylabel('\fontsize{12}N')
plot(xarr,N,'linewidth',1.5)

shg

Answer 1

[darova]

See this madness

lam = 0.1;
m = 100;
n = 20;
t = linspace(0,0.1,m);
x = linspace(0,1,n);
dt = t(2)-t(1);
dx = x(2)-x(1);
% dt/dx/dx

N = zeros(m,n);
N(1,:) = sin(pi*x)/2;
h = surf(x,t,N);
for i = 1:m-1
    for j = 2:n-1
        N(i+1,j) = N(i,j) + dt/dx^2*(N(i,j+1)-2*N(i,j)+N(i,j-1)) + ...
            dt*lam*N(i,j)*(1-N(i,j));
    end
    set(h,'zdata',N)
    pause(0.1)
end

Answer 2

[David Goodmanson]

Hi Patrick,

The problem appears to be normalization. The diffusion term, being linear in N, does not care about normalization and comes out the same, but the logistic term certainly cares. The statement

N(round(xlims/deltax+1)-xIC:round(xlims/deltax+1)+xIC) = 40;

has 21 sites, all of value 40, in which case the logistic term lambda*N*(1-N) is huge. Setting the 40 to 1 gives good results for reasonable values of lambda. You don't want N to be greater than 1, so setting the 40 to something a bit less than 1 might be safer.