Solving the system of ODEs and algebraic Equation

I am trying to solve the coupled ODE with dependent algebraic equation using RK4 method. I tried with my code below. I want to plot v, w, p, R, u_s and h_r, v/s, t. But I am getting error "Too many input arguments". I tried to figure out the problem for a day but could not able to do. Can anyone point my mistake? Your help will be much aappreciated.

t0      = deg2rad(0);
a       = deg2rad(60); 
m       = 1; 
d_j     = 0.6; 
We      = 558; 
Re      = 87; 
b       = (-0.13*m^3 + 0.263*m^2 + 0.039*m + 0.330)/tan(a);

%% Initial Conditions for differential equation
u_b0    = 0.1;  
p0      = pi/2; 
R0      = 0.002*We/b; 
s_b0    = 0.02712;

%% Estimating Sheet Velocity at theta = 0 (initial sheet velocity)
x       = b*cos(t0)*sin(a)^2; 
x01     = 4*(1 - (cos(t0)^2)*cos(a)^2); 
x02     = (b*sin(t0)*sin(a))^2;
q_j0    = - (x - (x01 - x02)^0.5)/(x01/4);

syms q 
d       = (b^2*sin(a)^2 + q^2*(1 - cos(t0)^2*cos(a)^2) + 2*b*q*cos(t0)*sin(a)^2)^0.5;
u_j0    = 2*(m-1)*(4*d^2 - 1) + m;
F0      = int(q*u_j0, q, 0, q_j0);
F0      = double(F0);
eq_e    = @(u_s) (R0-q_j0)*F0*sin(a)*u_s.^3 +...
                 (R0.^2+2*R0*sqrt(pi*s_b0)+(16/(sin(a)*3*q_j0.^3*R0*Re))*((R0-q_j0)*F0*sin(a)).^2)*u_s.^2 -...
                ((R0-q_j0)*F0*sin(a)-((F0.^3*sin(a).^3)/(9))*((1/q_j0.^3)-(1/R0.^3))+((R0-q_j0)*F0*sin(a))*(4/(q_j0*sin(a))))*u_s-...
                ((R0-q_j0)*F0*sin(a))*((4*sin(a).^2*F0.^3*(R0.^5-q_j0.^5))/(15*q_j0.^7*R0.^5*Re));
u_s0    = fsolve(eq_e, 0.5)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

u_s0 = 0.5589

%% DAE System
syms q t
q_j     = (-(b*cos(t)*sin(a)^2)+(1-cos(t)^2*cos(a)^2-(b*sin(t)*sin(a))^2)^0.5)/(1-cos(t)^2*cos(a)^2);
u_j     = 2*(m-1)*(4*((q*sin(t))^2 + sin(a)^2*(q*cos(t)+b)^2)-1) + m;
F       = matlabFunction(int(q*u_j,   q, 0, q_j)); 
G       = matlabFunction(int(q*u_j^3, q, 0, q_j));

t       = 0:0.01:pi;
y       = zeros(4, length(t));
y(:,1)  = [u_b0; s_b0; p0; R0];

%% Solve algebraic equation for u_s
for i = 1:length(t)-1  
    q_j     = (-(b*cos(t(i))*sin(a)^2)+(1-cos(t(i))^2*cos(a)^2-(b*sin(t(i))*sin(a))^2)^0.5)/(1-cos(t(i))^2*cos(a)^2);
    S1      = double(F(t(i))); 
    us_eqs  = @(u_s) (y(4,i) - q_j(i))*S1*sin(a)*u_s + (S1*sin(a)).^3*((1/q_j(i).^3) - (1/y(4,i).^3))/(9*u_s) + y(4,i).^2 + 2*y(4,i)*sqrt(pi*y(2,i)) - 4*(y(4,i) - q_j(i))*S1/(q_j(i)*u_s) + (16*sin(a)*S1.^2*(y(4,i) - q_j(i)).^2/(3*q_j(i)*y(4,i)*Re)) + 4*sin(a).^3*S1.^4*(y(4,i) - q_j(i))*(y(4,i).^5-q_j(i).^5)/(15*q_j(i).^7*y(4,i).^5*Re*u_s.^2) - (y(4,i) - q_j(i))*S1*sin(a)/u_s;
    u_s     = fsolve(us_eqs, 1);
    sol     = ode45(@(t, y) odesym(t, y, We, Re, F, G, u_s), [t(i) t(i+1)], y(:,i));
end

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

Index exceeds the number of array elements. Index must not exceed 1.

Error in solution>@(u_s)(y(4,i)-q_j(i))*S1*sin(a)*u_s+(S1*sin(a)).^3*((1/q_j(i).^3)-(1/y(4,i).^3))/(9*u_s)+y(4,i).^2+2*y(4,i)*sqrt(pi*y(2,i))-4*(y(4,i)-q_j(i))*S1/(q_j(i)*u_s)+(16*sin(a)*S1.^2*(y(4,i)-q_j(i)).^2/(3*q_j(i)*y(4,i)*Re))+4*sin(a).^3*S1.^4*(y(4,i)-q_j(i))*(y(4,i).^5-q_j(i).^5)/(15*q_j(i).^7*y(4,i).^5*Re*u_s.^2)-(y(4,i)-q_j(i))*S1*sin(a)/u_s (line 46)
us_eqs = @(u_s) (y(4,i) - q_j(i))*S1*sin(a)*u_s + (S1*sin(a)).^3*((1/q_j(i).^3) - (1/y(4,i).^3))/(9*u_s) + y(4,i).^2 + 2*y(4,i)*sqrt(pi*y(2,i)) - 4*(y(4,i) - q_j(i))*S1/(q_j(i)*u_s) + (16*sin(a)*S1.^2*(y(4,i) - q_j(i)).^2/(3*q_j(i)*y(4,i)*Re)) + 4*sin(a).^3*S1.^4*(y(4,i) - q_j(i))*(y(4,i).^5-q_j(i).^5)/(15*q_j(i).^7*y(4,i).^5*Re*u_s.^2) - (y(4,i) - q_j(i))*S1*sin(a)/u_s;

