You are getting a warning about the matrix condition.
How did you verify this isn't causing the problem?
NaN Error , while solving Newton Raphson method
6 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
% 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 Kommentare
Antworten (1)
Torsten
am 12 Okt. 2023
Bearbeitet: Torsten
am 12 Okt. 2023
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 Kommentare
MAYUR MARUTI DHIRDE
am 12 Okt. 2023
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. 
Torsten
am 12 Okt. 2023
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
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
Sie können auch eine Website aus der folgenden Liste auswählen:
Amerika
- América Latina (Español)
- Canada (English)
- United States (English)
Europa
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asien-Pazifik
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)