Unable to meet integration tolerances without reducing the step size

Nx = 50;         % Number of grid points
L = 5000;        % Length of the spatial domain
dx = L / (Nx-1); % Spatial step size
x = linspace(0, L, Nx); % Spatial points
A=0.1963;
f = 0.008;
D = 500e-3;
rho=1.2;
c = 348.5;
P_atm = 1.013e5; % Atmospheric pressure in Pa
tspan = [0, 3600]; % Time span

% Initial conditions
P0 = 5e6*ones(Nx,1); % Initial condition for pressure (P)
m0 = zeros(Nx, 1); % Initial condition for mass flow rate (m)

% Initial conditions vector
u0 = [P0; m0];
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);

% Solve the system of ODEs using ode15s
[t, u] = ode15s(@(t, u) MyODEs(t, u, Nx, dx, f, c, D), tspan, u0,options);

Warning: Failure at t=8.024469e+02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.818989e-12) at time t.

% Extract the solutions for P and m
P_solution = u(:, 1:Nx);
m_solution = u(:, Nx + 1:end);

%% when m =kg/s*m^2
Q_n= (m_solution./rho); 
Q=(Q_n .*A);

% nodes_to_plot = 1:7:Nx;
% 
% % Plotting P_solution for selected nodes with transparency
% figure;
% hold on;
% colors = rand(length(nodes_to_plot), 3);
% for idx = 1:length(nodes_to_plot)
%     i = nodes_to_plot(idx);
%     plot(t, P_solution(:, i), 'LineWidth', 0.5, 'DisplayName', ['Node ' num2str(i)], 'Color', colors(idx,:), 'LineStyle', '-');
% end
% hold off;
% xlabel('Time (t)');
% ylabel('Pressure (P)');
% title('Pressure (P) at Selected Nodes over Time');
% legend('show');

% % Choose Node 1 (x = 1) for plotting
 node1_idx = 25;

% % Extract the pressure and Q values at Node 1
Q_node1 = Q(:, node1_idx);
figure;
plot(t,  Q_node1 , 'LineWidth', 2);
xlabel('Time (t)');
ylabel('Q at Node 1');
title('Q at Node 1 over Time with covective term');

% Set up the system of ODEs including the boundary condition ODEs
function dudt = MyODEs(t, u, Nx, dx, f, c, D)
    P = u(1:Nx);
    m = u(Nx + 1:end);
    dPdt = zeros(Nx, 1);
    dmdt = zeros(Nx, 1);

    % Boundary Condition
    P(1) = 5e6; % P(1) boundary condition
    
    % Nozzle at i=25 (Algebraic relationships)
    p_hst = 503.51e3; % p_hst = rho * g * d (where d = 40m)
    Cd = 0.55;
    gamma = 1.4;
    rho1 = 1.2;
   
  % Mass flowrate of nozzle 
   m_nz = Cd * sqrt((2 * gamma / (gamma - 1)) * P(25) * rho1 * (1 - (P(25) / p_hst)^((gamma - 1) / gamma)) * (P(25) / p_hst)^(2 / gamma)); 

  % Equation
   dPdt(25) =-c^2 * (3*m(25) - 4*m(24) + m(23)) / (2*dx); % Backward discretization of dP/dt at i=25
   dmdt(25) =-(2*m(25)*c^2/P(25)*((- 3*m(25)+4*m(26)-m(27)) / (2 * dx)))+(c^2*m(25)^2/P(25)^2)*((- 3*P(25)+ 4*P(26)-P(27)) / (2*dx)) - ((-3*P(25)+4*P(26)-P(27)) / (2*dx)) - (f * m(25) * abs(m(25)) * c^2) / (2 *D* (P(25))); % Forward discretization for dm/dt at i=26
   dmdt(1) =-(2*m(1)*c^2/P(1)*((- 3*m(1)+4*m(2)-m(3)) / (2 * dx)))+(c^2*m(1)^2/P(1)^2)*((- 3*P(1)+ 4*P(2)-P(3)) / (2*dx)) - ((-3*P(1)+4*P(2)-P(3)) / (2*dx)) - (f * m(1) * abs(m(1)) * c^2) / (2 *D* (P(1))); % Forward discretization for dm/dt at i=1
    
   % Update m(Nx) based on the time t and the specified conditions
    for i = 1:length(t)
        if t(i) < 600
            m(Nx) = 0;
        elseif t(i) < 1800
            m(Nx) = 509.28; %509.28 kg/s*m^2;
       else 
            m(Nx) = 0;
        end
    end
   
    % Equations
    dPdt(Nx) =-c^2 * (3*m(Nx) - 4*m(Nx - 1) + m(Nx-2)) / (2*dx); % Backward discretization of dP/dt at i=Nx
    
 % Central discretization for in-between grid points (i=2 to 49)
