Solving ODE with explicit equations

I'm trying to solve a system consisting of two ODEs and two explicit algebraic equations. The required (main) variables are deltaP2(t), deltaP1(t), V1(t) and Rf(t) (actually i also need to output the derivative dRf/dt as well). t (time) - independent variable. The program says there are not enough input arguments. Therefore, I think that I somehow incorrectly specified the vector with the required variables or the function output. Will appreciate any help with starting this program!

Function file code:

function f = SpherODEs(t,Y)
DeltaP2 = Y(1);
Rf = Y(2);
V1 = Y(3);
DeltaP1 = Y(4);

%Constants
gamma1 = 1.4;
gamma2 = 1.27;
C0 = 340;
sigma0 = 7;
SL = 1.37*10^-3;
n = 1;

%Additioal conditions/parameters
Af = 1;
if and(t >= 0.05, t < 1.5)
    k = 0.05;
elseif and(t < 0.05, t >= 1.5)
    k = 0;
end
Ug = n + k * V1;
sigma = sigma0 * ((1 + DeltaP1)^(1/gamma1)) / ((1 + DeltaP2)^(1/gamma2));
bettaV = 0;
dbettaVdt = 0;

%Differential equations
dDeltaP2dt = ((3 * gamma2 * (1 + DeltaP2)) / Rf) * ((Af * Ug * sigma) / bettaV - bettaV * (Af * Ug + V1) - Rf / 3 * dbettaVdt);
dRfdt = Af * Ug + V1;

%Explicit equations
DeltaP1 = ((sigma - 1) * Rf * SL^2) / (sigma * (1 + dRfdt * SL / gamma1^1/2)) * ((3/2 * dRfdt^2) + Rf * 0) + 0;
DeltaP1 = DeltaP2 + (sigma - 1) * (1 + DeltaP1)^(1/gamma1) * Af^2 * Ug^2 * SL^2;

%Output
f = [dDeltaP2dt; dRfdt; V1; DeltaP1];
end

Script file code:

clc
tspan = [0 2]; %Integration time span
y0 = [0 0.001 0 0]; %Initial conditions 
[t y] = ode45 (@SpherODEs, tspan, y0);

Unrecognized function or variable 'k'.

Error in solution>SpherODEs (line 22)
Ug = n + k * V1;

Error in odearguments (line 93)
f0 = ode(t0,y0,args{:});

Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

plot (t, Y(:,1));
legend ('DeltaP2(t)')
ylabel ('DeltaP2');
xlabel ('t');
axis ([0 2 0 100]);
title ('Dimensionless overpressure by time')

Answer 1

[Walter Roberson]

if and(t >= 0.05, t < 1.5)
    k = 0.05;
elseif and(t < 0.05, t >= 1.5)
    k = 0;
end
Ug = n + k * V1;

That code only defines k within certain time ranges. If t does not match any of those time ranges then k will not be defined.

and(t < 0.05, t >= 1.5)

Because of the and the outcome of this statement is true if there is a time that is both less than 0.05 and greater than or equal to 1.5. There is no such time possible: a time that is less than 0.05 cannot also be greater than or equal to 1.5.

You wanted an or test there instead of an and test.

On the other hand... since t cannot be NaN, you might as well use a simple else there instead of an elseif

Fixing this reduces the problem. Now you are faced with the issue that your equations are discontinuous at time boundaries 0.05 and 1.5 . If you are lucky, ode45() will notice the discontinuity and will complain about failure to meet integration tolerances somewhere near time 0.05 . If you are not lucky, ode45() will not notice the discontinuity, and will silently give you the wrong solution.

The expressions calculated in the objective function must be continuous and have continuous first and second derivatives to meet the mathematics under which Runge-Kutta calculations are valid. This continuity applies for any one call to ode45()

What you need to do is run three different ode45() calls; the first one with tspan 0 to 0.05. Take the last output, Y(end,:) and use that as the boundary conditions for a second ode45() call, this one tspan 0.05 to 1.5, Take the last output of that, Y(end,:) and use that as the boundary conditons for the third ode45() call, this one tspan 1.5 to 2. In this way, any one ode45() call does not cross a boundary that changes k and everything works out.

[Maxim]

Thank you for your reply! If I use ode15s instead of ode45 (since, apparently, I work with the DAE system), will I need to divide the calculation into three time intervals (in relation to this case with k)? In the future, I will need to set various changes for variables depending on time, not only for k, to simulate the reaction of the system to sudden changes (flow over an obstacle). Thus, there can be quite a lot of time intervals...Or is there some more elegant approach to setting conditions for changing a variable over time or depending on another variable, not through logical operators?

Answer 2

[Torsten]

If you have a mixture of differential and algebraic equations, you need to define the mass matrix appropriately.

In your case - if the first two equations are differential and the last two are algebraic -

tspan = [0 2]; %Integration time span
y0 = [0 0.001 0 0]; %Initial conditions 
M = [1 0 0 0;0 1 0 0;0 0 0 0;0 0 0 0];
options = odeset('Mass',M)
[t, y] = ode15s (@SpherODEs, tspan, y0, options);

Note that you cannot solve a mixture of differential and algebraic equations with ode45 - you need to use ode15s.

[Maxim]

Thank you for your reply! I followed your advice. Now the program outputs the following, but I tried to solve this system in Mathcad Prime, with the current initial conditions, and still got some result . Do I understand correctly that my system is a so-called DAE, which should be designed and solved accordingly (for example, according to "https://ww2.mathworks.cn/help/symbolic/solve-differential-algebraic-equations.html")?

Error using daeic12 (line 148)
Initial conditions are inconsistent. Specify a better guess of consistent initial conditions.

Error in ode15s (line 298)
        [y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(ode,odeArgs,t,ICtype,Mt,y,yp0,f0,...

Error in ScriptForSpherODEs (line 6)
[t y] = ode15s (@SpherODEs, tspan, y0, options);

[Torsten]

When looking at your code, I don't think that your implementation of the algebraic equations is correct.

E.g.

%Output
f = [dDeltaP2dt; dRfdt; V1; DeltaP1];

would set V1 = 0 and DeltaP1 = 0 as algebraic equations. This does not seem correct since V1 and DeltaP1 are your algebraic solution variables, not your algebraic equations.

If your algebraic equations read

g3(deltaP2, Rf,deltaP1, V1) = 0
g4(deltaP2, Rf,deltaP1, V1) = 0

you must return

%Output
f = [dDeltaP2dt; dRfdt; g3(deltaP2, Rf,deltaP1, V1);g4(deltaP2, Rf,deltaP1, V1)];

to ode15s.

Take a look at

Solve Robertson Problem as Semi-Explicit Differential Algebraic Equations (DAEs)

under

https://uk.mathworks.com/help/matlab/math/solve-differential-algebraic-equations-daes.html

for an example.