Why are values in my loop filled with NaN?

8 visualizaciones (últimos 30 días)
Westin Messer
Westin Messer el 27 de Nov. de 2018
Editada: Stephen23 el 27 de Nov. de 2018
Hello,
Can anyone tell why 'u' and 'v' inside my loop are all being filled with NaN? I cant figure out what I'm doing wrong. Thanks!
close all;
clear all;
a = 0;
b = 20;
nel = 200; % number of elements
h = (b-a)/nel; % step size
nv = nel+1;% number of vertices
x = a:h:b; % mesh
dt = 0.001; % time steps
tend = 20;
e=ones(nv, 1);
%Discretization of the Laplace operator
Delta=1/(h^2)*spdiags([e -2*e e], -1:1, nv, nv);
Delta(1,2) = 2/h^2;
Delta(end,end-1) = -2/h^2;
% parameters of the model
eps = 0.02;
alpha = 0.1;
% initial guess
u = 0.5 * (x < 3);
v = 0 * x;
counter = 1;
% explicit Euler scheme
for t=dt:dt:tend
u = u + dt/eps*(u .*(1-u).*(u-alpha) - v - ((Delta * u')'));
v = v + dt * (u - v);
end
figcure(1)
plot(x,u,x,v,'r');
legend('u','v')
title('Solution of the Fitz-Hugh-Nagumo system')

Respuesta aceptada

Stephen23
Stephen23 el 27 de Nov. de 2018
Editada: Stephen23 el 27 de Nov. de 2018
Actually this has nothing to do with division by zero. When you read the NaN help it lists several ways that NaNs can occur:
"These operations produce NaN:"
  • "Any arithmetic operation on a NaN, such as sqrt(NaN)"
  • "Addition or subtraction, such as magnitude subtraction of infinities as (+Inf)+(-Inf)"
  • "Multiplication, such as 0*Inf"
  • "Division, such as 0/0 and Inf/Inf"
  • "Remainder, such as rem(x,y) where y is zero or x is infinity"
In your case your code has increasingly large values in both u and v until on the eighth iteration they include several Inf elements (which critically have the same locations in both arrays):
>> any(any(isinf(u)==isinf(v))) % colocated Inf values.
ans = 1
>> any(isnan(u(:)))
ans = 0
>> any(isnan(v(:)))
ans = 0
Thus on the ninth iteration you calculate this: v - ((Delta * u')'), which is thus a subtraction of infinities and (exactly as the documentation states) produces NaNs:
>> tmp = v - (Delta * u')';
>> any(isnan(tmp(:)))
ans = 1

Más respuestas (0)

Categorías

Más información sobre Loops and Conditional Statements en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by