for i = 2:49
    if i == 25
        % Skip the calculation for dPdt(25) since it's already calculated outside the loop
        % dPdt(25) = ...; 
        % dmdt(25) =.....;
    elseif i == 24
         m(i)= m(25)+m_nz;% Two mass flow rate enter m(25) and m_nz 
         
    else
        % Calculate dP/dt and dm/dt for in-between grid points using central discretization
        dPdt(i) = -c^2 * (m(i+1) - m(i-1)) / (2 * dx); % Central discretization of dP/dt
        dmdt(i) =-2*m(i)*c^2/P(i)*((m(i + 1) - m(i - 1)) / (2 * dx))+(c^2*m(i)^2/P(i)^2)*(P(i+1) - P(i-1)) / (2*dx) - (P(i+1) - P(i-1)) / (2*dx) - (f * m(i) * abs(m(i))* c^2 ) /( 2*D*P(i)); % Central discretization of dm/dt

    end
end

    % Combine dPdt and dmdt into a single vector
    dudt = [dPdt; dmdt];
end

%% Due to this condition ,i am getting warning that integration tolerance cannot be meet 
% The logic mass flowrate entering at node 24 , is summation of m(25)+ m_nz 
%elseif i == 24
%m(i)= m(25)+m_nz;% Two mass flow rate enter m(25) and m_nz 

Answer 1

[Walter Roberson]

[t, u] = ode15s(@(t, u) MyODEs(t, u, Nx, dx, f, c, D), tspan, u0,options);

ode15s is defined as passing a scalar value as the first parameter -- so t will be scalar.

    for i = 1:length(t)

length(t) is guaranteed to be 1 -- unless, that is, you are somehow calling MyODEs from your own custom driver.

        if t(i) < 600
            m(Nx) = 0;
        elseif t(i) < 1800
            m(Nx) = 509.28; %509.28 kg/s*m^2;
       else 
            m(Nx) = 0;
        end

All of the ode*() variable-step routines are based upon mathematics that requires that the ode function have continuous first two derivatives. But when you have that kind of if statement, you have an abrupt jump in your m(Nx) value, leading to discontinuity in your first derivatives. When you are fortunate, ode15s will notice the discontunity and give you an error message about failing on integration tolerances. When you are not fortunately, ode15s will fail to notice and will keep going and confidently give you the wrong answer.

Your transitions are at specific times. You need to divide your integration up into three parts: 0 to 600, 600 to 1800, then over 1800. You will need three different ode*() calls with three different tspan for this. You would parameterize the call to myODEs according to what the current m(Nx) value should be for that tspan.

I notice by the way that you loop over i = 1:length(t) but your assignments to m are not indexed by i or an expression dependent on i: they are indexed at the constant value Nx . Nx does not change inside the loop, so the value of m(Nx) is going to come out according to whatever was assigned for the last t value, t(end) . There is no point in checking the other t(i) values: you would get the same effect if you had for i = length(t) instead. You should rethink that logic... especially in view of the fact that as discussed above, t will only ever be scalar unless you are calling the function "manually".

[MAYUR MARUTI DHIRDE]

The m(Nx) is a time-varying boundary condition, it is implemented correctly. When I comment else if i==24, m(i)=m(25)+m_nz , the code runs and matches the result.

%% LOGIC WHICH I WANT TO INTRODUCE IS EXPLAINED BELOW.

First, we have two boundary conditions P(1) = 5e6 Pas and m(Nx). Next, we have a total of 50 Nodal points, out of which we have specified the pressure boundary condition as P(1) at the first nodal point and the mass flow rate boundary condition as m(Nx), at the last grid point. Now we need to have an equation to calculate dmdt(1) at the first grid point using forward discretization and an equation for dPdt(Nx) at the last grid point using Backward discretization.

At node point 25 we have a nozzle, but since the nozzle does not account for any space it is calculated using algebraic equations m_nz

we calculate dPdt(25) using backward discretization and dmdt(25) using forward discretization.

Now at node point i=24 , the mass entering will be m(25) +m_nz

To calculate in between grid points we use central difference discretization and we skip the calculation of i=25 from the loop because it is calculated outside of the loop.

Nx = 50;         % Number of grid points
L = 5000;        % Length of the spatial domain
dx = L / (Nx-1); % Spatial step size
x = linspace(0, L, Nx); % Spatial points
A=0.1963;
f = 0.008;
D = 500e-3;
rho=1.2;
c = 348.5;
P_atm = 1.013e5; % Atmospheric pressure in Pa
tspan = [0, 3600]; % Time span