Error in fsolve (line 270)
fuser = feval(funfcn{3},x,varargin{:});

Caused by:
Failure in initial objective function evaluation. FSOLVE cannot continue.

%% Test system if u_s is a constant
% u_s    = 1;
% y0     = [u_b0; s_b0; p0; R0];
% [t, y] = ode15s(@(t, y) odesym(t, y, We, Re, F, G, u_s), t, y0);
% for j = 1:4
%     subplot(2,2,j)
%     plot(t/pi, y(:,j)), grid on
%     xlabel('\theta/\pi ')
% end

%% Differential Equations
function dydt = odesym(t, y, We, Re, F, G, u_s)
    u_b  = y(1); 
    s_b  = y(2); 
    p    = y(3); 
    R    = y(4);
    F    = F(t);
    
    K    = (F^(3/2))/(G(t)^0.5);
    F_v  = (u_b - u_s*cos(p))/sqrt(s_b/pi)*K/(sin(p)*Re);
    
    dudt = (u_s^2*cos(p)*K - (u_b*u_s*K + F_v))/(u_b*s_b);
    dsdt = (2*u_b*u_s*K - (cos(p)*u_s^2*K - F_v))/(u_b^2);
    dpdt = (2*R/(We*sin(p)) - sin(p)*u_s^2*K)/(u_b^2*s_b) - 1;
    dRdt = R/tan(p);
    
    dydt = [dudt; dsdt; dpdt; dRdt];
end 

[Dyuman Joshi]

k1 = odesym(t(i,1),y(:,i),a,We,Re,F,u_s);

The function odesym is defined with 6 inputs variables, whereas you have provided 7 inputs in the above call. That is the cause of the error. 'a' is the extra input.

function dydt = odesym(t,y,We,Re,F,u_s)

[Surendra Ratnu]

Okay but after removing a i am getting new error "Array indices must be positive integers or logical values."

[Dyuman Joshi]

F is a scalar and you are trying to using t as an index to it in the function. Same with G, but, you have not passed G as a variable, yet you are using it in odesym().

There might be many more errors in your code, and rather than clearing them one-by-one, it will be better if you can describe what you are doing, so we have a better understanding of your code.

Provide the mathematical equations of the ODE you are trying to solve.

Also, you should incorporate putting comments in your code.

