ODE15s Unable to meet integration tolerances without reducing the step size below the smallest value allowed

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 (0 to 3600 seconds)

% 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=7.137117e+01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.273737e-13) at time t.

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

Q_n= (m_solution./rho); %m/sec
Q=(Q_n .*A);% m^3/sec

figure;
plot(t,Q , 'LineWidth', 2);
xlabel('Time (t)');
ylabel('Q');
title('Q over Time');

figure;
plot(t,P_solution , 'LineWidth', 2);
xlabel('Time (t)');
ylabel('P');
title('P over Time');

% Define the boundary conditions as separate ODEs
function dPdt = PBoundaryODE(t, P)
    % This ODE defines how P(1, t) evolves with time
    dPdt = 0; 
end

function dmNdt = mBoundaryODE(t, m)
    % This ODE defines how m(Nx, t) evolves with time
    if t < 60
        dmNdt = 0;
    elseif t < 1800
        dmNdt = 509.10;
    else
        dmNdt = 0;
    end
end

% 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);
    dudt = zeros(Nx * 2, 1); % Initialize the derivative vector

    % Continuity equation (dP/dt) for P= Nx 
    dudt(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
    dudt(Nx + 1) = -(4*P(2) - 3*P(1)-P(3)) / (2*dx) - (f * m(1)^2 * 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
    dudt(i) = -c^2 * (m(i+1) - m(i-1)) / (2 * dx);
    dudt(i + Nx) = -(P(i+1) - P(i-1)) / (2 * dx) - f * m(i) * m(i) * c^2 / (2 * D * P(i));
end

    % Pressure boundary condition ODE at i = 1
    dudt(1) = PBoundaryODE(t, u(1));

    % Mass flow rate boundary condition ODE at m(Nx, t)
    dudt(2*Nx) =mBoundaryODE(t, u(2*Nx));  % Mass flow rate boundary condition ODE at m(Nx, t)

end

%% why boundary condition should be include in main function?
% So, even though you have separate functions for the boundary conditions, 
% you still need to include their contributions within the main ODE function MyODEs. 
% The reason for this is that the ODE solvers in MATLAB expect a single function that defines the entire system of ODEs.

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

Answer 1

[Torsten]

If you store P for 1:Nx and m for Nx+1:2*Nx, you should also return dPdt for 1:Nx and dmdt for Nx+1:2*Nx. But you return dPdt for 1 and Nx+2:2*Nx-1 and dmdt for 2:Nx and 2*Nx, don't you ?

[MAYUR MARUTI DHIRDE]

could you please explain what changes should be made in above code to adabt your suggestion?

I tried but somehow, i am getting same error.

[Torsten]

dudt(1:Nx) must be dPdt(1:Nx), dudt(Nx+1:2*Nx) must be dmdt(1:Nx).

So maybe you should work with arrays dPdt and dmdt of size (Nx,1), and the last step in your function "MyODEs" would be to set

dudt = [dPdt;dmdt]

[Torsten]

I just looked at your equations here

https://uk.mathworks.com/matlabcentral/answers/2032419-nan-error-while-solving-newton-raphson-method?s_tid=srchtitle

and saw that they are

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

In this case, your discretization is (at least formally) correct.

So in the case above, you try to set P = 5 bar at the inlet and m = 509.1*(t-60) for 60 <=t<=1800 ? This will produce a strong negative pressure in the tube, doesn't it ? And so I think your solution is not that bad for this case.

Of course you must keep in mind that there is a maximum mass flow density given by the speed-of-sound.

[MAYUR MARUTI DHIRDE]

The equation are

continuity:

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

Momentum:

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

I have boundary condition as P(1,t)= 5MPa and m(Nx,t)= f(t) and f(t) is shown in graphs which varies in time. For 0<t<600 sec, m=0 For 600<t<1800sec,m=509.10 and For 1800<t<3600 sec, m=0

To solve in between grid points I have for loop from i=2:Nx-1 for dP/dt and dm/dt.Further to solve the last grid point of dP(Nx)/dt, i have Backward discretization equation and for dm(1)/dt i have forward discretization.

[Torsten]

For 600<t<1800sec,m=509.10 and For 1800<t<3600 sec, m=0

Your code sets at x = L

dmdt = 0 , thus m = 0, for 0 <= t < 60,

dmdt = 509.1, thus m = 509.1*(t-60) for 60 <= t < 1800 and

dmdt = 0, thus m = 509.1*(1800-60) for 1800 <= t <=3600.

[MAYUR MARUTI DHIRDE]

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 (0 to 3600 seconds)

% 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);

Q_n= (m_solution./rho); %m/sec
Q=(Q_n .*A);% m^3/sec

figure;
plot(t,Q , 'LineWidth', 2);
xlabel('Time (t)');
ylabel('Q');
title('Q over Time');

figure;
plot(t,P_solution , 'LineWidth', 2);
xlabel('Time (t)');
ylabel('P');
title('P over Time');

% Define the boundary conditions as separate ODEs
function dPdt = PBoundaryODE(t, P)
    % This ODE defines how P(1, t) evolves with time
    dPdt = 0; % You can specify the time evolution of P(1, t) here
end

function dmNdt = mBoundaryODE(t, m)
    % This ODE defines how m(Nx, t) evolves with time
    if t < 600
        dmNdt = 0;
    elseif t < 1800
        dmNdt = 509.1 * (t - 600);
    else
        dmNdt = 0;
    end
end

% 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);

    % 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)^2 * 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) * m(i) * c^2 / (2 * D * P(i));
    end

    % Pressure boundary condition ODE at i = 1
    dPdt(1) = PBoundaryODE(t, P);
   
    % Mass flow rate boundary condition ODE at m(Nx, t)
    dmdt(Nx) = mBoundaryODE(t, m);
    
    % Combine dPdt and dmdt into a single vector
    dudt = [dPdt; dmdt];
