Solving the system of ODEs and algebraic Equation

Question

Surendra Ratnu 2023-12-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2058409-solving-the-system-of-odes-and-algebraic-equation

编辑： Torsten 2024-1-18

Equation_system.jpg

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 

5 个评论
显示 3更早的评论隐藏 3更早的评论

John D'Errico 2023-12-9

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 2023-12-10

编辑：Surendra Ratnu 2023-12-10

Picture2.jpg

@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 2023-12-10

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2058409-solving-the-system-of-odes-and-algebraic-equation#answer_1368844

编辑：Torsten 2023-12-10

在 MATLAB Online 中打开

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

31 个评论
显示 29更早的评论隐藏 29更早的评论

Torsten 2023-12-12

在 MATLAB Online 中打开

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

Sam Chak 2023-12-13

编辑：Sam Chak 2023-12-14

在 MATLAB Online 中打开

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

Surendra Ratnu 2023-12-14

编辑：Torsten 2023-12-14

在 MATLAB Online 中打开

DAE_System.m

@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

Surendra Ratnu 2024-1-16

编辑：Torsten 2024-1-16

在 MATLAB Online 中打开

DAE_System.m

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 2024-1-17

编辑：Torsten 2024-1-17

在 MATLAB Online 中打开

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

Sam Chak 2024-1-18

在 MATLAB Online 中打开

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

Sam Chak 2024-1-18

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 2024-1-18

编辑：Torsten 2024-1-18

在 MATLAB Online 中打开

@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.

请先登录，再进行评论。

Solving the system of ODEs and algebraic Equation

5 个评论
显示 3更早的评论隐藏 3更早的评论

回答（1 个）

31 个评论
显示 29更早的评论隐藏 29更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Solving the system of ODEs and algebraic Equation

5 个评论 显示 3更早的评论隐藏 3更早的评论

回答（1 个）

31 个评论 显示 29更早的评论隐藏 29更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

WeChat

5 个评论
显示 3更早的评论隐藏 3更早的评论

31 个评论
显示 29更早的评论隐藏 29更早的评论