[John D'Errico]

I closed the duplicate question you posted. Note that really, you should be using a tool like ODE15i to solve problems like this, instead of writing your own code. Regardless, there is no need to keep on asking the same question.

[Surendra Ratnu]

@Dyuman Joshi I attached the equation system of ODE and algebraic which want to solve. Here u_b,s_b,phi, R and u_s are the dependent variable and θ is the independ variable and vary from 0 to 180 degree (θ = 0 to 180 degree). α, We, and Re are the constant ( α = 45 degree, We = 277, Re = 75). I want to plot the u_b,s_b,phi,u_s with θ. I tried using ODE15i but the solution is diverse after some value of theta. Any help is appreciate.

Answer 1

[Torsten]

As you can see, your algebraic equation does not seem to have a zero for your initial vector for the other variables.

Further tan(pi/2) = 0 so that you will get an immediate division-by-zero in the equation for R.

y0 = [0.0424*0.0042; 0.0424^2*0.0042; pi/2; 1.7440; 1.0];
us = 0:0.01:10;
for i = 1:numel(us)
  dydt = fun(0,[0.0424*0.0042; 0.0424^2*0.0042; pi/2; 1.7440;us(i)]);
  fus(i) = dydt(5);
end
plot(us,fus)

tspan = [0 pi];
M = eye(5);
M(5,5) = 0;
options = odeset('Mass',M);
%[T,Y] = ode15s(fun,tspan,y0,options)

function dydt = fun(t,y)
  alpha = deg2rad(45); 
  We = 277; 
  Re = 75; 

  ub = y(2)/y(1);
  sb = y(1)/ub;
  phi = y(3);
  R = y(4);
  us = y(5);

  q = -(0.3284-(3+0.671*sin(t)^2)^0.5)/(1-0.25*cos(t)^2);
  funF = @(Q)Q.*(-4.098*(0.4379+Q*cos(t)).^2+5.464*Q.^2*sin(t)^2+2.049);
  funG = @(Q)Q.*(-4.098*(0.4379+Q*cos(t)).^2+5.464*Q.^2*sin(t)^2+2.049).^3;
  F = integral(funF,0,q);
  G = integral(funG,0,q);
  hr = F^1.5/(G^0.5*R);
  Fnu = (ub-us*cos(phi))/sqrt(sb/pi) * hr*R/(sin(phi)*Re);

  dydt(1) = us*hr*R;
  dydt(2) = us^2*hr*R*cos(phi) + Fnu;
  dydt(3) = (2*R/(We*sin(phi))- us^2*hr*R*sin(phi))/(ub^2*sb) - 1;
  dydt(4) = R/tan(phi);
  dydt(5) = (R-q)*us*F*sin(alpha)+F^3*sin(alpha)^3/(9*us)*(1/q-1/R)+(R^2+2*R*sqrt(pi*sb))-...
            4*(R-q)*F/(q*us) + 16*sin(alpha)*F^2*(R-q)^2/(3*q^3*R*Re) + ...
            4*sin(alpha)^3*F^4*(R-q)*(R^5-q^5)/(15*us^2*q^7*R^5*Re) - (R-q)*sin(alpha)*F/us;
  dydt = dydt(:);
end

[Surendra Ratnu]

Yes, For initial vector u_s is 0.3475 and tha(pi/2) is infinity. When using ode15s, it gives error "Not enough input arguments." (Error in untitled05>fun (line 17)

v = y(1)) why??

[Torsten]

In the code above, replace

[T,Y] = ode15s(fun,tspan,y0,options)

by

[T,Y] = ode15s(@fun,tspan,y0,options)

But before the two problems (missing zero of the algebraic equation for the initial conditions and division by zero in the fourth equation) are not solved, you don't need to start computing. You will immediately get an error from ode15s.

And I don't see a line

v = y(1)

in the code for which you got the error message.

[Surendra Ratnu]

Sorry there is slight change in the expression of q, F and G. Now i did some change in code and also update the initial condition according to problem. I also attached the modified code. But I get solution upto theta = 0.0119 deg and getting warning " Warning: Failure at t=2.083673e-04. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (4.336809e-19) at t". Can anyone helping out for this problem.

[Torsten]

It seems that the reason for failure of the integrator is that ub becomes 0.

And here, you divide by ub:

dydt(3) = (2*R/(We*sin(phi))- us^2*hr*R*sin(phi))/(ub^2*sb) - 1;

clc
clear
close all

us = 1.4:0.001:10;
for i = 1:numel(us)
  dydt = fun(0,[0.1*0.0042; 0.1^2*0.6541; pi/2; 1.1965;us(i)]);
  fus(i) = dydt(5);
end
%plot(us,fus)

[~,idx] = min(abs(fus));
us0 = us(idx)

us0 = 1.5160

y0 = [0.1*0.0042; 0.1^2*0.6541; pi/2; 1.1965; us0];

%tspan = [0 pi];
tspan = [0 3e-4];
M = eye(5);
M(5,5) = 0;
options = odeset('Mass',M);
[T,Y] = ode15s(@fun,tspan,y0,options);

Warning: Failure at t=2.084712e-04. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (4.336809e-19) at time t.

UB = Y(:,2)./Y(:,1);
SB = Y(:,1)./UB;
PHI = Y(:,3);
R = Y(:,4);
US = Y(:,5);
figure(1)
plot(T,UB)
grid on

figure(2)
plot(T,SB)
grid on

figure(3)
plot(T,PHI)
grid on

figure(4)
plot(T,R)
grid on

figure(5)
plot(T,US)
grid on

function dydt = fun(t,y)
  alpha = deg2rad(45); We = 277; Re = 75;     
  ub = y(2)/y(1); sb = y(1)/ub; phi = y(3); R = y(4); us = y(5);

  q = (0.4630*cos(t)-(8-2*cos(t).^2-0.4287*sin(t).^2)^0.5)./(cos(t).^2-2);
  funF = @(Q)Q.*(1.0285-0.7408*cos(t)*Q - 0.8*Q.^2*(2-cos(t).^2));
  funG = @(Q)Q.*(1.0285-0.7408*cos(t)*Q - 0.8*Q.^2*(2-cos(t).^2)).^3;
  F = integral(funF,0,q);
  G = integral(funG,0,q);
  hr = F^1.5/(G^0.5*R);
  Fnu = (ub-us*cos(phi))/sqrt(sb/pi) * hr*R/(sin(phi)*Re);

  dydt(1) = us*hr*R;
  dydt(2) = us^2*hr*R*cos(phi) + Fnu;
  dydt(3) = (2*R/(We*sin(phi))- us^2*hr*R*sin(phi))/(ub^2*sb) - 1;
  dydt(4) = R/tan(phi);
  dydt(5) = (R-q)*us*F*sin(alpha)+F^3*sin(alpha)^3/(9*us)*(1/q-1/R)+(R^2+2*R*sqrt(pi*sb))-...
            4*(R-q)*F/(q*us) + 16*sin(alpha)*F^2*(R-q)^2/(3*q^3*R*Re) + ...
            4*sin(alpha)^3*F^4*(R-q)*(R^5-q^5)/(15*us^2*q^7*R^5*Re) - (R-q)*sin(alpha)*F/us;
  dydt = dydt(:);
end

[Surendra Ratnu]

Thank you for your effort. I observed the value of t start repeating after t = 2.0838e-04 and the value of ub is constant after t = 2.0838e-04. The value of equation u_s is not zero at us0. Is this issue come due to initial condition or non-zero value value of equation for initial condition or anything else ??

[Torsten]

The issue is due to your differential equations that make ub decrease and reach zero very fast. I don't know where you see that ub is constant after t = 2.0838e-04, I only see that the solver gives up at this time.

So you should check your equations and parameters to see if they are as you want them to be.

[Surendra Ratnu]

The differential equations are correct as per the problem. ub should increase with t as per the problem. I checked all equations and initial condition, but it gave the failure warning at some t. If i changed initial condition, solver gives up at different t. I am trying to solve this problem last few weeks. Can you help ??

Thanks

[Torsten]

As you can see from the result curve for ub, your equations result in a decreasing ub. So the error must lie in your equations or the parameters you use. I cannot help you in this respect.

Why do you prescribe

y0 = [0.1*0.0042; 0.1^2*0.6541; pi/2; 1.1965; us0];

This means that

ub * sb = 0.1*0.0042

and

ub^2 * sb = 0.1^2*0.6541

which gives a contradictory choice for sb (sb = 0.0042 in the first case, sb = 0.6541 in the second case).

[Surendra Ratnu]

Actual i change sb in my code. Is this issue occured due to solver or equations ?? Ok Thanks for your input.

[Torsten]

I didn't differentiate out the expressions for d(ub*sb)/dt and d(ub^2*sb)/dt, but solved in ub*sb and ub^2*sb instead of ub and sb. Thus at the first two positions of the y0 vector, you have to prescribe ub0*sb0 and ub0^2*sb0, not ub0 and sb0. But you should have noticed that when you looked over the code ...

[Sam Chak]

I executed the code from 'DAE_System.m,' and at some point, complex values emerged in 'ub'.

Is this behavior expected in the physical DAE System? The display is quite lengthy. I hope we can incorporate a 'hide and expand' or 'spoiler' feature in the forum.

Update: I removed the lengthy display of 'ub' values to make scrolling more user-friendly.

tspan   = [0 pi];
y0      = [0.1*0.0042; 0.1^2*0.6541; pi/2; 1.1965; 1.1423];
M       = eye(5);
M(5, 5) = 0;
options = odeset('Mass', M);
[T, Y]  = ode15s(@fun, tspan, y0, options);
plot(T, Y(:,2)./Y(:,1))

function dydt = fun(t,y)
    alpha   = deg2rad(45); 
    We      = 277; 
    Re      = 75;     
    ub      = y(2)/y(1);    % <-- view ub by removing the semicolon ';' 
    sb      = y(1)/ub; 
    phi     = y(3); 
    R       = y(4); 
    us      = y(5);
    
    q       = (0.4630*cos(t)-(8-2*cos(t).^2-0.4287*sin(t).^2)^0.5)./(cos(t).^2-2);
    funF    = @(Q)Q.*(1.0285-0.7408*cos(t)*Q - 0.8*Q.^2*(2-cos(t).^2));
    funG    = @(Q)Q.*(1.0285-0.7408*cos(t)*Q - 0.8*Q.^2*(2-cos(t).^2)).^3;
    F       = integral(funF,0,q);
    G       = integral(funG,0,q);
    hr      = F^1.5/(G^0.5*R);
    Fnu     = (ub-us*cos(phi))/sqrt(sb/pi) * hr*R/(sin(phi)*Re);
    
    dydt(1) = us*hr*R;
    dydt(2) = us^2*hr*R*cos(phi) + Fnu;
    dydt(3) = (2*R/(We*sin(phi))- us^2*hr*R*sin(phi))/(ub^2*sb) - 1;
    dydt(4) = R/tan(phi);
    dydt(5) = (R-q)*us*F*sin(alpha)+F^3*sin(alpha)^3/(9*us)*(1/q-1/R)+(R^2+2*R*sqrt(pi*sb))-...
            4*(R-q)*F/(q*us) + 16*sin(alpha)*F^2*(R-q)^2/(3*q^3*R*Re) + ...
            4*sin(alpha)^3*F^4*(R-q)*(R^5-q^5)/(15*us^2*q^7*R^5*Re) - (R-q)*sin(alpha)*F/us;
    dydt    = dydt(:);
end

[Sam Chak]

Hi @Surendra Ratnu

In the differential equation:

dydt(4) = R/tan(phi);

The initial value of phi is , thus making . During the integration process, phi approaches π causing . This leads to approaching a division-by-zero singularity event. What does this actually imply in the true physical system?

[Surendra Ratnu]

@Torsten @Sam Chak ub should not have any complex solution. It should increases with theta and phi is decreases with theta. The initial conditions are fixed for ub, sb, phi, R. us at t=0, i used 2.9701 (which i calculate using algebraic equation in matlab using solve function. I checked the equations which are correct. At t = pi, phi has value nearely 30 degree. So eq 4 it should not be tending to division by zero.I modified the initial condition in code and attached. Please help regarding this issues.

clc
clear
close all
us = 1.4:0.001:10;
for i = 1:numel(us)
  dydt = fun(0,[0.1*0.6541; 0.1^2*0.6541; pi/2; 1.1965;us(i)]);
  fus(i) = dydt(5);
end
%plot(us,fus)
[~,idx] = min(abs(fus));
us0 = us(idx)

us0 = 1.4000

y0 = [0.1*0.6541; 0.1^2*0.6541; pi/2; 1.1965; us0];
tspan = [0 pi];
M = eye(5);
M(5,5) = 0;
options = odeset('Mass',M);
[T,Y] = ode15s(@fun,tspan,y0,options);

Warning: Failure at t=3.146910e-02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.110223e-16) at time t.

UB = Y(:,2)./Y(:,1);
SB = Y(:,1)./UB;
PHI = Y(:,3);
R = Y(:,4);
US = Y(:,5);
figure(1)
plot(T,UB)
grid on

figure(2)
plot(T,SB)
grid on

figure(3)
plot(T,PHI)
grid on

figure(4)
plot(T,R)
grid on

figure(5)
plot(T,US)
grid on

function dydt = fun(t,y)
  alpha = deg2rad(45); We = 277; Re = 75;     
  ub = y(2)/y(1); sb = y(1)/ub; phi = y(3); R = y(4); us = y(5);
  q = (0.4630*cos(t)-(8-2*cos(t).^2-0.4287*sin(t).^2)^0.5)./(cos(t).^2-2);
  funF = @(Q)Q.*(1.0285-0.7408*cos(t)*Q - 0.8*Q.^2*(2-cos(t).^2));
  funG = @(Q)Q.*(1.0285-0.7408*cos(t)*Q - 0.8*Q.^2*(2-cos(t).^2)).^3;
  F = integral(funF,0,q);
  G = integral(funG,0,q);
  hr = F^1.5/(G^0.5*R);
  Fnu = (ub-us*cos(phi))/sqrt(sb/pi) * hr*R/(sin(phi)*Re);

  dydt(1) = us*hr*R;
  dydt(2) = us^2*hr*R*cos(phi) + Fnu;
  dydt(3) = (2*R/(We*sin(phi))- us^2*hr*R*sin(phi))/(ub^2*sb) - 1;
  dydt(4) = R/tan(phi);
  dydt(5) = (R-q)*us*F*sin(alpha)+F^3*sin(alpha)^3/(9*us)*(1/q-1/R)+(R^2+2*R*sqrt(pi*sb))-...
            4*(R-q)*F/(q*us) + 16*sin(alpha)*F^2*(R-q)^2/(3*q^3*R*Re) + ...
            4*sin(alpha)^3*F^4*(R-q)*(R^5-q^5)/(15*us^2*q^7*R^5*Re) - (R-q)*sin(alpha)*F/us;
  dydt = dydt(:);
end

[Torsten]

As you can see in the graphics for PHI, PHI reaches pi and the division by 0 in dR/dt = R/tan(phi) makes MATLAB give up.

[Surendra Ratnu]

Hi @Torsten I changed my code to solve these equations. First i solved for u_s then use u_s value in ode system and solve ode system by using ode45 in for loop. I tried ode15s but it gives wrong result. but while solving system i am getting error "Array indices must be positive integers or logical values". I also attched my code. Can you resolve my problem?

clc
clear
close all

a = deg2rad(45); m = 1; d_j = 0.6; 
We = 558; Re = 87; t0 = deg2rad(0);
b = (-0.13*m.^3+0.263*m.^2 +0.039*m +0.330)/tan(a);

%% Initial Conditions for differential equation
u_b0 = 0.1;  p0 = pi/2; R0 = 0.002*We/b; s_b0 = 0.01163;

%% Estimating Sheet Velocity at theta = 0 (initial sheet velocity)
x = b*cos(t0)*sin(a).^2; x01 = 4*((1-cos(t0).^2*cos(a).^2)); x02 = (b*sin(t0)*sin(a)).^2;
q_j0 = -((x -(x01-x02).^0.5)/(x01/4));
syms q 
d = (b.^2*sin(a).^2 + q.^2*(1-cos(t0).^2*cos(a).^2) + 2*b*q*cos(t0)*sin(a).^2).^0.5;
u_j0 = 2*(m-1)*(4*d.^2-1) + m;
F0 = int(q*u_j0, q, 0, q_j0);
F0 = double(F0);
eq_e = @(u_s) (R0-q_j0)*F0*sin(a)*u_s.^3 +...
              (R0.^2+2*R0*sqrt(pi*s_b0)+(16/(sin(a)*3*q_j0.^3*R0*Re))*((R0-q_j0)*F0*sin(a)).^2)*u_s.^2 -...
              ((R0-q_j0)*F0*sin(a)-((F0.^3*sin(a).^3)/(9))*((1/q_j0.^3)-(1/R0.^3))+((R0-q_j0)*F0*sin(a))*(4/(q_j0*sin(a))))*u_s-...
              ((R0-q_j0)*F0*sin(a))*((4*sin(a).^2*F0.^3*(R0.^5-q_j0.^5))/(15*q_j0.^7*R0.^5*Re));
u_s0 = fsolve(eq_e,0.5);

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

%% DAE System
syms q t
q_j = (-(b*cos(t)*sin(a)^2)+(1-cos(t)^2*cos(a)^2-(b*sin(t)*sin(a))^2)^0.5)/(1-cos(t)^2*cos(a)^2);
u_j = 2*(m-1)*(4*((q*sin(t))^2 + sin(a)^2*(q*cos(t)+b)^2)-1) + m;
F = matlabFunction(int(q*u_j, q, 0, q_j)); 
G = matlabFunction(int(q*u_j^3,q,0,q_j));
t = 0:0.01:pi;
y = zeros(5,length(t));
y(:,1)  = [u_b0; s_b0; p0; R0; u_s0];

for i = 1:length(t)-1  
    q_j = (-(b*cos(t(i))*sin(a)^2)+(1-cos(t(i))^2*cos(a)^2-(b*sin(t(i))*sin(a))^2)^0.5)/(1-cos(t(i))^2*cos(a)^2);
    S1 = double(F(t(i))) ;
    us_eqs = @(u_s)(y(4,i)-q_j(i))*S1*sin(a)*u_s + (S1*sin(a)).^3*((1/q_j(i).^3)-(1/y(4,i).^3))/(9*u_s)+...
             y(4,i).^2 + 2*y(4,i)*sqrt(pi*y(2,i)) - 4*(y(4,i)-q_j(i))*S1/(q_j(i)*u_s)+...
             (16*sin(a)*S1.^2*(y(4,i)-q_j(i)).^2/(3*q_j(i)*y(4,i)*Re)) + 4*sin(a).^3*S1.^4*(y(4,i)-q_j(i))*(y(4,i).^5-q_j(i).^5)/(15*q_j(i).^7*y(4,i).^5*Re*u_s.^2)-...
             (y(4,i)-q_j(i))*S1*sin(a)/u_s;
    u_s = fsolve(us_eqs, 1);
    sol = ode45(@(t,y) odesym(t,y,We,F,G,u_s),[0 pi],y(:,i));
end

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

Array indices must be positive integers or logical values.

Error in solution>odesym (line 50)
K = ((F^(3/2))/(G(t)^0.5));

Error in solution>@(t,y)odesym(t,y,We,F,G,u_s) (line 44)
sol = ode45(@(t,y) odesym(t,y,We,F,G,u_s),[0 pi],y(:,i));

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

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

function dydt = odesym(t,y,We,Re,F,G,u_s)
  u_b = y(1); s_b = y(2); p = y(3); R = y(4);
  F = F(t);
  K = ((F^(3/2))/(G(t)^0.5));
  F_v = (((u_b-u_s.*cos(p))/sqrt(s_b./pi)).*((K)./(sin(p).*Re)));
  dudt = ((u_s^2*cos(p)*K) - (u_b*u_s*K +F_v))./(u_b*s_b);
  dsdt = ((2*u_b*u_s*K) - (cos(p).*u_s^2*K-F_v))./(u_b^2);
  dpdt = (((2*R/(We*sin(p))) -sin(p)*u_s^2*K)./(u_b^2*s_b))-1;
  dRdt = R ./ tan(p);
  
  dydt = [dudt; dsdt; dpdt; dRdt];
end 

[Torsten]

Your function odesym has 7 input arguments

function dydt = odesym(t,y,We,Re,F,G,u_s)

but you let ode45 call it with 6:

sol = ode45(@(t,y) odesym(t,y,We,F,G,u_s),[0 pi],y(:,i));

Further in the last call, [0 pi] has to be replaced by [t(i) t(i+1)], and y(:,i) is not defined for i>1.

But all these changes to the code won't alter the problem that your equations produce a division by zero because PHI reaches pi. You must work on the equations themselves - the solution method is irrelevant in this case.

[Surendra Ratnu]

Actually i tried with solver ode15s. I am not getting PHI reaches to pi. I attched the graph. But this trend is not correct that why i want to first calculate the u_s value from algebraic equation and then solve the ode system using ode45. Can you tell me how to approach it ??

[Sam Chak]

Hi @Surendra Ratnu,

I have included the image of your differential-algebraic equations (DAE system) in your question. This way, when users view your post, they can better understand what you are trying to accomplish. Initially, I did not see the attached image, and your code contained 5 ODEs.

Now, the system has been reduced to only 4 ODEs, along with an algebraic equation that you are attempting to solve for . I also noticed that you made changes to the ODEs in the code. Please consider updating the ODEs in the image accordingly.

Based on your updated code, it seems you are trying to solve the algebraic equation at each time step and then pass it to the odesym() function for ode15s to solve the system. Am I correct?

[Surendra Ratnu]

Yes, i pass it to odesym() function for ode45 in updated code.

[Sam Chak]

Hi @Surendra Ratnu,

I appreciate your efforts in updating the code. However, we haven't had the chance to review the updated contents of the new ODEs, as requested earlier. It's possible there might be a slight oversight. Could you please take a moment to revisit my previous message and share the updated ODEs? Thank you for your attention.

[Surendra Ratnu]

i attched the an image for equation system

[Torsten]

Nothing has changed concerning the error.

p approaches pi, and you divide by tan(p) in dRdt. So ode15s gives up.

You won't establish a different result by just using a different solution method.

But I wonder why you use different equations to compute u_s:

First expression:

eq_e = @(u_s) (R0-q_j0)*F0*sin(a)*u_s.^3 +...
              (R0.^2+2*R0*sqrt(pi*s_b0)+(16/(sin(a)*3*q_j0.^3*R0*Re))*((R0-q_j0)*F0*sin(a)).^2)*u_s.^2 -...
              ((R0-q_j0)*F0*sin(a)-((F0.^3*sin(a).^3)/(9))*((1/q_j0.^3)-(1/R0.^3))+((R0-q_j0)*F0*sin(a))*(4/(q_j0*sin(a))))*u_s-...
              ((R0-q_j0)*F0*sin(a))*((4*sin(a).^2*F0.^3*(R0.^5-q_j0.^5))/(15*q_j0.^7*R0.^5*Re));

Second expression:

    us_eqs = @(u_s)(y(4,i)-q_j(i))*S1*sin(a)*u_s + (S1*sin(a)).^3*((1/q_j(i).^3)-(1/y(4,i).^3))/(9*u_s)+...
             y(4,i).^2 + 2*y(4,i)*sqrt(pi*y(2,i)) - 4*(y(4,i)-q_j(i))*S1/(q_j(i)*u_s)+...
             (16*sin(a)*S1.^2*(y(4,i)-q_j(i)).^2/(3*q_j(i)*y(4,i)*Re)) + 4*sin(a).^3*S1.^4*(y(4,i)-q_j(i))*(y(4,i).^5-q_j(i).^5)/(15*q_j(i).^7*y(4,i).^5*Re*u_s.^2)-...
             (y(4,i)-q_j(i))*S1*sin(a)/u_s;   

Which of the two is the correct one ?

clc
clear
close all
a = deg2rad(45); m = 1; d_j = 0.6; 
We = 558; Re = 87; t0 = deg2rad(0);
b = (-0.13*m.^3+0.263*m.^2 +0.039*m +0.330)/tan(a);
%% Initial Conditions for differential equation
u_b0 = 0.1;  p0 = pi/2; R0 = 0.002*We/b; s_b0 = 0.01163;
%% Estimating Sheet Velocity at theta = 0 (initial sheet velocity)
x = b*cos(t0)*sin(a)^2; x01 = 4*((1-cos(t0)^2*cos(a)^2)); x02 = (b*sin(t0)*sin(a))^2;
q_j0 = -((x -(x01-x02)^0.5)/(x01/4));
syms q 
d = (b.^2*sin(a).^2 + q.^2*(1-cos(t0).^2*cos(a).^2) + 2*b*q*cos(t0)*sin(a).^2).^0.5;
u_j0 = 2*(m-1)*(4*d.^2-1) + m;
F0 = int(q*u_j0, q, 0, q_j0);
F0 = double(F0);
eq_e = @(u_s) (R0-q_j0)*F0*sin(a)*u_s.^3 +...
              (R0.^2+2*R0*sqrt(pi*s_b0)+(16/(sin(a)*3*q_j0.^3*R0*Re))*((R0-q_j0)*F0*sin(a)).^2)*u_s.^2 -...
              ((R0-q_j0)*F0*sin(a)-((F0.^3*sin(a).^3)/(9))*((1/q_j0.^3)-(1/R0.^3))+((R0-q_j0)*F0*sin(a))*(4/(q_j0*sin(a))))*u_s-...
              ((R0-q_j0)*F0*sin(a))*((4*sin(a).^2*F0.^3*(R0.^5-q_j0.^5))/(15*q_j0.^7*R0.^5*Re));
u_s0 = fsolve(eq_e,0.5,optimset('Display','off'))

u_s0 = 6.3568e-06

%% DAE System
syms q t
q_j = (-(b*cos(t)*sin(a)^2)+(1-cos(t)^2*cos(a)^2-(b*sin(t)*sin(a))^2)^0.5)/(1-cos(t)^2*cos(a)^2);
u_j = 2*(m-1)*(4*((q*sin(t))^2 + sin(a)^2*(q*cos(t)+b)^2)-1) + m;
F = matlabFunction(int(q*u_j, q, 0, q_j)); 
G = matlabFunction(int(q*u_j^3,q,0,q_j));
y0 = [u_b0; s_b0; p0; R0; u_s0];
M = eye(5);
M(5,5) = 0;
options = odeset('Mass',M);
[T,Y] = ode15s(@(t,y) odesym(t,y,We,Re,a,b,F,G),[0 pi],y0,options);

Warning: Failure at t=1.022748e-02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.775558e-17) at time t.