% Initial conditions
P0 = 5e6*ones(Nx,1); % Initial condition for pressure (P)
m0 = zeros(Nx, 1); % Initial condition for mass flow rate (m)

% Initial conditions vector
u0 = [P0; m0];
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);

% Solve the system of ODEs using ode15s
[t, u] = ode15s(@(t, u) MyODEs(t, u, Nx, dx, f, c, D), tspan, u0,options);

% Extract the solutions for P and m
P_solution = u(:, 1:Nx);
m_solution = u(:, Nx + 1:end);

%% when m =kg/s*m^2
Q_n= (m_solution./rho); 
Q=(Q_n .*A);

% % Choose Node 1 (x = 1) for plotting
 node1_idx = 1;

% % Extract the pressure and Q values at Node 1
Q_node1 = Q(:, node1_idx);
figure;
plot(t,  Q_node1 , 'LineWidth', 2);
xlabel('Time (t)');
ylabel('Q at Node 1');
title('Q at Node 1 over Time with covective term');

% Set up the system of ODEs including the boundary condition ODEs
% Set up the system of ODEs including the boundary condition ODEs
function dudt = MyODEs(t, u, Nx, dx, f, c, D)
    P = u(1:Nx);
    m = u(Nx + 1:end);
    dPdt = zeros(Nx, 1);
    dmdt = zeros(Nx, 1);

    % Boundary Condition
    P(1) = 5e6; % Define P(1) boundary condition
    % Nozzle at i=25 (Algebraic relationships)
    p_hst = 503.51e3; % p_hst = rho * g * d (where d = 40m)
    Cd = 0.55;
    gamma = 1.4;
    rho1 = 1.2;

  % Mass flowrate of nozzle 
   m_nz = Cd * sqrt((2 * gamma / (gamma - 1)) * P(25) * rho1 * (1 - (P(25) / p_hst)^((gamma - 1) / gamma)) * (P(25) / p_hst)^(2 / gamma)); 

  % Equation
   dPdt(25) =-c^2 * (3*m(25) - 4*m(24) + m(23)) / (2*dx);
   dmdt(25) =-(2*m(25)*c^2/P(25)*((- 3*m(25)+4*m(26)-m(27)) / (2 * dx)))+(c^2*m(25)^2/P(25)^2)*((- 3*P(25)+ 4*P(26)-P(27)) / (2*dx)) - ((-3*P(25)+4*P(26)-P(27)) / (2*dx)) - (f * m(25) * abs(m(25)) * c^2) / (2 *D* (P(25))); % Forward discretization for dm/dt at i=26
   dmdt(1) =-(2*m(1)*c^2/P(1)*((- 3*m(1)+4*m(2)-m(3)) / (2 * dx)))+(c^2*m(1)^2/P(1)^2)*((- 3*P(1)+ 4*P(2)-P(3)) / (2*dx)) - ((-3*P(1)+4*P(2)-P(3)) / (2*dx)) - (f * m(1) * abs(m(1)) * c^2) / (2 *D* (P(1))); % Forward discretization for dm/dt at i=1
    
   % Update m(Nx) based on the time t and the specified conditions
    for i = 1:length(t)
        if t(i) < 600
            m(Nx) = 0;
        elseif t(i) < 1800
            m(Nx) = 509.28; %509.28 kg/s*m^2;
       else 
            m(Nx) = 0;
        end
    end
   
    % Equations
    dPdt(Nx) =-c^2 * (3*m(Nx) - 4*m(Nx - 1) + m(Nx-2)) / (2*dx); % Backward discretization of dP/dt at i=Nx
    
 % Central discretization for in-between grid points (i=2 to 49)
for i = 2:49
    if i == 25
        % Skip the calculation for dPdt(25) since it's already calculated outside the loop
        % dPdt(25) = ...; 
        % dmdt(25) =.....;
%     elseif i == 24
%           m(i)= m(25)+m_nz;% Two mass flow rate enter m(25) and m_nz 
    else
        % Calculate dP/dt and dm/dt for in-between grid points using central discretization
        dPdt(i) = -c^2 * (m(i+1) - m(i-1)) / (2 * dx); % Central discretization of dP/dt
        dmdt(i) =-2*m(i)*c^2/P(i)*((m(i + 1) - m(i - 1)) / (2 * dx))+(c^2*m(i)^2/P(i)^2)*(P(i+1) - P(i-1)) / (2*dx) - (P(i+1) - P(i-1)) / (2*dx) - (f * m(i) * abs(m(i))* c^2 ) /( 2*D*P(i)); % Central discretization of dm/dt

    end
