Why is time steps changing the graphic? (diffusion equation)

This is a model to know how the heat will pass through the soil (diffusion equation). It's working good with this time steps, but my problem is when I put the dt bigger than 500 (even if I change n_steps). The graphic becomes crazy, it's totally unphysical. Why is it happening? Thank you.

% set model parameters & constants dt = 240; % time step (s)

n_steps = 2*15; % number of time steps (15 steps = 1h)

J = 10; % number of soil levels

D = 1e-7; % diffusion coefficient (m^2 s^{-1})

soil_depth = 0.1; % soil depth (m)

dz = soil_depth/J; % depth increment (m)

DD = D*dt/(dz*dz); % calculate coefficient for FTCS scheme

% set up vertical z grid (j index) and time grid (n index)

ja = [1:J]; % set up an array for vertical levels

ja2 = [2:J-1]; % set up an array for the interior vertical levels

z(ja) = (ja - 0.5)*dz; % define z grid with levels at midpoint of soil layers

na = [1:n_steps]; % set up an array for time grid

% define initial conditions, T_i is initial temperature, i.e. T_i(ja)=T(0,ja)

Trock = 10;

T_i(ja) = Trock;

Tair= 40;

% The function Tair was defined in other script and called

% function t = Tair(n)

% t=10+10*sin(2*pi*n/360) obs: dt=240 or 4 minutes, so 1/360 of day

% end

% Problem becomes calculation of T(n,j) , where n index is for time, j index is for depth

T(na,ja) = NaN %set up array of correct dimensions

% Calculate first time step ----------------------------------------

% Note this has to be done separately, as Matlab indices must be integers

% Use top boundary condition derived in lectures (section 6.2)

T(1,1) = T_i(1)+DD*(T_i(2)-3*T_i(1)+2*Tair(1)); % Top BC

T(1,ja2) = T_i(ja2)+DD*(T_i(ja2+1)-2*T_i(ja2)+T_i(ja2-1)); % Interior points

T(1,J) = T_i(J)+DD*(2*Trock-3*T_i(J)+T_i(J-1)); % Bottom BC

% Calculate all subsequent time steps----------------------------------

% Use for loop over time (n), filling in T array as we go

for n=2:n_steps

T(n,1) = T(n-1,1)+DD*(T(n-1,2)-3*T(n-1,1)+2*Tair);

T(n,ja2) = T(n-1,ja2)+DD*(T(n-1,ja2+1)-2*T(n-1,ja2)+T(n-1,ja2-1));

T(n,J) = T(n-1,J)+DD*(2*Trock-3*T(n-1,J)+T(n-1,J));

end

% plot results-----------------------------------------

plot(T(30,:),z,'b-o')

hold on;

xlabel('Temperature (^o C)')

ylabel('Depth (m)')

set(gca,'ydir','reverse')

axis([0 40 0 0.1])

Answer 1

[Mohammad Abouali]

haven't read through your code, but you change your time step it is highly possible that your scheme becoming unstable and it is just producing garbage.

There is a limit for time step for each numerical scheme. If you go beyond that your solution becomes unstable and as you mentioned unphysical results.

[Rodrigo]

But is there any common reason for it?

Thank you.

[Mohammad Abouali]

Yes, it depends on how the error is amplified. That depends on the numerical scheme that was used to discretize the equation in space and the method that is used to discretize the equation in time (or integrate in time).

Also being an explicit or implicit scheme plays a major role.

Check "Stability Analysis" and you would find many article on the subject. One of the easiest method to be used is Von Neumann Stability Analysis. Although that requires periodic boundary condition. The other approach is to look at the eigen values of the discretized system and see if it is going to amplify the error or damp the error.

It is indeed a big topic, this stability analysis. As I said, I have not checked your code, but the sign that you mentioned are the classic sign of instability in solving PDEs numerically.

[Mohammad Abouali]

If you are solving the following equation:

dU/dt=alpha * d^2U/dx^2

then Von Nuemann Stability analysis using FTCS discretization gives you the following condition

k<=1/2 where K is the diffusion number and is defined as

k= alpha * dt / dx^2.

This means that dt and dx should be chosen in such a way that k stays less than 0.5 otherwise, as the simulation time progresses, the error would grow to finally take over the solution.

As you can see two things can cause your scheme become unstable, (1) increasing dt (or time step) or (2) decreasing spacing, i.e. dx or refining your grid in space.

It appears the diffusion number in your code is DD i.e. DD=D*dt/(dz*dz) and D is alpha in above equation. Now:

D=1e-7 dz=0.1/10=0.01

Now if dt=240 you get

DD=1e-7*240/0.01^2=0.25<= 0.5 So it should be stable

but dt=500

DD=1e-7*240/0.01^2=0.5<=0.5 so anything more than 500 shouldn't be stable, (dt=500 would have too much error too)

Note the K<=0.5 is for periodic boundary condition. If you are using another boundary condition the limit might be different.