UB = Y(:,2)./Y(:,1);
SB = Y(:,1)./UB;
PHI = Y(:,3);
R = Y(:,4);
US = Y(:,5);
figure(1)
plot(T,UB)
grid on

figure(2)
plot(T,SB)
grid on

figure(3)
plot(T,PHI)
grid on

figure(4)
plot(T,R)
grid on

figure(5)
plot(T,US)
grid on

function dydt = odesym(t,y,We,Re,a,b,F,G)
  u_b = y(1); s_b = y(2); p = y(3); R = y(4); u_s = y(5);
  q_j = (-(b*cos(t)*sin(a)^2)+(1-cos(t)^2*cos(a)^2-(b*sin(t)*sin(a))^2)^0.5)/(1-cos(t)^2*cos(a)^2);
  F = F(t);
  G = G(t);
  K = ((F^(3/2))/(G^0.5));
  F_v = (((u_b-u_s*cos(p))/sqrt(s_b/pi))*(K/(sin(p)*Re)));
  dudt = ((u_s^2*cos(p)*K) - (u_b*u_s*K +F_v))/(u_b*s_b);
  dsdt = ((2*u_b*u_s*K) - (cos(p).*u_s^2*K-F_v))/(u_b^2);
  dpdt = (((2*R/(We*sin(p))) - sin(p)*u_s^2*K)/(u_b^2*s_b))-1;
  dRdt = R / tan(p);
  dusdt = (R-q_j)*F*sin(a)*u_s + (F*sin(a))^3*((1/q_j^3)-(1/R^3))/(9*u_s)+...
          R^2 + 2*R*sqrt(pi*s_b) - 4*(R-q_j)*F/(q_j*u_s)+...
          (16*sin(a)*F^2*(R-q_j)^2/(3*q_j*R*Re)) + 4*sin(a)^3*F^4*(R-q_j)*(R^5-q_j^5)/(15*q_j^7*R^5*Re*u_s^2)-...
          (R-q_j)*F*sin(a)/u_s;
  dydt = [dudt; dsdt; dpdt; dRdt; dusdt];
