Results are NaN or absurd if a certain parameter changes.

Hello everyone,

The code I have given below is the explicit 1-D finite difference solution of a heat transfer question. The problem I'm having is that depending on the number of nodes (n), the results differ significantly. Between 2 and 46 nodes, the results are reasonable where accuracy, if my solution is correct, improves. However, between nodes 47 and 50, the results are erratic and I can't make any sense of them. After 50, results are given as NaN.

I'm asuuming that beyond 50, "dx" becomes too small, however I don't know what happens between 46 and 50. I can say that lowering "dx" results in a similar behaviour.

I was wondering if there was a problem with my poorly written code or am I missing something else?

clear, clc;
n = 6; %number of nodes 
L = 0.05; 
dx = L/(n-1); 
k = 28;
h = 60;
q = 6*10^5;
alp = 12.5*10^-6;
dt = .1;
t = 0:dt:100;
Ta = 30;
Ti = 200;
tau = alp*dt/dx^2; 
c = zeros(length(t),1);
for j = 1:length(t)
     if j == 1
        T(1) = tau*(q*dx^2/(2*k) + Ti) + Ti*(1-tau);  
        T(n) = tau*(Ti + q*dx^2/(2*k) + h*Ta*dx/k) + Ti*(1-tau-tau*h*dx/k);
     else
        T(1) = tau*(q*dx^2/(2*k) + Ti) + T(1)*(1-tau);  
        T(n) = tau*(T(n-1) + q*dx^2/(2*k) + h*Ta*dx/k) + T(n)*(1-tau-tau*h*dx/k);
     end
    for i = 2:1:n-1 
     if j == 1
        T(i) = tau*(Ti+Ti) + (1 - 2*tau)*Ti + (tau*q*dx^2)/(2*k);
     else
        T(i) = tau*(T(i-1)+T(i+1)) + (1 - 2*tau)*T(i) + (tau*q*dx^2)/(2*k); 
     end
    end
end 

m = (1:n);
mt = transpose(m);
Tt = transpose(T);
table(mt,Tt) 
figure
plot(mt,T)

[Jan]

Just a hint: Simplify

tau*(Ti+Ti) + (1 - 2*tau)*Ti + (tau * q * dx^2) / (2 * k)   % ???

to:

Ti + (tau * q * dx^2) / (2 * k)

[Furkan Erdem]

I wouldn't ever have noticed that, I should be more orgnaised. Unfortunately, in this instance, the code is for a class so I'm going to have to leave it.

Answer 1

[Tanmay Das]

Hi,

By default, all numeric values are stored as double data type. Once the number of nodes cross 50, the values in T array exceed the bounds of double data type. MATLAB constructs the double data type according to IEEE® Standard 754 for double precision. The range for a negative number of type double is between -1.79769 x 10^308 and -2.22507 x 10^-308, and the range for positive numbers is between 2.22507 x 10^-308 and 1.79769 x 10^308. You can look into the double Documentation for further information.

Answer 2

[darova]

I suggest you to visualize the data in 3D. See how values change each step

T = zeros(20);
a = 0.04;                   % constant
t = linspace(0,pi,20);
T(1,:) = sin(t);            % initial condition
h = surf(T);
for i = 1:19
    for j = 2:19
         T(i+1,j) = T(i,j) - a*(T(i,j+1)-T(i,j-1)+2*T(i,j));
    end
    set(h,'zdata',T)
    pause(0.5)
end