end

    % Combine dPdt and dmdt into a single vector
    dudt = [dPdt; dmdt];
end

[Walter Roberson]

I guarantee that your m(Nx) is not implemented correctly.

[MAYUR MARUTI DHIRDE]

Could you please show me right way to implement m(Nx)?

If you do so it will be great help.

[Torsten]

I tried to explain several times now that your discretization for the acoustic equations is sure to fail. Why do you still waste your time with it and don't concentrate on the upwind schemes ?

And this loop is possible, but it doesn't make sense since the parameter "t" as input to your function is one single value, namely the time at which the differential solver has actually arrived in the integration process. Thus length(t) is always 1.

[Walter Roberson]

I repeat:

You need to divide your integration up into three parts: 0 to 600, 600 to 1800, then over 1800. You will need three different ode*() calls with three different tspan for this. You would parameterize the call to myODEs according to what the current m(Nx) value should be for that tspan.

[t1, u1] = ode15s(@(t, u) MyODEs(t, u, Nx, dx, f, c, D, 0), [0 600], u0,options);
[t2, u2] = ode15s(@(t, u) MyODEs(t, u, Nx, dx, f, c, D, 509.28), [600 1800], u1, options);
[t3, u3] = ode15s(@(t, u) MyODEs(t, u, Nx, dx, f, c, D, 0), [1800 tspan(end)], u2, options);
t = [t1; t2(2:end); t3(2:end)];
u = [u1; u2(2:end,:); u3(2:end,:)];


function dudt = MyODEs(t, u, Nx, dx, f, c, D, mNx)

...
    
   %% Update m(Nx) based on the time t and the specified conditions
   % for i = 1:length(t)
   %     if t(i) < 600
   %         m(Nx) = 0;
   %     elseif t(i) < 1800
   %         m(Nx) = 509.28; %509.28 kg/s*m^2;
   %    else 
   %         m(Nx) = 0;
   %     end
   % end
   
   m(Nx) = mNx;

Any time you have a single call to one of the ode*() routines, and your code has an if or interp1 in it, you should asking yourself whether the code has been carefully constructed to match second derivatives across the different branches of the code; if the answer is No, then the code is quite likely faulty (unless you have reason to know that only one of the code branches can be activated in context.)

[Torsten]

You need to divide your integration up into three parts: 0 to 600, 600 to 1800, then over 1800.

It's advisable to do this (I even provided the code for this splitting, but OP sticks to his/her coding). But the main problem is the unstable discretization of the underlying PDEs.

[MAYUR MARUTI DHIRDE]

Dear torsten,

The equation that, I was trying to solve was not related to acoustic, they are related to the unsteady compressible flow of the pipeline. In an earlier case in which you suggested using an upwind scheme instead of central discretization, there was a typo error in the central discretized scheme of the equation. After correction, the result has shown 100% matches with the result of research from which the equation and boundary condition were taken.

[Torsten]

The equation that, I was trying to solve was not related to acoustic, they are related to the unsteady compressible flow of the pipeline.

In my opinion, your equations model unsteady weakly compressible (barotropic) flow in pipelines. Here, p is a function of rho, but dp/dT is considered as 0.

And maybe these are different applications for the same kind of equations.

continuity:

dP/dt + c^2*dm/dx = 0 

Momentum:

dm/dt + dP/dx + f*c^2*m^2/(2*P*D)=0 

are the equations (4) under

http://www.clawpack.org/riemann_book/html/Acoustics.html

and are called "linear acoustics equations" there (the source term f*c^2*m^2/(2*P*D) is missing, but it's not contributing to instability).

I'd be interested in the test case that gives a stable solution with the discretization scheme you use.

[Walter Roberson]

Could you remind me what the difference is between a pressure wave and an acoustic wave?

[MAYUR MARUTI DHIRDE]

Below is the correct form of code, the catch was in this term f*c^2*m^2/(2*P*D), it was not m^2 , it was m * abs(m)

Nx = 50;         % Number of grid points
L = 5000;        % Length of the spatial domain
dx = L / (Nx-1); % Spatial step size
x = linspace(0, L, Nx); % Spatial points
A=0.1963;
f = 0.008;
D = 500e-3;
rho=1.2;
c = 348.5;
tspan = [0, 3600]; % Time span

% Initial conditions
P0 = 5e6*ones(Nx,1); % Initial condition for pressure (P)
m0 = zeros(Nx, 1); % Initial condition for mass flow rate (m)

% Initial conditions vector
u0 = [P0; m0];
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);

% Solve the system of ODEs using ode15s
[t, u] = ode15s(@(t, u) MyODEs(t, u, Nx, dx, f, c, D), tspan, u0,options);