end 

[Surendra Ratnu]

Both equations are same. one equation is rearranged.

[Sam Chak]

Hi @Surendra Ratnu

As a test, I've artfully infused a guiding line in 'phi' to thwart the inevitable descent of both and towards the abyss of zero. This mathematical spell permits me to glimpse into the future, unraveling the interplay of how the five states evolve across the expanse of .

Alas, despite my effort to prevent the order of the system from collapsing, defiantly converges to 0, setting ablaze the singularity in some divisional terms. The system gracefully succumbs at rad. Only you possess the power to quell the screaming of these equations and halt the impending madness.

a    = deg2rad(45); 
m    = 1; 
d_j  = 0.6; 
We   = 558; 
Re   = 87; 
t0   = deg2rad(0);
b    = (-0.13*m.^3+0.263*m.^2 +0.039*m +0.330)/tan(a);

%% Initial Conditions for differential equation
u_b0 = 0.1; 
p0   = pi/2; 
R0   = 0.002*We/b; 
s_b0 = 0.01163;

%% Estimating Sheet Velocity at theta = 0 (initial sheet velocity)
x    = b*cos(t0)*sin(a)^2; x01 = 4*((1-cos(t0)^2*cos(a)^2)); x02 = (b*sin(t0)*sin(a))^2;
q_j0 = -((x -(x01-x02)^0.5)/(x01/4));

