Jumps in plotted data

Hi! So i know that the plotted graph should be a smooth curve, but no matter what I do, the plot does this on a small scale:

Any help appriciated! I have tried taking every nth value and interpolating but wierdly this make the data look even worse. I',m aware that the noise is likely due to my way of making the datapoints unique for interp1, but if I decrease the noise any further, interp1 says the sample points are no longer unique.

In case it is important, please also note that the sliv value in question is going to be an input variable for another function.

At this point, I've spent to long on this and have so much more left to do that if a proper solution can't be found, then I'm happy to just smooth the curve if there is some way to do so?

Here is my code currently, with the best outcome i have so far.:

clear all;

% Material Paramaters
% Thermal conductivity
ks      = 27.3; % solid W/m/K
kl      = 32.8; % liquid W/m/K
km      = 100;  % mould W/m/K

% Density
rhos   = 8110; % solid [kg/m^3]
rhol   = 8058; % liquid [kg/m^3]
rhom   = 2200; % mould [kg/m^3]

% Specific heat capacity
cps        = 711;  % solid [J/Kg/K]
cpl        = 700;  % liquid [J/Kg/K]
cpm        = 1700; % mould [J/Kg/K]

% Thermal Diffisivity
alphal = kl/(rhol*cpl); % Thermal diffusivity
alphas = ks/(rhos*cps); % Thermal diffusivity

% Solidification Range
T_solidus    = 1230; % [celcius]
T_liquidus   = 1270; % [celcius] 

% Process Paramaters
L = 2.5; % speration distance between hot and cold fields (cm)
mzt = 1350; % Melt zone temp in Celcius
czt = 900; % Cold zone temp in Celcius

% Spatial Discretisation
nx = 250; % Number of Spatial Nodes (cm)
dx = L/nx; % Spatial Node Length
x = 1:1:nx+1; % Initialise Position Vector

% Temporal Discretisation
nt = 100000; % Number of Temporal Nodes (deciseconds)
t = 100000; % Total Simulation Time
dt = t/nt; % Temporal Node Length
t = 1:1:nt; % Initialise Temporal Vector

% Initialise Temperature Vector
T = ones(1,nx+1); % Initialise Vector
T(1) = czt; % Inital Temp at cold zone
T(2:nx) = mzt; % Initial Temp in Liquid
T1 = T; % Initialise Update Vector

% Initialise S/L interface position (slip) and velocity (sliv) vectors
slip = zeros(1, nt);
sliv = zeros(1, nt);

% Solver Controls
F0s = alphas*dt/dx^2; % Solid Constant
F0l = alphal*dt/dx^2; % Liquid Constant

% 2D Grpahing Controls
num = 1;
check = round((logspace(1, log10(nt), 10)), 2, 'significant' ); % Takes logarithmically spaced sample points for evenly spaced 2D-graph plotting

% Solver
for j=2:nt % Temporal Loop
    for i=2:nx % Spatial Loop

        if T1(i) <= T_solidus       % Properties of Mould
            F0 = F0s;
        else                        % Properties of liquid
            F0 = F0l; 
        end

         T1(i) = F0*T(i+1) + (1-2*F0)*T(i) + F0*T(i-1); 
    end

    % Boundary Conditions
    T1(1) = czt;
    T1(nx) = mzt;

    %Temperature Update
    T=T1;
   
    % Grpahing
    num = num + 1;
        if ismember(num, check) == 1

            % 2D plot
            plot(x(1:nx),T1(1:nx)); hold on
            set(gca, 'YDir','reverse')
            xlabel('Position (cm)', "FontSize", 20);
            ylabel('Temperature (C)', "FontSize", 20);
            title('Temperature Gradient across mold and melt', "FontSize", 20);
        end

        % make slip vector
        slip_t = T + rand(1, nx+1)*1e-4;

        slip_idx = interp1(slip_t, 1:numel(slip_t), T_solidus);

        slip(j) = slip_idx;

end

Xnz = slip(slip ~= 0);                             % Non-Zero Elements
vnz = find(slip ~= 0);                             % Indices Of Non-Zero Elements
iv = 1:length(slip);                               % Interpolation Vector
slip = interp1(vnz, Xnz, iv, 'linear');            % Interpolate To Get Desired Output

hold off % seperating the plots

%Calculate sliv (S/L Interface Velocity)
for s = 2:1:nt
    sliv(s) = (slip(s)-slip(s-1))/dt;
end

% 2-D plot of S/L interface position approximation with time
figure;
plot(t, slip)
xlabel('Time (s)');
ylabel('S/L interface position (cm/s)');
title('S/L interface position evolution', "FontSize", 20);

% 2-D plot of S/L interface velocity approximation with time
figure;
plot(t, sliv)
xlabel('Time (s)');
ylabel('S/L interface Velocity (cm/s)');
title('S/L interface velocity evolution', "FontSize", 20);

