Trying to create code which solves a BVP using a second-order finite difference scheme and Newton-Raphson method

% Code to calculate the solution to a 
% nonlinear BVP with Dirichlet and Robin BCs

% Some physical parameter.
mu=10;

% Define the interval over which solution is calculated.
a=0; b=1;

% Define N.
N=500;

% Define the grid spacing.
h=(b-a)/N;

% Define the grid. Note that there is no need to solve at x=a.
x = reshape(linspace(a+h,b,N),N,1);

% Define an initial guess for the solution 
% of the BVP (make sure it is a column vector)
U=1.0 - 0.2*x
;

% Define a tolerance for the termination of Newton-Raphson.
tol=10^(-8);

% Ensure that F is such that at least one iteration is done.
F=ones(N,1);
J=zeros(N,N);

% Store the initial guess in SOL.
SOL= U;

% Loop while the norm(F) is greater than tol.
while (norm(F)>tol)

%% Define F(u).
% Boundary conditions.
F(N)=U(N-2)-4*U(N-1)+3*U(N)+2*h*(U(N))^3;

% Finite difference approximation to ODE at interior nodes.
F(1)=2*(U(2)+1-2*U(1)+h*exp(U(1))*(U(2)-1)-2*h^(2)*mu*sin(2*pi*x(1)));
for i=2:N-1
   F(i)=2*(U(i+1)+U(i-1)-2*U(i)+h*exp(U(i))*(U(i+1)-U(i-1))-2*h^(2)*mu*sin(2*pi*x(i)));
end

%% Define the Jacobian J.
% First row corresponds to BC at x=0.
J(1,1)=(-4)+h*exp(U(1))*(U(2)-1);
J(1,2)=2+h*exp(U(1));

% Last row corresponds to BC at x=1.
J(N,N-2)=1;
J(N,N-1)=-4;
J(N,N)=3+6*h*(U(N))^2;

% Intermediate rows correspond to F(i)=...
for i=2:(N-1)
% Diagonal entries.
  J(i,i)=(-4)+h*exp(U(i))*(U(i+1)-U(i-1));
% There may be off-diagonal entries, check your calculation.
  J(i,i-1)=2-h*exp(U(i));
  J(i,i+1)=2+h*exp(U(i));
end

%% Having defined F and J, update the approximate solution to
% the difference equations.
U=U-J\F ;

%% Store the new approximation.
SOL=[SOL,U];

end

%% Insert your plot commands here.
figure;
plot(x, SOL(:, 1), 'r--', 'LineWidth', 1.5); hold on;
for k = 2:size(SOL,2)
    plot(x, SOL(:, k), 'LineWidth', 1.5);
end
xlabel('x');
ylabel('U(x)');
title('Convergence of Newton-Raphson Iterations');
legend('Initial Guess', 'Iterations');
grid on;

This is my code for trying to solve the BVP

u'' + exp(u)u' = mu*sin(2*pi*x), u(0) = 1, u'(1) + (u(1))^3

with a second-order finite difference scheme and the Newton-Raphson method, I was given unfinished code and this is my finished version. When I run the code and it plots the iterations I can clearly see something isn't right as I know what the end result should approximately look like, but I can't see where I've gone wrong. Also, I cannot use certain inbuilt MATLAB functions, only what is already in the code.

Any help?

[Ole]

I believe that my jacobian matrix J should be Tridiagonal but isn't, is that perhaps the issue?

[Torsten]

For comparison with your approximations. It's always good to know what you are aiming at.

mu = 10;
xmesh = linspace(0,1,20);
solinit = bvpinit(xmesh, @(x)[1-0.2*x;-0.2]);
fun = @(x,y)[y(2);-exp(y(1))*y(2)+mu*sin(2*pi*x)];
bc = @(ya,yb)[ya(1)-1;yb(2)+yb(1)^3];
sol = bvp4c(fun, bc, solinit);
plot(sol.x,sol.y(1,:))

Answer 1

[Torsten]