syms q 
d    = (b.^2*sin(a).^2 + q.^2*(1-cos(t0).^2*cos(a).^2) + 2*b*q*cos(t0)*sin(a).^2).^0.5;
u_j0 = 2*(m-1)*(4*d.^2-1) + m;
F0   = int(q*u_j0, q, 0, q_j0);
F0   = double(F0);
eq_e = @(u_s) (R0-q_j0)*F0*sin(a)*u_s.^3 +...
              (R0.^2+2*R0*sqrt(pi*s_b0)+(16/(sin(a)*3*q_j0.^3*R0*Re))*((R0-q_j0)*F0*sin(a)).^2)*u_s.^2 -...
              ((R0-q_j0)*F0*sin(a)-((F0.^3*sin(a).^3)/(9))*((1/q_j0.^3)-(1/R0.^3))+((R0-q_j0)*F0*sin(a))*(4/(q_j0*sin(a))))*u_s-...
              ((R0-q_j0)*F0*sin(a))*((4*sin(a).^2*F0.^3*(R0.^5-q_j0.^5))/(15*q_j0.^7*R0.^5*Re));
u_s0 = fsolve(eq_e,0.5,optimset('Display','off'));

%% DAE System
syms q t
q_j     = (-(b*cos(t)*sin(a)^2)+(1-cos(t)^2*cos(a)^2-(b*sin(t)*sin(a))^2)^0.5)/(1-cos(t)^2*cos(a)^2);
u_j     = 2*(m-1)*(4*((q*sin(t))^2 + sin(a)^2*(q*cos(t)+b)^2)-1) + m;