[Mathieu NOE]

hello

I tried quickly a few things - at least here below you have now a code that works whataever the noise amplitude you give (0 is fine too) - but even w/o the noise the peaks remain and now the next question is where do those peaks come from and what can we learn from their x spacing - which does not sems to match with some initilization parameters - need more time to find out ...

to see my modifications , look for comments with % MN

clear all;
% Material Paramaters
% Thermal conductivity
ks      = 27.3; % solid W/m/K
kl      = 32.8; % liquid W/m/K
km      = 100;  % mould W/m/K
% Density
rhos   = 8110; % solid [kg/m^3]
rhol   = 8058; % liquid [kg/m^3]
rhom   = 2200; % mould [kg/m^3]
% Specific heat capacity
cps        = 711;  % solid [J/Kg/K]
cpl        = 700;  % liquid [J/Kg/K]
cpm        = 1700; % mould [J/Kg/K]
% Thermal Diffisivity
alphal = kl/(rhol*cpl); % Thermal diffusivity
alphas = ks/(rhos*cps); % Thermal diffusivity
% Solidification Range
T_solidus    = 1230; % [celcius]
T_liquidus   = 1270; % [celcius] 
% Process Paramaters
L = 2.5; % speration distance between hot and cold fields (cm)
mzt = 1350; % Melt zone temp in Celcius
czt = 900; % Cold zone temp in Celcius
% Spatial Discretisation
nx = 250; % Number of Spatial Nodes (cm)
dx = L/nx; % Spatial Node Length
x = 1:1:nx+1; % Initialise Position Vector
% Temporal Discretisation
nt = 100000; % Number of Temporal Nodes (deciseconds)
t = 100000; % Total Simulation Time
dt = t/nt; % Temporal Node Length
t = 1:1:nt; % Initialise Temporal Vector
% Initialise Temperature Vector
T = ones(1,nx+1); % Initialise Vector
T(1) = czt; % Inital Temp at cold zone
T(2:nx) = mzt; % Initial Temp in Liquid
T1 = T; % Initialise Update Vector
% Initialise S/L interface position (slip) and velocity (sliv) vectors
slip = zeros(1, nt);
sliv = zeros(1, nt);
% Solver Controls
F0s = alphas*dt/dx^2; % Solid Constant
F0l = alphal*dt/dx^2; % Liquid Constant
% 2D Grpahing Controls
num = 1;
check = round((logspace(1, log10(nt), 10)), 2, 'significant' ); % Takes logarithmically spaced sample points for evenly spaced 2D-graph plotting
% Solver
for j=2:nt % Temporal Loop
    for i=2:nx % Spatial Loop
        if T1(i) <= T_solidus       % Properties of Mould
            F0 = F0s;
        else                        % Properties of liquid
            F0 = F0l; 
        end
         T1(i) = F0*T(i+1) + (1-2*F0)*T(i) + F0*T(i-1); 
    end
    % Boundary Conditions
    T1(1) = czt;
    T1(nx) = mzt;
%     T1(nx+1) = mzt; % MN
    %Temperature Update
    T=T1;
   
    % Grpahing
    num = num + 1;
        if ismember(num, check) == 1
            % 2D plot
            plot(x(1:nx),T1(1:nx)); hold on
            set(gca, 'YDir','reverse')
            xlabel('Position (cm)', "FontSize", 20);
            ylabel('Temperature (C)', "FontSize", 20);
            title('Temperature Gradient across mold and melt', "FontSize", 20);
        end
        % make slip vector
        slip_t = T + rand(1, nx+1)*0;
%         slip_idx = interp1(slip_t, 1:numel(slip_t), T_solidus);
        slip_idx = find_zc(1:numel(slip_t),slip_t,T_solidus);  % MN
        slip(j) = slip_idx;
end

% MN - why all this below ?? 
% Xnz = slip(slip ~= 0);                             % Non-Zero Elements
% vnz = find(slip ~= 0);                             % Indices Of Non-Zero Elements
% iv = 1:length(slip);                               % Interpolation Vector
% slip = interp1(vnz, Xnz, iv, 'linear');            % Interpolate To Get Desired Output
% slip = [Xnz(1) Xnz];  % MN

hold off % seperating the plots

%Calculate sliv (S/L Interface Velocity)
for s = 2:1:nt
    sliv(s) = (slip(s)-slip(s-1))/dt;
end

sliv2 = dt*gradient(slip); % MN - why not simply this instead of the for loop above ?? 

% 2-D plot of S/L interface position approximation with time
figure;
plot(t, slip)
xlabel('Time (s)');
ylabel('S/L interface position (cm/s)');
title('S/L interface position evolution', "FontSize", 20);