% Extract the solutions for P and m
P_solution = u(:, 1:Nx);
m_solution = u(:, Nx + 1:end);


%% when m =kg/s*m^2
Q_n= (m_solution./rho); 
Q=(Q_n .*A); 

nodes_to_plot = 1:7:Nx;

% Plotting P_solution for selected nodes with transparency
figure;
hold on;
colors = rand(length(nodes_to_plot), 3);
for idx = 1:length(nodes_to_plot)
    i = nodes_to_plot(idx);
    plot(t, P_solution(:, i), 'LineWidth', 0.5, 'DisplayName', ['Node ' num2str(i)], 'Color', colors(idx,:), 'LineStyle', '-');
end
hold off;
xlabel('Time (t)');
ylabel('Pressure (P)');
title('Pressure (P) at Selected Nodes over Time');
legend('show');

% Choose Node 1 (x = 1) for plotting
node1_idx =1;

% % Extract the pressure and Q values at Node 1
Q_node1 = Q(:, node1_idx);
figure;
plot(t,  Q_node1 , 'LineWidth', 2);
xlabel('Time (t)');
ylabel('Q at Node 1');
title('Q at Node 1 over Time without covective term');

% Set up the system of ODEs including the boundary condition ODEs
function dudt = MyODEs(t, u, Nx, dx, f, c, D)
    P = u(1:Nx);
    m = u(Nx+1:end);
    dPdt = zeros(Nx, 1);
    dmdt = zeros(Nx, 1);
  
    % Apply the specified boundary conditions directly within the ODE function
    P(1) = 5e6; % Set P(1) = 5 MPa for all time steps
    
    % Update m(Nx) based on the time t and the specified conditions
    for i = 1:length(t)
        if t(i) < 600
            m(Nx) = 0;
        elseif t(i) < 1800
            m(Nx) = 509.28; %509.28 kg/s*m^2; %100.008 kg/s 
       else 
            m(Nx) = 0;
        end
    end

     %Continuity equation (dP/dt) for P = Nx 
     dPdt(Nx) = -c^2 * (3*m(Nx) - 4*m(Nx - 1) + m(Nx-2)) / (2*dx);  % Eq at P(Nx) - Backward discretization
          
     %Momentum equation (dm/dt) for m = 1
     dmdt(1) = -(-3*P(1) + 4*P(2) - P(3)) / (2*dx) - (f * m(1) * abs(m(1)) * c^2) / (2*D* P(1));% Eq at m(1) forward discretization ;
    
    % Interior points for Pressure and Mass Flow Rate
    for i = 2:Nx - 1
        dPdt(i) = -c^2 * (m(i+1) - m(i-1)) / (2 * dx);
        dmdt(i) = -(P(i+1) - P(i-1)) / (2 * dx) - f * m(i) * abs(m(i)) * c^2 / ( 2*D * P(i));   
    end
  
    % Combine dPdt and dmdt into a single vector
    dudt = [dPdt; dmdt];
end

[MAYUR MARUTI DHIRDE]

Dear walter,

Pressure wave : Any disturbance that results in variations in the pressure within a medium is referred to as a pressure wave in general.

Acoustic wave : An acoustic wave is a particular kind of pressure wave that is associated with sound waves traveling through a medium, especially in the audible frequency range.

[Torsten]

And you think this zigzagging in the results is physical ?

[MAYUR MARUTI DHIRDE]

That zigzagging has a physical meaning , those are fluctuation amplitude and damping prior to reaching the flow to its steady state condition.

If you are intersted to know more in detail you can look in to this paper "DYNAMIC MODELING OF NON-ISOTHERMAL GAS PIPELINE SYSTEMS"

[Walter Roberson]

    for i = 1:length(t)
        if t(i) < 600
            m(Nx) = 0;
        elseif t(i) < 1800
            m(Nx) = 509.28; %509.28 kg/s*m^2; %100.008 kg/s 
       else 
            m(Nx) = 0;
        end
    end

You posted asking why you get an error message about not being able to meet integration tolerances.

You are most directly getting that error message because of the above code, and I already showed you how to fix that problem.

If you do not fix that code to account for the mathematics of Runge Kutta solvers, then you are destined to fail, even if your PDE discretization were fine.

Your code cannot work as long as that problem continues to exist in the code.

Runge Kutta methods are based upon successive piecewise approximation to ODEs, with the weighting coefficients used having been carefully chosen for ease of calculation. However, the internal weighting coefficients that are used are only valid for equations that are at least . Your method of assigning m(Nx) is not so your equations are not valid for use with any of the mathworks ode*() variable-step solver routines.