F       = matlabFunction(int(q*u_j,   q, 0, q_j)); 
G       = matlabFunction(int(q*u_j^3, q, 0, q_j));

y0      = [u_b0; s_b0; p0; R0; u_s0];
M       = eye(5);
M(5,5)  = 0;
options = odeset('Mass', M);
[T, Y]  = ode15s(@(t,y) odesym(t,y,We,Re,a,b,F,G),[0 pi],y0,options);

Warning: Failure at t=1.523761e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552714e-15) at time t.

UB  = Y(:,2)./Y(:,1);
SB  = Y(:,1)./UB;
PHI = Y(:,3);
R   = Y(:,4);
US  = Y(:,5);

subplot(2,3,1), plot(T, UB),  grid on, title('u_{b}')
subplot(2,3,2), plot(T, SB),  grid on, title('s_{b}')
subplot(2,3,3), plot(T, PHI), grid on, title('\phi')
subplot(2,3,4), plot(T, R),   grid on, title('R')
subplot(2,3,5), plot(T, US),  grid on, title('u_{s}')

function dydt = odesym(t,y,We,Re,a,b,F,G)
    u_b   = y(1); 
    s_b   = y(2); 
    p     = y(3);
    R     = y(4); 
    u_s   = y(5);
    
    %% ----------
    % phi = p;                      % original phi
    phi   = pi/3*tanh(p) + pi/3;    % prevent sin(phi) and tan(phi) from going to 0
    %% ----------
    
    q_j   = (-(b*cos(t)*sin(a)^2)+(1-cos(t)^2*cos(a)^2-(b*sin(t)*sin(a))^2)^0.5)/(1-cos(t)^2*cos(a)^2); % 0.9122
    F     = F(t);                       % 0.4161
    G     = G(t);                       % 0.4161
    K     = ((F^(3/2))/(G^0.5));        % 0.4161
    F_v   = (((u_b - u_s*cos(p))/sqrt(s_b/pi))*(K/(sin(phi)*Re)));
    
    dudt  = ((u_s^2*cos(p)*K) - (u_b*u_s*K +F_v))/(u_b*s_b);
    dsdt  = ((2*u_b*u_s*K) - (cos(p).*u_s^2*K-F_v))/(u_b^2);
    dpdt  = (((2*R/(We*sin(phi))) - sin(p)*u_s^2*K)/(u_b^2*s_b))-1;    % singularity
    dRdt  = R/tan(phi);
    dusdt = (R-q_j)*F*sin(a)*u_s + (F*sin(a))^3*((1/q_j^3)-(1/R^3))/(9*u_s) + R^2 + 2*R*sqrt(pi*s_b) - 4*(R-q_j)*F/(q_j*u_s) + (16*sin(a)*F^2*(R-q_j)^2/(3*q_j*R*Re)) + 4*sin(a)^3*F^4*(R-q_j)*(R^5-q_j^5)/(15*q_j^7*R^5*Re*u_s^2) - (R-q_j)*F*sin(a)/u_s;
    dydt  = [dudt; dsdt; dpdt; dRdt; dusdt];