% Code to calculate the solution to a 
% nonlinear BVP with Dirichlet and Robin BCs
% Some physical parameter.
mu = 10;
% Define the interval over which solution is calculated.
a = 0; b = 1;
% Define N.
N = 50;
% Define the grid spacing.
h = (b-a)/(N-1);
% Define the grid. Note that there is no need to solve at x=a.
x = linspace(a,b,N).';
% Define an initial guess for the solution 
% of the BVP (make sure it is a column vector)
U = 1.0 - 0.2*x;
% Define a tolerance for the termination of Newton-Raphson.
tol = 1e-8;
% Ensure that F is such that at least one iteration is done.
F = ones(N,1);
J = zeros(N,N);
% Store the initial guess in SOL.
SOL = U;
% Loop while the norm(F) is greater than tol.
while (norm(F)>tol)
%% Define F(u).
% Boundary conditions.
F(1) = U(1)-1.0;
F(N) = (-2*h*U(N)^3-2*U(N)+2*U(N-1))/h^2-exp(U(N))*U(N)^3-mu*sin(2*pi*x(N));
% Finite difference approximation to ODE at interior nodes.
for i=2:N-1
   F(i) = (U(i+1)-2*U(i)+U(i-1))/h^2 + exp(U(i))*(U(i+1)-U(i-1))/(2*h)-mu*sin(2*pi*x(i));
end
%% Define the Jacobian J.
% First row corresponds to BC at x=0.
J(1,1) = 1.0;
% Last row corresponds to BC at x=1.
J(N,N-1) = 2/h^2;
J(N,N) = -6*U(N)^2/h-2/h^2-exp(U(N))*U(N)^3-3*exp(U(N))*U(N)^2;
% Intermediate rows correspond to F(i)=...
for i=2:N-1
% Diagonal entries.
  J(i,i-1) = 1/h^2-exp(U(i))/(2*h);
% There may be off-diagonal entries, check your calculation.
  J(i,i) = -2/h^2+exp(U(i))*(U(i+1)-U(i-1))/(2*h);
  J(i,i+1) = 1/h^2+exp(U(i))/(2*h);
end
%% Having defined F and J, update the approximate solution to
% the difference equations.
U = U-J\F ;
%% Store the new approximation.
SOL = [SOL,U];
end
%% Insert your plot commands here.
figure;
plot(x, SOL(:, 1), 'r--', 'LineWidth', 1.5); hold on;
for k = 2:size(SOL,2)
    plot(x, SOL(:, k), 'LineWidth', 1.5);
end
xlabel('x');
ylabel('U(x)');
title('Convergence of Newton-Raphson Iterations');
legend('Initial Guess', 'Iterations');
grid on;

[Ole]

How did you get your F(N) from the boundary condition?

[Torsten]

Boundary condition:

(U(N+1)-U(N-1))/(2*h) + U(N)^3 = 0           (1)

Discretized equation:

F(N) = (U(N+1)-2*U(N)+U(N-1))/h^2 + exp(U(N))*(U(N+1)-U(N-1))/(2*h) - mu*sin(2*pi*x(N)) = 0

Now insert in F(N)

U(N+1) as U(N-1) - 2*h*U(N)^3

from (1) in the approximation of the 2nd derivative at x = 1 and

(U(N+1)-U(N-1))/(2*h) as -U(N)^3

from (1)

[Ole]

Ah I see you used the first order backwards-difference equation for the derivative of u, I need to ensure second-order accuracy so have to use the second-order backwards difference equation.

[Torsten]

All approximations of derivatives in my code are second-order accurate. I don't know where you see a first-order backwards difference.

(U(i+1)-U(i-1))/(2*h) is second-order accurate for the first derivative of U in x(i)

(U(i+1)-2*U(i)+U(i-1))/h^2 is second-order accurate for the second derivative of U in x(i)

[Ole]

In my notes it tells me that the second-order backwards difference formula, so for the right hand boundary 1 in this case, is

(U(N-2)-4*U(N-1)+3U(N))/2h for the first derivative of U at x(N)

I have attached my notes, please look at page 122, Remark 6.32

[Torsten]

You can also use the formula from your notes, but introducing the ghostpoint x(N)+h with value U(N+1) and approximating the first-order derivative in x(N) as (U(N+1)-U(N-1))/(2*h) as I did doesn't destroy the tridiagonal Jacobian and is also second-order accurate.

But make an attempt to use your formula - I think you already did it correctly in your former code.