% 2-D plot of S/L interface velocity approximation with time
figure;
semilogy(t, sliv, t, sliv2)
xlabel('Time (s)');
ylabel('S/L interface Velocity (cm/s)');
title('S/L interface velocity evolution', "FontSize", 20);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Zx] = find_zc(x,y,threshold)

% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
Zx = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end

[dpb]

@Mathieu NOE, I didn't have time to delve, but

T = ones(1,nx+1); % Initialise Vector

catches my eye because it doesn't look as though that extra last element gets updated in the spatial loop? My hunch is it's got to do with that...

[Mathieu NOE]

yep that's what I also tested but it doesn't seem to have any impact on the peaks in the last plot .

still I think it should be done this way

[Charlie Rideout]

Thanks for having a look at least. Is there a particular reason for defining a new function instead of using interp1? I can't really follow what's going on in your code haha

[Mathieu NOE]

hello again

I wanted to try another approch especially when you said that interp1 would fail here when the noise amplitude is too small - because the data is not anymore monotonic and unique

you can see now that the new function avoids this situation , because IMHO it's a better way to find why index of slip_t correspond to slip_t = T_solidus.

there is nothing wrong with the interp1 method but it fails when data is not unique as you experienced -- so you must be aware of this limitation of correct the data accordingly (if you can)

%         slip_idx = interp1(slip_t, 1:numel(slip_t), T_solidus);
slip_idx = find_zc(1:numel(slip_t),slip_t,T_solidus);  % MN

[dpb]

Again, I don't have sufficient time to spare to dig really deeply, but here's a way to see what's going on...I ran the code with MN's updates and then did

dS=diff(slip);
plot(t,slip)
xlim([1 100])
ylabel('Slip')
vyyaxis right
plot(t,dS,'r-')
ylim([0 0.2])
xlabel('Time')
ylabel('delta Slip')

This amplified scale illustrates that there is indeed a fluctuation in the Slip trace even though the previous plots look completely smooth because of the scale encompassing the entire range. Over a zoomed-in range, one can see the discontinuity that is introduced every so often; as noted earlier it is not a fixed time interval between the discontinuities which does indicate it isn't associated with the time vector reset but within the restart of the state change values.

There is a very large discontinuity at the outset; the initial Slip value is zero which jumps to 1.7786 at the first timestep; hence the first delta is huge in comparison to the rest being that same magnitude where as it then is uniformly under 0.2 from there on. This is a problem in the initialization to be corrected, but won't cure the other discontinuity.

I then also did the second derivative and poked around at it some...

K>> dS2=diff(dS);                   % now look at second difference
K>> plot(t(3:end),dS2)
K>> xlim([0 1E2])
K>> ylim(0.015*[-1 1])
K>> xlabel('Time')
K>> ylabel('delta(delta Slip)')     % that's interesting...let's find the peaks we see locations
K>> ix=find(dS2>0.001);             % just arbitrary for now by inspection of plot
K>> dix=diff(ix);                   % and how far apart are they?
K>> dix=dix(dix>1);                 % eliminate the one that returned two points, not just one
K>> diff(dix)                       % and now the second difference here
ans =
    15    15    16    15    15
K>>

which is most interesting, indeed! This shows the differential in the discontinuities is increasing by 15/16 points each time, at least in the early part of the trace. Should be able to find the cause of that constant increasing difference between the jumps.

Hopefully, this will help to uncover where the problem(s) lie...

[Mathieu NOE]

I just played again with the code and found something intriguing... I naively thought that increasing the number of spatial points (nx) , the results would improve, get's smoother

so I plotted T for a 1000 time iteration duration (to make it faster) and was surprised to find some strange behaviour

up to nx = 500 the T profiles are always nice looking and smooth like :

the horizontal red line correspond to T_solidus (and the intersection with the T profiles are supposed to give us slip_idx)

(plots are for time index : 2 10 17 28 46 77 130 220 360 600 1000)

if I increase nx at 750 I get this :

and it gets even worse if you continue to increase nx

so I wonder if there is a numerical issue in the spatial for loop

code a bit modified to look at this phenomenon :

clear all;
close all