end 

[Surendra Ratnu]

but i tried with ode15s and not getting any failure from 0 to pi. please see attached code.

[Sam Chak]

@Surendra Ratnu, kudos! The system has danced its way triumphantly to the grand finale at .

run('DAE_System.m')

[Surendra Ratnu]

@Sam Chak It works well, but the trade of these parameter is not correct. The equations are correct but i didnot know why is gives this solution ??

[Sam Chak]

@Surendra Ratnu, with utmost respect, could you illuminate the scientific foundation behind your reservations about the accuracy of the plots, especially when you've previously affirmed the correctness of all equations with unwavering certainty?

[Surendra Ratnu]

Here phi should decreases with t, R should increases with t and increment is higher at large t, u should increases with t, s should increases with t.

[Sam Chak]

Hi @Surendra Ratnu

Your insights are valued, yet they appear more as subjective interpretations rather than scientifically substantiated claims. Delving deeper into the analysis of the DAE system through exploration of published technical papers would undoubtedly enrich our understanding.

I look forward to the potential unveiling of new knowledge or breakthroughs from your continued contributions.

[Torsten]

@Surendra Ratnu

Since one of your DAE_system.m files work (at least technically, I cannot interprete the results) and the other does not, there should be a difference in your equations, shouldn't it ?

And check again the two functions for u_s: they are not equal. I just checked the first term in both

(R0-q_j0)*F0*sin(a)*u_s.^3
(R-q_j)*F*sin(a)*u_s

and already found a discrepancy.