[MAYUR MARUTI DHIRDE]

I don't understand one thing, that if the problem is in implementation of m (Nx) boundary condition, then it should show an inconsistency error if we comment this part of the code.

elseif i == 24
m(i)= m(25)+m_nz;% Two mass flow rate enter m(25) and m_nz 

If we are commenting this part, then why the code is running?

[Walter Roberson]

As I wrote above,

When you are fortunate, ode15s will notice the discontunity and give you an error message about failing on integration tolerances. When you are not fortunately, ode15s will fail to notice and will keep going and confidently give you the wrong answer.

[Torsten]

Here is an implementation of @Walter Roberson idea.

You should start making further changes from this code because the discontinuities at t = 600 and t = 1800 are now taken out of the continuous integration.

A problem is that you already use m(24) before you define it as m(25)+m_nz.

And after you come back from the integrator, you must change all values in the solution vector u you explicitly set in the function "MyODEs" because for the solver, their rate of change is 0 and they thus remain at their initial values.

 
Nx = 50;         % Number of grid points
L = 5000;        % Length of the spatial domain
dx = L / (Nx-1); % Spatial step size
x = linspace(0, L, Nx); % Spatial points
A=0.1963;
f = 0.008;
D = 500e-3;
rho=1.2;
c = 348.5;
P_atm = 1.013e5; % Atmospheric pressure in Pa
tspan = [0, 3600]; % Time span
tcut = [tspan(1) 600 1800 tspan(end)];
mcut = [0 509.28 0];

options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);

t_solution = [];
P_solution = [];
m_solution = [];

for i = 1:numel(mcut)
  tspan = [tcut(i) tcut(i+1)] ; 
  
  if i == 1
     % Initial conditions
     P0 = 5e6*ones(Nx,1); % Initial condition for pressure (P)
     m0 = zeros(Nx, 1); % Initial condition for mass flow rate (m)
     m0(end) = mcut(1);
  else
     P0 = u(end, 1:Nx).';
     m0 = u(end, Nx + 1:end).';
     m0(end) = mcut(i); 
  end
  % Initial conditions vector
  u0 = [P0;m0];
  
  % Solve the system of ODEs using ode15s
  [t, u] = ode15s(@(t, u) MyODEs(t, u, Nx, dx, f, c, D), tspan, u0,options);

  % Extract the solutions for P and m
  t_solution = [t_solution;t];
  P_solution = [P_solution;u(:, 1:Nx)];
  m_solution = [m_solution;u(:, Nx + 1:end)];
  
  if t(end) < tspan(2) 
      break
  end
end

Warning: Failure at t=7.977123e+02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.818989e-12) at time t.

%% when m =kg/s*m^2
Q_n= (m_solution./rho); 
Q=(Q_n .*A);

% nodes_to_plot = 1:7:Nx;
% 
% % Plotting P_solution for selected nodes with transparency
% figure;
% hold on;
% colors = rand(length(nodes_to_plot), 3);
% for idx = 1:length(nodes_to_plot)
%     i = nodes_to_plot(idx);
%     plot(t, P_solution(:, i), 'LineWidth', 0.5, 'DisplayName', ['Node ' num2str(i)], 'Color', colors(idx,:), 'LineStyle', '-');
% end
% hold off;
% xlabel('Time (t)');
% ylabel('Pressure (P)');
% title('Pressure (P) at Selected Nodes over Time');
% legend('show');

% % Choose Node 1 (x = 1) for plotting
 node1_idx = 25;

% % Extract the pressure and Q values at Node 1
Q_node1 = Q(:, node1_idx);
figure;
plot(t_solution,  Q_node1 , 'LineWidth', 2);
xlabel('Time (t)');
ylabel('Q at Node 1');
title('Q at Node 1 over Time with covective term');

