NaN Error , while solving Newton Raphson method
显示 更早的评论
% Parameters
Nx = 6; % Number of grid points
dx = 2; % Grid spacing
initial_guess = 1.4; % Initial guess for m(1)
boundary_condition = 4.16; % Boundary condition at m(Nx)
% Initialize m, including m(1)
m = zeros(Nx, 1);
m(1) = initial_guess;
% Convergence tolerance
tolerance = 1e-6;
% Maximum number of iterations
max_iterations = 100;
% Initialize variables
iteration = 0;
converged = false;
% Newton-Raphson iteration
while ~converged && iteration < max_iterations
% Calculate the residual vector
R = continuity_residual(m, Nx, dx, boundary_condition);
% Calculate the Jacobian matrix
J = continuity_jacobian(m, Nx, dx);
% Solve for the update using the linear system
update = -J \ R;
% Update m
m = m + update;
% Check for convergence
if max(abs(update)) < tolerance
converged = true;
end
% Increment the iteration counter
iteration = iteration + 1;
end
% Display the result
if converged
disp('Converged to a solution:');
disp(['m = ', num2str(m')])
else
disp('Newton-Raphson did not converge');
end
% Continuity Residual Function
function R = continuity_residual(m, Nx, dx, boundary_condition)
R = zeros(Nx, 1);
% Calculate the residual for interior points
for i = 2:Nx - 1
R(i) = (m(i + 1) - m(i - 1)) / (2 * dx);%% Central discretization
end
% Handle residual at i = 1
R(1) = (m(2) - m(1)) / dx; %% Forward discretization to avoid ghost point i=0
% Handle boundary condition at i = Nx
R(Nx) = (m(Nx) - m(Nx - 1)) / dx + boundary_condition; %% Backward discretization to avoid ghost point i=Nx+1
end
% Continuity Jacobian Function
function J = continuity_jacobian(m, Nx, dx)
J = zeros(Nx, Nx);
% Calculate the Jacobian matrix for interior points
for i = 2:Nx - 1
J(i, i - 1) = -1 / (2 * dx);
J(i, i + 1) = 1 / (2 * dx);
end
end
ERROR : WHEN I RUN THIS CODE , THE RESIDUAL VECTOR SOMEHOW GETTING VALUES AS NAN AT ALL 6 GRID POINT.
2 个评论
Rik
2023-10-12
You are getting a warning about the matrix condition.
How did you verify this isn't causing the problem?
MAYUR MARUTI DHIRDE
2023-10-12
回答（1 个）
for i = 2:Nx - 1
J(i, i - 1) = -1 / (2 * dx);
J(i, i + 1) = 1 / (2 * dx);
end
Your loop starts with i = 2, thus the first row of the Jacobian matrix J is a zero row. So J is singular.
Even if you add
J(1,1) = -1/dx;
J(1,2) = 1/dx;
J(Nx,Nx) = 1/dx;
J(Nx,Nx-1) = -1/dx;
which can be deduced from R(1) and R(Nx), your system remains underdetermined.
Which boundary value problem do you try to solve with your code ?
6 个评论
MAYUR MARUTI DHIRDE
2023-10-12
I am trying to solve this two steady state compressible flow equation using Newton Raphson method. I tried solving both equation simultaneoulsly but i was facing issue , so now i am trying to apply Newton raphson method on continuity equation to check the algorithm but i am getting error in above code.
Solution is simple:
m(x) = m(x=0) = constant
and for P, see below.
Why do you need Newton-Raphson here ?
syms a b x P(x) P0
eqn = diff(P,x) + a/P + b*P == 0
cond = P(0) == P0;
dsolve(eqn,cond)
MAYUR MARUTI DHIRDE
2023-10-12
I have two boundary condition P(1) = 1.13*10^6 and m(Nx) = 4.16 , In the above equation both m and P are unknown. I discretize above equation using cental difference method. Consider Nx=6 , The unknow variables are P2,P3,P4,P5,P6 and m1,m2,m3,m4,m5. I was trying to solve both equation simultaneoulsy using newton Raphason method , by forming four function, function 1 for residual continuity vector , function 2 residual momentum vector , function 3 jacobian continuity matrix and function 4 jacobian momentum matrix. At the end delta_combined = - combined_Jacobian \ combined_residual;
Further to avoid ghost point at boundary i= for m(1) , i discretize continuity equation using forward discretization which gave me equation as m(1) =m(2) and to avoid ghost point at i=Nx+1 for P(Nx) , i discretize momentum equation using backward discretization , which gave me equation P(Nx) = sqrt(dx / 8) * sqrt(abs(8 * P(Nx) * P(Nx - 1) - 8 * f * m(Nx)^2 * c^2 * pi * D - 8 * P(Nx)^2 * g / c^2));
Torsten
2023-10-12
Are you forced to use this overcomplicated way of solving the equations ?
As said, the solution is a one-liner (see above).
MAYUR MARUTI DHIRDE
2023-10-12
Thing is i have to solve above equation for both steady and unsteady , in steady state equation, one is linear so m will have constant value, it is like whatever comes in goes out , m =4.16 for all grid point and then we can only solve for P in second equation as per the way you suggested. RIGHT?
In unsteady state , the equation changes over here , i need to use Newton raphson method somehow. 
Ok. Here is a code that works for the stationary continuity equation.
If you have further specific questions, you should open a new request.
% Parameters
Nx = 6; % Number of grid points
dx = 2; % Grid spacing
initial_guess = 1.4; % Initial guess for m(1)
boundary_condition = 4.16; % Boundary condition at m(Nx)
% Initialize m, including m(1)
m = zeros(Nx, 1);
m(1) = initial_guess;
% Convergence tolerance
tolerance = 1e-6;
% Maximum number of iterations
max_iterations = 100;
% Initialize variables
iteration = 0;
converged = false;
% Newton-Raphson iteration
while ~converged && iteration < max_iterations
% Calculate the residual vector
R = continuity_residual(m, Nx, dx, boundary_condition);
% Calculate the Jacobian matrix
J = continuity_jacobian(m, Nx, dx);
% Solve for the update using the linear system
update = -J \ R;
% Update m
m = m + update;
% Check for convergence
if max(abs(update)) < tolerance
converged = true;
end
% Increment the iteration counter
iteration = iteration + 1;
end
% Display the result
if converged
disp('Converged to a solution:');
disp(['m = ', num2str(m')])
else
disp('Newton-Raphson did not converge');
end
% Continuity Residual Function
function R = continuity_residual(m, Nx, dx, boundary_condition)
R = zeros(Nx, 1);
% Calculate the residual for interior points
for i = 2:Nx - 1
R(i) = (m(i + 1) - m(i - 1)) / (2 * dx);%% Central discretization
end
% Handle residual at i = 1
R(1) = (m(2) - m(1)) / dx; %% Forward discretization to avoid ghost point i=0
% Handle boundary condition at i = Nx
R(Nx) = m(Nx) - boundary_condition; %% Backward discretization to avoid ghost point i=Nx+1
end
% Continuity Jacobian Function
function J = continuity_jacobian(m, Nx, dx)
J = zeros(Nx, Nx);
% Calculate the Jacobian matrix for interior points
for i = 2:Nx - 1
J(i, i - 1) = -1 / (2 * dx);
J(i, i + 1) = 1 / (2 * dx);
end
J(1,1) = -1/dx;
J(1,2) = 1/dx;
J(Nx,Nx) = 1.0;
end
类别
在 帮助中心 和 File Exchange 中查找有关 Newton-Raphson Method 的更多信息
2023-10-12
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