% Material Paramaters
% Thermal conductivity
ks      = 27.3; % solid W/m/K
kl      = 32.8; % liquid W/m/K
km      = 100;  % mould W/m/K
% Density
rhos   = 8110; % solid [kg/m^3]
rhol   = 8058; % liquid [kg/m^3]
rhom   = 2200; % mould [kg/m^3]
% Specific heat capacity
cps        = 711;  % solid [J/Kg/K]
cpl        = 700;  % liquid [J/Kg/K]
cpm        = 1700; % mould [J/Kg/K]
% Thermal Diffisivity
alphal = kl/(rhol*cpl); % Thermal diffusivity
alphas = ks/(rhos*cps); % Thermal diffusivity
% Solidification Range
T_solidus    = 1230; % [celcius]
T_liquidus   = 1270; % [celcius] 
% Process Paramaters
L = 2.5; % speration distance between hot and cold fields (cm)
mzt = 1350; % Melt zone temp in Celcius
czt = 900; % Cold zone temp in Celcius
% Spatial Discretisation
nx = 50*10; % Number of Spatial Nodes (cm) (250) % seems above 500 the results (T) becomes unstable
dx = L/nx; % Spatial Node Length
x = 1:1:nx+1; % Initialise Position Vector
% Temporal Discretisation
nt = 1000; % Number of Temporal Nodes (deciseconds)
t = nt; % Total Simulation Time
dt = t/nt; % Temporal Node Length
t = 1:dt:t; % Initialise Temporal Vector
% Initialise Temperature Vector
T = ones(1,nx+1); % Initialise Vector
T(1) = czt; % Inital Temp at cold zone
T(2:nx) = mzt; % Initial Temp in Liquid
T(nx+1) = mzt; % MN

T1 = T; % Initialise Update Vector
% Initialise S/L interface position (slip) and velocity (sliv) vectors
slip = zeros(1, nt);
sliv = zeros(1, nt);
% Solver Controls
F0s = alphas*dt/dx^2; % Solid Constant
F0l = alphal*dt/dx^2; % Liquid Constant
% 2D Grpahing Controls
num = 1;
check = round((logspace(1, log10(nt), 10)), 2, 'significant' ); % Takes logarithmically spaced sample points for evenly spaced 2D-graph plotting

check = [2 check]; % MN add first time iteration

% Solver
    figure(1) % MN
    hold on   % MN
    
for j=2:nt % Temporal Loop

    for i=2:nx % Spatial Loop
        if T1(i) <= T_solidus       % Properties of Mould
            F0 = F0s;
        else                        % Properties of liquid
            F0 = F0l; 
        end
         T1(i) = F0*T(i+1) + (1-2*F0)*T(i) + F0*T(i-1); 
    end
    % Boundary Conditions
    T1(1) = czt;
    T1(nx) = mzt;
    T1(nx+1) = mzt; % MN
    %Temperature Update
    T=T1;
   
    % Graphing
    num = num + 1;
%         if ismember(num, check) == 1
%             % 2D plot
%             plot(x(1:nx),T1(1:nx)); hold on
%             set(gca, 'YDir','reverse')
%             xlabel('Position (cm)', "FontSize", 20);
%             ylabel('Temperature (C)', "FontSize", 20);
%             title('Temperature Gradient across mold and melt', "FontSize", 20);
%         end
%         % make slip vector
%         slip_t = T + rand(1, nx+1)*0;
% %         slip_idx = interp1(slip_t, 1:numel(slip_t), T_solidus);
%         slip_idx = find_zc(1:numel(slip_t),slip_t,T_solidus);  % MN
%         slip(j) = slip_idx;
        
        if ismember(num, check) == 1
        plot(T(1:15),'*-')   % MN
        end
        
end
plot(T_solidus*ones(1,15),'r--');  % MN
title(['T plot vs time index for nx = ' num2str(nx)])
xlabel('Position')
ylabel('T (°C)')

[dpb]

I don't know what is trying to be modelled well enough to derive the actual physical system equations, but I think the spatial iteration definitely has issues--

for i=2:nx % Spatial Loop
    if T1(i) <= T_solidus       % Properties of Mould
        F0 = F0s;
    else                        % Properties of liquid
        F0 = F0l; 
    end
     T1(i) = F0*T(i+1) + (1-2*F0)*T(i) + F0*T(i-1)

What are the units of F? Unless it is unitless, then T has units of degrees*units(F), not degrees.

And, I don't follow the (1-2*F0) factor weighting T(i)?

And the second, F0 is based on T1(i), but used for all i-1, i+1 locations that could be across the boundary.

I think this needs a set of diff EQs to be solved by numerical integration or a lot more detail in linearizing and finite differencing a numeric solution.

But, as noted, I've no idea what the problem actually is other than something to do about a liquid/solidus interface...

ADDENDUM:

OK, I checked the units and the alpha are in m^2/s so multiplying by dt/dx^2 is indeed unitless. So, if one assumes unity in y, then dx^2 would be proportional to an area.

However, what appears to be totally missing in the model is the latent heat of the phase change; endothermic for solid-liquid, exothermic the other way (for ordinary materials, anyway; don't see that the material is identified here). When the material changes phase, there's a period with zero temperature change until the latent heat is absorbed; meanwhile in the other direction there is exothermic release of additonal energy/heat without a bulk temperature change until the material state change is complete. These will need to be addressed and ignoring it may possibly be the root cause of the spikes although I've not delved more deeply.