% Set up the system of ODEs including the boundary condition ODEs
function dudt = MyODEs(t, u, Nx, dx, f, c, D)
    P = u(1:Nx);
    m = u(Nx + 1:end);
    dPdt = zeros(Nx, 1);
    dmdt = zeros(Nx, 1);

    % Boundary Condition
    P(1) = 5e6; % P(1) boundary condition
    
    % Nozzle at i=25 (Algebraic relationships)
    p_hst = 503.51e3; % p_hst = rho * g * d (where d = 40m)
    Cd = 0.55;
    gamma = 1.4;
    rho1 = 1.2;
   
  % Mass flowrate of nozzle 
   m_nz = Cd * sqrt((2 * gamma / (gamma - 1)) * P(25) * rho1 * (1 - (P(25) / p_hst)^((gamma - 1) / gamma)) * (P(25) / p_hst)^(2 / gamma)); 

  % Equation
   dPdt(25) =-c^2 * (3*m(25) - 4*m(24) + m(23)) / (2*dx); % Backward discretization of dP/dt at i=25
   dmdt(25) =-(2*m(25)*c^2/P(25)*((- 3*m(25)+4*m(26)-m(27)) / (2 * dx)))+(c^2*m(25)^2/P(25)^2)*((- 3*P(25)+ 4*P(26)-P(27)) / (2*dx)) - ((-3*P(25)+4*P(26)-P(27)) / (2*dx)) - (f * m(25) * abs(m(25)) * c^2) / (2 *D* (P(25))); % Forward discretization for dm/dt at i=26
   dmdt(1) =-(2*m(1)*c^2/P(1)*((- 3*m(1)+4*m(2)-m(3)) / (2 * dx)))+(c^2*m(1)^2/P(1)^2)*((- 3*P(1)+ 4*P(2)-P(3)) / (2*dx)) - ((-3*P(1)+4*P(2)-P(3)) / (2*dx)) - (f * m(1) * abs(m(1)) * c^2) / (2 *D* (P(1))); % Forward discretization for dm/dt at i=1
    
   % Update m(Nx) based on the time t and the specified conditions
   % for i = 1:length(t)
   %     if t(i) < 600
   %         m(Nx) = 0;
   %     elseif t(i) < 1800
   %         m(Nx) = 509.28; %509.28 kg/s*m^2;
   %    else 
   %         m(Nx) = 0;
   %     end
   % end
   
    % Equations
    dPdt(Nx) =-c^2 * (3*m(Nx) - 4*m(Nx - 1) + m(Nx-2)) / (2*dx); % Backward discretization of dP/dt at i=Nx
    
 % Central discretization for in-between grid points (i=2 to 49)
for i = 2:49
    if i == 25
        % Skip the calculation for dPdt(25) since it's already calculated outside the loop
        % dPdt(25) = ...; 
        % dmdt(25) =.....;
    elseif i == 24
         m(i)= m(25)+m_nz;% Two mass flow rate enter m(25) and m_nz 
         
    else
        % Calculate dP/dt and dm/dt for in-between grid points using central discretization
        dPdt(i) = -c^2 * (m(i+1) - m(i-1)) / (2 * dx); % Central discretization of dP/dt
        dmdt(i) =-2*m(i)*c^2/P(i)*((m(i + 1) - m(i - 1)) / (2 * dx))+(c^2*m(i)^2/P(i)^2)*(P(i+1) - P(i-1)) / (2*dx) - (P(i+1) - P(i-1)) / (2*dx) - (f * m(i) * abs(m(i))* c^2 ) /( 2*D*P(i)); % Central discretization of dm/dt

    end
end

    % Combine dPdt and dmdt into a single vector
    dudt = [dPdt; dmdt];
end

[Torsten]

I don't understand one thing, that if the problem is in implementation of m (Nx) boundary condition, then it should show an inconsistency error if we comment this part of the code.

Why do you think so ?

If you comment this part, dmdt(24) and dPdt(24) are computed from the expressions in the "else" statement:

        % Calculate dP/dt and dm/dt for in-between grid points using central discretization
        dPdt(i) = -c^2 * (m(i+1) - m(i-1)) / (2 * dx); % Central discretization of dP/dt
        dmdt(i) =-2*m(i)*c^2/P(i)*((m(i + 1) - m(i - 1)) / (2 * dx))+(c^2*m(i)^2/P(i)^2)*(P(i+1) - P(i-1)) / (2*dx) - (P(i+1) - P(i-1)) / (2*dx) - (f * m(i) * abs(m(i))* c^2 ) /( 2*D*P(i)); % Central discretization of dm/dt

If you don't comment this part, dmdt(24) and dPdt(24) are taken from their initialization

    dPdt = zeros(Nx, 1);
    dmdt = zeros(Nx, 1);    

thus remain unchanged in the course of the computation, and the setting

    m(i)= m(25)+m_nz;    

doesn't influence anything.

[MAYUR MARUTI DHIRDE]

See in the calculation of dPdt(25) =-c^2 * (3*m(25) - 4*m(24) + m(23)) / (2*dx); the m(25) is current value and m(24) and m(23) are previously calculated value mainly using initial condition. But while performing calculation for dPdt(24) and dmdt(24) , the m(24) should make use of this condition m(24)=m(25)+m_nz , in which m(25) and m_nz are previously calculated value.

[Torsten]

After the lines

P = u(1:Nx);
m = u(Nx + 1:end);

manipulate the values of P and m at the node positions where needed.

