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

1 次查看（过去 30 天）
Ole
Ole 2024-12-2
编辑： Torsten 2024-12-4
% 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?
  2 个评论
Torsten
Torsten 2024-12-2
Ran in:
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,:))

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2024-12-2
编辑：Torsten 2024-12-2
Ran in:
% 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;
  6 个评论
显示 4更早的评论
Ole
Ole 2024-12-4
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
Torsten 2024-12-4
编辑：Torsten 2024-12-4
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2024b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by