Temperature matrix giving NaN

clc
clear
close all



%% initialisation
alpha1=0.02;
a=0.6;
b=0.05;
c=0.01;
q=0.05;
g=20;
L=4;
tmax=100;
dx=0.01;
nx=L/dx+1;
dt=0.005;
nt=tmax/dt+1;
x=0:dx:L;
t=0:dt:tmax;
ap=1/dt+2*alpha1/dx^2;
T=zeros(nx,nt);
f=zeros(nx,1);
v=zeros(nt,1);
S=zeros(nx,1);

That is my code. I want to know what is wrong with my intialisation because my Temperature matrix keeps giving me NaN.

[Torsten]

What's the PDE you are trying to solve with your code ?

[Izk]

The PDE is above. This is the task.

Answer 1

[Sam Chak]

@Mohammad Hassan

The T matrix gives NaN because the temperature T is an array of zeros to begin with

T = zeros(nx, nt);  

and f is a scalar that returns a zero value at the end of the loop

f = T(i, k-1)/dt + S(i);

Thus, causing the following to give some Inf and NaN.

T = A/f

T =

Inf Inf NaN -Inf NaN

-Inf Inf Inf NaN NaN

NaN -Inf Inf Inf NaN

NaN NaN -Inf Inf Inf

NaN Inf NaN -Inf Inf

[Izk]

So how can I change it so the matrix shows how the temperature changes? Do I get rid of the T = zero(nx, nt); ?

[Torsten]

You were told to set T0=g=20. So just do it:

T(:,1) = g;

Further,

f(i)=T(i,k-1)/dt+S(i);

instead of

f = T(i,k-1)/dt+S(i);

Further,

T(:,k) = A\f;

instead of

T=A/f;

and this before the last end-statement.

Didn't check in detail if the periodic boundary condition is correctly implemented.

[Izk]

So I did what you said, however got an error saying Unrecognized function or variable 'T'.

clc
clear
close all



%% initialisation
alpha1=0.02;
a=0.6;
b=0.05;
c=0.01;
q=0.05;
g=20;
L=4;
tmax=100;
dx=0.01;
nx=L/dx+1;
dt=0.005;
nt=tmax/dt+1;
x=0:dx:L;
t=0:dt:tmax;
ap=1/dt+2*alpha1/dx^2;
T=zeros(nx,nt);
T(:,1)=g;
f=zeros(nx,1);
v=zeros(nt,1);
S=zeros(nx,1);

T(:,k)=A\f;

Could the problem be with the line 'v=a+b*cos(0.5*pi*dt*(k-1))+c*cos(4*pi*dt*(k-1));' ? It says error using Evalin.

[Torsten]

  ...
  for i=1:nx
    S(i)=q*exp(-0.5*((dx*i-L/2)/0.05).^2);
    f(i)=T(i,k-1)/dt+S(i);
  end
  T(:,k)=A\f;
end

[Izk]

clc
clear
close all

%% initialisation
alpha1=0.02;
a=0.6;
b=0.05;
c=0.01;
q=0.05;
g=20;
L=4;
tmax=100;
dx=0.01;
nx=L/dx+1;
dt=0.005;
nt=tmax/dt+1;
x=0:dx:L;
t=0:dt:tmax;
ap=1/dt+2*alpha1/dx^2;
T=zeros(nx,nt);
T(:,1)=g;
f=zeros(nx,1);
v=zeros(nt,1);
S=zeros(nx,1);

The matrix disappears after this and I get message saying:

'NaN/Inf breakpoint hit for workspacefunc.m on line 310.

310 v = evalin('caller', varargin{1});

I'm not sure what part is causing this.

[Torsten]

Remove the line

v = zeros(nt,1)

These were the obvious errors in your code.

Now you will have to check whether the matrix A is correct and the boundary conditions are correctly implemented.

Answer 2

[Torsten]

%% initialisation
alpha1=0.02;
a=0.6;
b=0.05;
c=0.01;
q=0.05;
g=20;
L=4;
tmax=100;
dx=0.01;
nx=L/dx+1;
dt=0.005;
nt=tmax/dt+1;
x=(0:dx:L).';
t=0:dt:tmax;
T=zeros(nx,nt);
T(:,1)=g;
A=zeros(nx,nx);
f=zeros(nx,1);
%% main body
for k=2:nt
    v=a+b*cos(0.5*pi*t(k))+c*cos(4*pi*t(k));
    A = diag((1 + 2*alpha1*dt/dx^2)*ones(1,nx))+...
        diag((v*dt/(2*dx) - alpha1*dt/dx^2)*ones(1,nx-1),1)+...
        diag((-v*dt/(2*dx) - alpha1*dt/dx^2)*ones(1,nx-1),-1);
    A(1,nx-1) = -v*dt/(2*dx) - alpha1*dt/dx^2;
    A(nx,2) = v*dt/(2*dx) - alpha1*dt/dx^2;
    f(1:nx) = T(1:nx,k-1)+dt*q*exp(-0.5*((x(1:nx)-L/2)/0.05).^2);
    T(:,k) = A\f;
end
plot(x,[T(:,100),T(:,200),T(:,2000),T(:,10000)])
end