finite difference implicit method

Hi everyone, I have written this code but I do not know why Matlab does not read the if condition. It suppose to use different variable for (alfa) when it is reach N= 33, 66.

Could you please help me with that.

clear all
clc

ho=50;
hi=10;
N=99; %Number of space intervals
M=2000; %number of time steps
Th=196; %thickness of wall in meter;;

t_0=0; %initial time
t_f=3600; %final time (sec)
To=34; %exterior wall temperature
k_b=0.7; %thermal conductivity of bricks
roh_b=1600; %density of bricks (kg/m^3)
Cp_b=840; %Cp of bricks (J/kg-K)
T_in=22; %Inner wall temperature
abs=0.7;
dx = (Th/(N-1))*0.001; %size of the mesh
% dx=0.002;
x = 0:dx:N-1;
x = x'; %inverse of mesh
bi=hi*dx/k_b; %coefficient of terms
bo=ho*dx/k_b;

%n-Octadecane PCM
k_ps=0.358;
k_pl=0.148;
Cp_ps=1.934;
Cp_pl=2.196;
roh_ps=865;
roh_pl=780;
lamda=243.5;
dT=1;
Tm=25;

Gt=solarradiationcode;
%initial condition
u(:,1)=ones(N,1)*15;
for index=1:length(Gt)
% index=1;
for j=1:M  %time step
    V_Te=To*(abs*Gt/ho);
    Te=V_Te(index);

   dt =3;
  t = t_0:dt:t_f;

    alfa= k_b/(Cp_b*roh_b);
      % define the ratio r

      r=alfa*dt/dx^2; %Fourier number

      Beta=[2*r*bo*Te; zeros(N-2,1); 2*r*bi*T_in];

      A = zeros(N,N);
      A(1,1) = 1+2*r+2*r*bo;
      A(1,2) = -2*r;

      for i=2:N-1

          if N/3> i || i>2*N/3
              Cp=Cp_b;
              roh=roh_b;
              k=k_b;

              alfa= k/(Cp*roh);
              r=alfa*dt/dx^2;
          end

          if N/3==i && i==2*N/3

          Cp=2*Cp_b*Cp_pl/(Cp_pl+Cp_b);

              roh=2*roh_pl*roh_b/(roh_pl+roh_b);

              k=2*k_b*k_pl/(k_pl+k_b);

              alfa= k/(Cp*roh);
              r=alfa*dt/dx^2;
          end

          if N/3<i && i<2*N/3
              T=u(i,j);
              if T<Tm

                  Cp=Cp_ps;
                  roh=roh_ps;
                  k=k_ps;

              elseif Tm<T && T<Tm+dT

                  Cp= ((Cp_ps+Cp_pl)/2)+((roh_ps+roh_pl)/2)*(lamda/dT);

              elseif T>Tm+dT

                  Cp=Cp_pl;
                  roh=roh_pl;
                  k=k_pl;

              end
              alfa=k/(Cp*roh);
              r=alfa*dt/dx^2;
          end

          A(i,i-1) = -r;
          A(i,i) = 1+2*r;
          A(i,i+1) = -r;
      end

      alfa= k_b/(Cp_b*roh_b);
      % define the ratio r
      r=alfa*dt/dx^2; %Fourier number

      %A(N,N+1) = 0;
      A(N,N-1) = -2*r;
      A(N,N) = 1+2*r+2*r*bi;

      [L,U] = tridia(A);

      [v] = lower_solver(L,u(:,j)+Beta);
      u(:,j+1)= uper_solve(U,v);

end

end
plot(u(:,j+1))

Answer 1

[the cyclist]

We can't run your code, because you did not supply us with solarradiationcode. In which line do you expect there to be a branching?

My guess is that your problem is in one of your lines that is like

if N/3==i && i==2*N/3

You seem to be doing equality comparisons between floating-point numbers, which is not appropriate because not all numbers can be expressed exactly in floating-point arithmetic.

See, for example, some of the answers here for a fuller explanation.