So you keep control about which values of P and m you use in the sequel to compute the time derivatives.

Since you don't use the P and m values taken from the u-array at the manipulated positions, the time derivatives at these positions should be set to zero in order not to possibly hinder convergence of the time integrator.

[MAYUR MARUTI DHIRDE]

% Set up the system of ODEs including the boundary condition ODEs
function dudt = MyODEs(t, u, Nx, dx, f, c, D)
    P = u(1:Nx);
    m = u(Nx + 1:end);
    
    dPdt = zeros(Nx, 1);
    dmdt = zeros(Nx, 1);

    % Calculate dPdt(25) using previous known m(24)
    dPdt(25) = -c^2 * (3 * m(25) - 4 * m(24) + m(23)) / (2 * dx);   
    
    % Boundary Condition
    P(1) = 5e6; % P(1) boundary condition
 
    % Nozzle at i=25 (Algebraic relationships)
    p_hst = 503.51e3; % p_hst = rho * g * d (where d = 40m)
    Cd = 0.55;
    gamma = 1.4;
    rho1 = 1.2;
    % Mass flowrate of nozzle 
    m_nz = Cd * sqrt((2 * gamma / (gamma - 1)) * P(25) * rho1 * (1 - (P(25) / p_hst)^((gamma - 1) / gamma)) * (P(25) / p_hst)^(2 / gamma)); 

    % Manipulate m at node position 24
    manipulated_m_24 = m(25) + m_nz; 
    m(24) = manipulated_m_24; 
    
    % Setting derivatives at manipulated positions to zero
    dPdt(24) = 0;
    dmdt(24) = 0;

    % Calculation for dmdt(25) and dmdt(1)
    dmdt(25) = -(2 * m(25) * c^2 / P(25) * ((-3 * m(25) + 4 * m(26) - m(27)) / (2 * dx))) + (c^2 * m(25)^2 / P(25)^2) * ((-3 * P(25) + 4 * P(26) - P(27)) / (2 * dx)) - ((-3 * P(25) + 4 * P(26) - P(27)) / (2 * dx)) - (f * m(25) * abs(m(25)) * c^2) / (2 * D * (P(25))); % Forward discretization for dm/dt at i=25
    dmdt(1) = -(2 * m(1) * c^2 / P(1) * ((-3 * m(1) + 4 * m(2) - m(3)) / (2 * dx))) + (c^2 * m(1)^2 / P(1)^2) * ((-3 * P(1) + 4 * P(2) - P(3)) / (2 * dx)) - ((-3 * P(1) + 4 * P(2) - P(3)) / (2 * dx)) - (f * m(1) * abs(m(1)) * c^2) / (2 * D * (P(1))); % Forward discretization for dm/dt at i=1
    
    % Calculation for dPdt(Nx)
    dPdt(Nx) = -c^2 * (3 * m(Nx) - 4 * m(Nx - 1) + m(Nx - 2)) / (2 * dx); % Backward discretization of dP/dt at i=Nx

    % Loop for calculating derivatives from i=2 to i=49
    for i = 2:49
        if i == 24 || i == 25
            % Skip calculation or set to zero if needed for i=24 and i=25
            dPdt(i) = 0; 
            dmdt(i) = 0;
        else
            % Calculate dP/dt and dm/dt for other grid points using your scheme
            dPdt(i) = -c^2 * (m(i + 1) - m(i - 1)) / (2 * dx);  % Calculation for dP/dt at other points
            dmdt(i) = -2 * m(i) * c^2 / P(i) * ((m(i + 1) - m(i - 1)) / (2 * dx)) + (c^2 * m(i)^2 / P(i)^2) * (P(i + 1) - P(i - 1)) / (2 * dx) - (P(i + 1) - P(i - 1)) / (2 * dx) - (f * m(i) * abs(m(i)) * c^2) / (2 * D * P(i)); % Calculation for dm/dt at other points
        end
    end

    % Combine dPdt and dmdt into a single vector
    dudt = [dPdt; dmdt];
end

%%
Warning: Failure at t=3.776956e+02.  Unable to meet integration tolerances without reducing the step size below the
smallest value allowed (9.094947e-13) at time t. 

[Torsten]

You first set

dPdt(25) = -c^2 * (3 * m(25) - 4 * m(24) + m(23)) / (2 * dx)

, but then reset it to 0 in

if i == 24 || i == 25
    % Skip calculation or set to zero if needed for i=24 and i=25
    dPdt(i) = 0; 
    dmdt(i) = 0;
else
    ...       

If you set dPdt(25) = 0, you must fix a value for P(25) in your code since P(25) won't remain at its initial value of 5e6.