end

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

[Torsten]

It doesn't work (most probably because of the discretization), but the setup is correct.

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;
tend = [0,600,1800,3600];
mN = [0,509.1,0];
% 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) = mN(1);
for i = 1:numel(mN)
  tspan = [tend(i),tend(i+1)]; % Time span (0 to 3600 seconds)
  % 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
  t_solution{i} = t;
  P_solution{i} = u(:, 1:Nx);
  m_solution{i} = u(:, Nx + 1:end);
  Q_n{i} = (m_solution{i}./rho); %m/sec
  Q{i} = Q_n{i}*A;% m^3/sec
  P0 = u(end, 1:Nx);
  m0 = u(end, Nx + 1:end)
  m0(end) = mN(i+1);
end
figure(1)
hold on
for i = 1:numel(mN)
  plot(t_solution{i},Q{i} , 'LineWidth', 2);
  xlabel('Time (t)');
  ylabel('Q');
  title('Q over Time');
end
hold off
figure(2)
hold on
for i = 1:numel(mN)
  plot(t_solution{i},P_solution{i} , 'LineWidth', 2);
  xlabel('Time (t)');
  ylabel('P');
  title('P over Time');
end
hold off

% Define the boundary conditions as separate ODEs
function dPdt = PBoundaryODE(t, P)
    % This ODE defines how P(1, t) evolves with time
    dPdt = 0; % You can specify the time evolution of P(1, t) here
end

% 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);
    % 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)^2 * 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) * m(i) * c^2 / (2 * D * P(i));
    end
    % Pressure boundary condition ODE at i = 1
    dPdt(1) = PBoundaryODE(t, P);
    
    % Combine dPdt and dmdt into a single vector
    dudt = [dPdt; dmdt];
end

[MAYUR MARUTI DHIRDE]

%% I have made changes in the code, the code is almost correct the only problem which i am getting is,
Warning: Failure at t=1.935312e+03.  Unable to meet integration tolerances without reducing the step size below the
smallest value allowed (3.637979e-12) at time t. 
How Should i fix this ? 




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 (0 to 3600 seconds)

% 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-9, 'AbsTol', 1e-9);

% 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);

Q_n= (m_solution./rho); %m/sec
Q=(Q_n .*A);% m^3/sec

figure;
plot(t,Q , 'LineWidth', 2);
xlabel('Time (t)');
ylabel('Q');
title('Q over Time');

figure;
plot(t,P_solution , 'LineWidth', 2);
xlabel('Time (t)');
ylabel('P');
title('P over Time');

% 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.3;
        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)^2 * 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) * m(i) * c^2 / (2 * D * P(i));
    end
    % Combine dPdt and dmdt into a single vector
    dudt = [dPdt; dmdt];
end

[Torsten]

Did you look at the result curves at t = 1.935312e+03 ? They usually give an indication why the integrator gives up and can no longer keep the promised error tolerances.

You won't get a physical solution for this problem with the discretization scheme you use. A mathematical treatment of your equation and a code to solve it is given here:

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

[MAYUR MARUTI DHIRDE]

Yes, you are correct with ode15s ,I cannot get physical solution for this problem. I need to use another method " Fully Implicit Finite Difference Method (FIFDM), "

[Torsten]

ode15s (and thus the time discretization) is not the problem, it's the spatial discretization.

Did you hear about upwind methods for hyperbolic systems of partial differential equations ? They are explained in detail in the link I gave you (and the book by LeVeque that is mentionned on that page).

Try the "clawpack" package for your problem. The acoustic equations are already part of the set of examples.

https://www.clawpack.org/

[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);
  
    % 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.3; %
        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(3) + 4*P(2) - P(1)) / (2*dx) - (f * m(1)^2 * c^2) / (2*D * P(1)); % Eq at m(1) forward discretization ;
    
 % Implement upwind scheme for spatial derivatives
 %% second order accuracy    
 for i = 2:Nx - 1
    % Handling boundary points differently
    if i == 2
        % Forward differencing for i=2
        dPdt(i) = -c^2 * (-3*m(i + 2) + 4*m(i + 1) - m(i)) / (2 * dx);
        dmdt(i) = -(-3*P(i+2) + 4*P(i+1) - P(i)) / (2*dx) - (f * m(i)^2 * c^2) / (2*D * P(i));
    else
        % Backward differencing for other interior points
        dPdt(i) = -c^2 * (3*m(i) - 4*m(i - 1) + m(i - 2)) / (2 * dx);
        dmdt(i) = -(3*P(i) - 4*P(i-1) - P(i-2)) / (2*dx) - (f * m(i)^2 * c^2) / (2*D * P(i));
    end
 end

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


Is this correct way of implementing upwind scheme for spatial part? 
With this i wont get warning and solution runs for tspan from 0 t0 3600
But results are varying. 

[Torsten]

Is this correct way of implementing upwind scheme for spatial part?

No, it's not correct. Your hyperbolic system has two eigenvalues: c and -c. For these eigenvalues, you have to determine the corresponding eigenvectors. Then you have to diagonalize your PDE system with the help of the eigenvectors. Then the characteristic variable corresponding to the negative eigenvalue has to be approximated from the right while the one corresponding to the positive eigenvalue has to be approximated from the left when forming the spatial derivative. I don't remember the whole process in detail, but you must invest a lot of effort if you are not used to this process. You cannot shake a scheme out of your sleeve and hope it is correct. That's why I adviced you to use the CLAWPACK software in order to prevent you from coding when better solutions already exist.

[MAYUR MARUTI DHIRDE]

I appreciate your assistance. I will take into account your suggestion and work on it.