Solve differential and algebraic equation system

I tried using the ode15s solver but i am getting this error message 'Index exceeds the number of array elements. Index must not exceed 4.'How can i solve this type of equation system numerically. Here i have 6 variable that depend on the t (in radian). I have equations for 5 variable but not for h_r. h_r is also depend on the t.

Please provide a solution to solve these equations.

Answer 1

[Torsten]

S_1 and q_j are symbolic functions. You have to convert them to numerical functions first (by "MatlabFunction") to use them within ode15s.

After having done this, you could try to determine u_s within the function "ode_vis" by

eq5 = @(u_s)((R-q_j)*u_s*S_1*sin(a)) + ((S_1.^3*sin(a).^3)/(9*u_s))*((1/q_j.^3)-(1/R.^3)) + (R.^ + 2*R*sqrt(pi*s_b)) -((4*S_1*(R-q_j))/(q_j.^2*u_s)) + (16*sin(a)*S_1.^2*(R-q_j).^2/(3*q_j.^3*R*Re)) + ((4*sin(a).^3*S_1.^4*(R-q_j)*(R.^5-q_j.^5))/(15*u_s*q_j.^7*R.^5*Re)) -( (sin(a)*S_1.^2*(R-q_j))/u_s);
u_s = fsolve(eq5,1)

at the beginning of the function instead of setting u_s = y(6).

Or you could add u_s to the solution variables by defining

dydt(6) = ((R-q_j)*u_s*S_1*sin(a)) + ((S_1.^3*sin(a).^3)/(9*u_s))*((1/q_j.^3)-(1/R.^3)) + (R.^ + 2*R*sqrt(pi*s_b)) -((4*S_1*(R-q_j))/(q_j.^2*u_s)) + (16*sin(a)*S_1.^2*(R-q_j).^2/(3*q_j.^3*R*Re)) + ((4*sin(a).^3*S_1.^4*(R-q_j)*(R.^5-q_j.^5))/(15*u_s*q_j.^7*R.^5*Re)) -( (sin(a)*S_1.^2*(R-q_j))/u_s);

and calling ode15s as

M = eye(6);
M(6,6) = 0;
y0=[0.1; 0.3024; 1.5725; 0.10; 0; 0];
[t,y] = ode15s(@(t,y) ode_vis1(t,y,We,Re), tspan, y0,optimset('Mass',M);

[William Rose]

@Surendra Ratnu, when @Torsten makes a recommendation, I do it, because @Torsten is very insightful and good at Matlab.

[Surendra Ratnu]

Thank you for your answer. I tried the code but it still not working for problem. But i am getting error to solve equation for u_s by using fsolve. Here is u_s is the function of R, s_b and t. I am not able to solve this. Error message is "Operator '-' is not supported for operands of type 'function_handle''. fsolve used only for numeric data.

Can i solve symbolic equation 5 for u_s in terms of R,s_b,t then i used the expression of u_s in differential equations ??

Can you provide me support to rectify the problem.

[Torsten]

Did you change S_1 and q_j to numerical function handles using "matlabFunction" ? I cannot see it because you only included "ode_vis", not the call to ode15s.

Can i solve symbolic equation 5 for u_s in terms of R,s_b,t then i used the expression of u_s in differential equations ??

Of course. But is it possible to explicitly solve for u_s ?

[Surendra Ratnu]

I changed but it gives error " Operator '-' is not supported for operands of type 'function_handle".

No, it not possible to explicitly solve for u_s. if i can differtiate eq5 with respect to t and get an equation with du_s/dt, ds_b/dt and dR/dt in the equation then it may work to solve the system of differtial equations. Can give me any hint to express eq5 in terms of differtial with u_s, s_b, R w.r.t t ??

[Torsten]

Include your actual code so that we can run it and see how far you got.

[Surendra Ratnu]

I attached my code. Can you provide me solution for this ??

clc
clear
close all

a = deg2rad(60); d_j = 0.6;u_0 = 2.63; st = 0.036; d_l = 997.3; vis = 0.0089; m = 0.636;u_s = 0.23;
We = 1e-03*d_l*u_0.^2*d_j/st; Re = d_l*u_0*d_j/vis;

t = [0 pi];

syms q t
b = (-0.13*m.^3+0.263*m.^2 +0.039*m +0.330)/tan(a);
x = 2*b*cos(t)*sin(a).^2; x_01 = 1-cos(t).^2.*cos(a).^2; x_03 = 4*b.^2*sin(t).^2*sin(a).^2;
q_j = -(((x)-(x_01 - x_03).^0.5)/x_01);
d = (b.^2*sin(a).^2 + q.^2*(1-cos(t).^2*cos(a).^2) + 2*b*q*cos(t)*sin(a).^2).^0.5;
u_j = 2*(m-1)*(4*d.^2-1) + m;
S_1 = int(q*u_j, q, 0, q_j);
S_1 = matlabFunction(S_1);

M = eye(6);
M(6,6) = 0;
y0=[0.1; 0.3024; 1.5725; 0.10; 0; 0];
[t,y] = ode15s(@(t,y) ode_vis(t,y,We,Re,m,a,S_1), [0 pi], y0, odeset('Mass',M));

Unrecognized function or variable 'dh_rdt'.

Error in solution>ode_vis (line 42)
dydt = [dudt; dsdt; dpdt; drdt;dh_rdt; du_sdt];

Error in solution>@(t,y)ode_vis(t,y,We,Re,m,a,S_1) (line 22)
[t,y] = ode15s(@(t,y) ode_vis(t,y,We,Re,m,a,S_1), [0 pi], y0, odeset('Mass',M));

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

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

%% Function 
function dydt = ode_vis(t,y,We,Re,m,a,S_1)
u_b = y(1); s_b = y(2); p = y(3); R = y(4); h_r = y(5); u_s = y(6);
b = (-0.13*m.^3+0.263*m.^2 +0.039*m +0.330)/tan(a);
x = 2*b*cos(t)*sin(a).^2; x_01 = 1-cos(t).^2.*cos(a).^2; x_03 = 4*b.^2*sin(t).^2*sin(a).^2;
q_j = -(((x)-(x_01 - x_03).^0.5)/x_01);
F_v = (((u_b-u_s*cos(p))/sqrt(s_b/pi))*((h_r*R)./(sin(p)*Re)));
dudt = (u_s^2*cos(p)*h_r*R - u_b*u_s*h_r*R +F_v)./(u_b*s_b);
dsdt = (2*u_b*u_s*h_r*R - cos(p)*u_s^2*h_r*R-F_v)./(u_b^2);
dpdt = (((2*R/(We*sin(p))) -sin(p)*u_s^2*h_r*R)./(u_b^2*s_b))-1;
drdt = R ./ tan(p);
du_sdt = ((R-q_j)*u_s*S_1(t)*sin(a)) + ((S_1(t).^3*sin(a).^3)/...
         (9*u_s))*((1/q_j.^3)-(1/R.^3)) + (R.^ + 2*R*sqrt(pi*s_b))...
         -((4*S_1(t)*(R-q_j))/(q_j.^2*u_s)) + ...
         (16*sin(a)*S_1(t).^2*(R-q_j).^2/(3*q_j.^3*R*Re)) + ...
         ((4*sin(a).^3*S_1(t).^4*(R-q_j)*(R.^5-q_j.^5))/(15*u_s*q_j.^7*R.^5*Re))...
         -( (sin(a)*S_1(t).^2*(R-q_j))/u_s);

dydt = [dudt; dsdt; dpdt; drdt;dh_rdt; du_sdt];
end

[Torsten]

I made changes to your code, but I cannot test it because you didn't define dh_rdt (see above).

[Surendra Ratnu]

Thanks for your reply. Here you use u_s as differtent differential equation but originally it is a algebraic equation. So Does it not change the solution? .

Here i defined dh_r/dt as a numerical equation because i dont have any seperate equation for dh_r/dt. So change in function is:-

delta_t = 1e-06;
h_r_next = h_r + delta_t;
dh_rdt = (h_r_next - h_r) / delta_t;

But i am getting error while solving and error is "Need a better guess y0 for consistent initial conditions``. and i can not change my initial conditions. They are fixed according to problem.

I dont know how much i correct to define the dh_r/dt like this.

I only have 4 differential and 1 algebraic equation. but does not have aby seperate equation for h_r. Can you give me some suggestion to solve this system of equation.

[Torsten]

Here you use u_s as differtent differential equation but originally it is a algebraic equation.

No, I use it as an algebraic equation because I set the corresponding row in the mass matrix to zero:

M = eye(6);
M(6,6) = 0;
[t,y] = ode15s(@(t,y) ode_vis(t,y,We,Re,m,a,S_1), [0 pi], y0, odeset('Mass',M));

Here i defined dh_r/dt as a numerical equation because i dont have any seperate equation for dh_r/dt. So change in function is:-

delta_t = 1e-06;
h_r_next = h_r + delta_t;
dh_rdt = (h_r_next - h_r) / delta_t;

That's nonsense since you get dh_rdt = (h_r + delta_t - h_r) / delta_t = 1.

You will have to search for a different way to get h_r.

[Surendra Ratnu]

Thank you for your suggestion. Can you tell me any other way to solve these equation??

[Torsten]

If you have no algebraic or differential equation for h_r, your system is incomplete and thus cannot be solved.

[Surendra Ratnu]

Thank you for your time @Torsten.

if i set the h_r value as constant (h_r = 0.30) and try to solve the equation then it gives error " Error using * " and "Error in daeic12 (line 63) F = UM' * f;". I checked all multiplication operator (* and changed to .*) but still this code given error. Can you give any hint to solve this ??

[Torsten]

clc
clear
close all

a = deg2rad(60); d_j = 0.6;u_0 = 2.63; st = 0.036; d_l = 997.3; vis = 0.0089; m = 0.636;u_s = 0.23;
We = 1e-03*d_l*u_0.^2*d_j/st; Re = d_l*u_0*d_j/vis;

syms q t
b = (-0.13*m.^3+0.263*m.^2 +0.039*m +0.330)/tan(a);
x = 2*b*cos(t)*sin(a).^2; x_01 = 1-cos(t).^2.*cos(a).^2; x_03 = 4*b.^2*sin(t).^2*sin(a).^2;
q_j = -(((x)-(x_01 - x_03).^0.5)/x_01);
d = (b.^2*sin(a).^2 + q.^2*(1-cos(t).^2*cos(a).^2) + 2*b*q*cos(t)*sin(a).^2).^0.5;
u_j = 2*(m-1)*(4*d.^2-1) + m;
S_1 = int(q*u_j, q, 0, q_j);
S_1 = matlabFunction(S_1);
q_j = matlabFunction(q_j);

y0=[0.1; 0.3024; 1.5725; 0.10];
tspan = [0 0.16];
iflag = 0;
[t,y] = ode15s(@(t,y) ode_vis(t,y,We,Re,m,a,S_1,iflag), tspan, y0);
for i = 1:numel(t)
  if i == 1
    [~,u_s(i)] = ode_vis(t(i),y(i,:),We,Re,m,a,S_1,1);
  else
    [~,u_s(i)] = ode_vis(t(i),y(i,:),We,Re,m,a,S_1,0); 
  end
end
figure(1)
plot(t,u_s)

figure(2)
plot(t,y(:,1))

figure(3)
plot(t,y(:,2))

figure(4)
plot(t,y(:,3))

figure(5)
plot(t,y(:,4))

figure(6)
plot(t,S_1(t))

figure(7)
plot(t,q_j(t))

%% Function 
function [dydt,u_s_set] = ode_vis(t,y,We,Re,m,a,S_1,iflag)

persistent u_s
if isempty(u_s) | iflag == 1
  u_s = 3.1931;
end

u_b = y(1); s_b = y(2); p = y(3); R = y(4); h_r = 0.3; 

b = (-0.13*m.^3+0.263*m.^2 +0.039*m +0.330)/tan(a);
x = 2*b*cos(t)*sin(a).^2; x_01 = 1-cos(t).^2.*cos(a).^2; x_03 = 4*b.^2*sin(t).^2*sin(a).^2;
q_j = -(((x)-(x_01 - x_03).^0.5)/x_01);
eq5 = @(u_s)((R-q_j)*u_s*S_1(t)*sin(a)) + ((S_1(t).^3*sin(a).^3)/...
         (9*u_s))*((1/q_j.^3)-(1/R.^3)) + (R.^ + 2*R*sqrt(pi*s_b))...
         -((4*S_1(t)*(R-q_j))/(q_j.^2*u_s)) + ...
         (16*sin(a)*S_1(t).^2*(R-q_j).^2/(3*q_j.^3*R*Re)) + ...
         ((4*sin(a).^3*S_1(t).^4*(R-q_j)*(R.^5-q_j.^5))/(15*u_s*q_j.^7*R.^5*Re))...
         -( (sin(a)*S_1(t).^2*(R-q_j))/u_s);
u_s = fsolve(eq5,u_s,optimset('Display','none'));
u_s_set = u_s;

F_v = (((u_b-u_s*cos(p))/sqrt(s_b/pi))*((h_r*R)./(sin(p)*Re)));
dudt = (u_s^2*cos(p)*h_r*R - u_b*u_s*h_r*R +F_v)./(u_b*s_b);
dsdt = (2*u_b*u_s*h_r*R - cos(p)*u_s^2*h_r*R-F_v)./(u_b^2);
dpdt = (((2*R/(We*sin(p))) -sin(p)*u_s^2*h_r*R)./(u_b^2*s_b))-1;
drdt = R ./ tan(p);


dydt = [dudt; dsdt; dpdt; drdt];
end

[Surendra Ratnu]

Can you explain why are you set the value of u_s = 0.23 in this code ?? Can you explain why is above mention error coming after running the code??

[Dyuman Joshi]

@Surendra Ratnu, you have yourself used u_s = 0.23 in your code in this comment - https://in.mathworks.com/matlabcentral/answers/2019736-solve-differential-and-algebraic-equation-system#comment_2885211

[Surendra Ratnu]

Soory @Dyuman Joshi that is by mistake. u_s is the function of R,s_b and t and an unknown variable.

[Dyuman Joshi]

Then change the definition of u_s accordingly.

[Torsten]

Can you explain why are you set the value of u_s = 0.23 in this code ??

The value u_s = 0.23 is not used. I use u_s = 3.1931 at the start, and this is the value I computed in advance and which satisfies your equation for u_s at t = 0.

Can you explain why is above mention error coming after running the code??

The error message from the integrator comes because it has difficulties integrating your equations further than t = 0.175. So check the results up to this time and see what could be wrong with your equations.

[Dyuman Joshi]

@Torsten, OP is talking about the 5th line of script, the last variable defined is u_s = 0.23.

[Torsten]

@Dyuman Joshi

I know, but this value is nowhere used in the code, I think. So the line could be deleted, and the result would be the same. I just didn't recognize it when I changed the code.

[Surendra Ratnu]

@Torsten I run code upto 0.175. I observed that the s_b v/s t is showing wrong result (s_b should increase with t but it decreses with t ) and u_s should change with t but it gives constant value. So u_s we can use a constant, it should vary with t. and u_s values increasing with t.

Can you tell why is this happening to solve.

[Torsten]

I included the plots for u_s, S_1, q_j and y(1)-y(4) up to 0.16 s. u_s changes with time, as you can see.

I cannot interprete the curves - it's up to you to find what is wrong with your equations. The solver just uses your inputs.

[Surendra Ratnu]

Is there any way to solve these equation without using fsolve??

[Torsten]

clc
clear
close all

a = deg2rad(60); d_j = 0.6;u_0 = 2.63; st = 0.036; d_l = 997.3; vis = 0.0089; m = 0.636;u_s = 0.23;
We = 1e-03*d_l*u_0.^2*d_j/st; Re = d_l*u_0*d_j/vis;

syms q t
b = (-0.13*m.^3+0.263*m.^2 +0.039*m +0.330)/tan(a);
x = 2*b*cos(t)*sin(a).^2; x_01 = 1-cos(t).^2.*cos(a).^2; x_03 = 4*b.^2*sin(t).^2*sin(a).^2;
q_j = -(((x)-(x_01 - x_03).^0.5)/x_01);
d = (b.^2*sin(a).^2 + q.^2*(1-cos(t).^2*cos(a).^2) + 2*b*q*cos(t)*sin(a).^2).^0.5;
u_j = 2*(m-1)*(4*d.^2-1) + m;
S_1 = int(q*u_j, q, 0, q_j);
S_1 = matlabFunction(S_1);
q_j = matlabFunction(q_j);

M = eye(5);
M(5,5) = 0;
y0=[0.1; 0.3024; 1.5725; 0.10; 3.1931];
tspan = [0 0.16];
[t,y] = ode15s(@(t,y) ode_vis(t,y,We,Re,m,a,S_1), tspan, y0, odeset('Mass',M));
figure(1)
plot(t,y(:,5))

figure(2)
plot(t,y(:,1))

figure(3)
plot(t,y(:,2))

figure(4)
plot(t,y(:,3))

figure(5)
plot(t,y(:,4))

figure(6)
plot(t,S_1(t))

figure(7)
plot(t,q_j(t))

%% Function 
function dydt = ode_vis(t,y,We,Re,m,a,S_1)
u_b = y(1); s_b = y(2); p = y(3); R = y(4); h_r = 0.3; u_s = y(5);
b = (-0.13*m.^3+0.263*m.^2 +0.039*m +0.330)/tan(a);
x = 2*b*cos(t)*sin(a).^2; x_01 = 1-cos(t).^2.*cos(a).^2; x_03 = 4*b.^2*sin(t).^2*sin(a).^2;
q_j = -(((x)-(x_01 - x_03).^0.5)/x_01);
F_v = (((u_b-u_s*cos(p))/sqrt(s_b/pi))*((h_r*R)./(sin(p)*Re)));
dudt = (u_s^2*cos(p)*h_r*R - u_b*u_s*h_r*R +F_v)./(u_b*s_b);
dsdt = (2*u_b*u_s*h_r*R - cos(p)*u_s^2*h_r*R-F_v)./(u_b^2);
dpdt = (((2*R/(We*sin(p))) -sin(p)*u_s^2*h_r*R)./(u_b^2*s_b))-1;
drdt = R ./ tan(p);
du_sdt = ((R-q_j)*u_s*S_1(t)*sin(a)) + ((S_1(t).^3*sin(a).^3)/...
         (9*u_s))*((1/q_j.^3)-(1/R.^3)) + (R.^ + 2*R*sqrt(pi*s_b))...
         -((4*S_1(t)*(R-q_j))/(q_j.^2*u_s)) + ...
         (16*sin(a)*S_1(t).^2*(R-q_j).^2/(3*q_j.^3*R*Re)) + ...
         ((4*sin(a).^3*S_1(t).^4*(R-q_j)*(R.^5-q_j.^5))/(15*u_s*q_j.^7*R.^5*Re))...
         -( (sin(a)*S_1(t).^2*(R-q_j))/u_s);

dydt = [dudt; dsdt; dpdt; drdt;du_sdt];
end

[Surendra Ratnu]

Thank you. Now i defined the h_r in the system of equation but i am getting error "Too many input arguments" in line du_sdt. I also tried to defined the G function outside function but i getting error and now i defined G inside the function still getting error in the line du_sdt. I attached my code. Please help me out for this error.

[Torsten]

For your conditions at t = 0, there is no value that satisfies the algebraic equation for u_s because the below plot doesn't cross the x-axis.

Please check the expression to calculate u_s. Especially the below term looks wrong:

(R^ + 2*R*sqrt(pi*s_b))

Use the below code to get a good approximation for y0(5). If this approximation is bad, you will get an error from the integrator.

a = deg2rad(45); m = 0.683; We = 277; Re = 75; 
b = (-0.13*m.^3+0.263*m.^2 +0.039*m +0.330)/tan(a);
y0=[0.10; 0.65; pi/2; 1.1036; 0.5022];

syms q 
t = 0;
x = b*cos(t)*sin(a)^2; x01 = 4*((1-cos(t)^2*cos(a)^2)); x02 = (b*sin(t)*sin(a))^2;
q_j = -((x -(x01-x02)^0.5)/(x01/4));
d = (b^2*sin(a)^2 + q^2*(1-cos(t)^2*cos(a)^2) + 2*b*q*cos(t)*sin(a)^2)^0.5;
u_j = 2*(m-1)*(4*d^2-1) + m;
s_1 = double(int(q*u_j, q, 0, q_j));
R = y0(4);
s_b = y0(2);

f = @(u_s)((R-q_j)*u_s*s_1*sin(a)) + ((s_1^3*sin(a).^3)/...
    (9*u_s))*((1/q_j^3)-(1/R^3)) + (R^ + 2*R*sqrt(pi*s_b))...
    -((4*s_1*(R-q_j))/(q_j^2*u_s)) + ...
    (16*sin(a)*s_1^2*(R-q_j)^2/(3*q_j^3*R*Re)) + ...
    ((4*sin(a)^3*s_1^4*(R-q_j)*(R^5-q_j^5))/(15*u_s*q_j^7*R^5*Re))...
    -( (sin(a)*s_1^2*(R-q_j))/u_s);
plot(0.1:0.01:7,arrayfun(@(u_s)f(u_s),0.1:0.01:7))

And leave the symbolic calculations and the conversion to MATLAB functions in the calling program. Just pass S_1 and G as function handles to the ode integrator:

a = deg2rad(45); m = 0.683; We = 277; Re = 75; 
b = (-0.13*m.^3+0.263*m.^2 +0.039*m +0.330)/tan(a);

syms q t
x = b*cos(t)*sin(a)^2; x01 = 4*((1-cos(t)^2*cos(a)^2)); x02 = (b*sin(t)*sin(a))^2;
q_j = -((x -(x01-x02)^0.5)/(x01/4));
Q_j = matlabFunction(q_j);
d = (b^2*sin(a)^2 + q^2*(1-cos(t)^2*cos(a)^2) + 2*b*q*cos(t)*sin(a)^2)^0.5;
u_j = 2*(m-1)*(4*d^2-1) + m;
S_1 = int(q*u_j, q, 0, q_j);
G = int(q*u_j^3,q,0,q_j);
S_1 = matlabFunction(S_1);
G = matlabFunction(G);
M = eye(5); M(5,5) = 0;
y0=[0.10; 0.65; pi/2; 1.1036; 0.5022]; tspan = [0 pi];
[t,y] = ode15s(@(t,y) ode_vis(t,y,We,Re,b,a,m,Q_j,S_1,G), tspan, y0, odeset('Mass',M));

Error using daeic12
Need a better guess y0 for consistent initial conditions.

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

figure(1)
plot(t,y(:,1),'r',t,y(:,2),'g',t,y(:,3),'k',t,y(:,4),'b',t,y(:,5),'m','LineWidth',1.5);

function dydt = ode_vis(t,y,We,Re,b,a,m,Q_j,S_1,G)
  u_b = y(1); s_b = y(2); p = y(3); R = y(4); u_s = y(5); 
  q_j = Q_j(t);
  s_1 = S_1(t);
  g = G(t);
  h_r = ((s_1^(3/2))/(g^0.5))/R;
  F_v = (((u_b-u_s*cos(p))/sqrt(s_b/pi))*((h_r*R)./(sin(p)*Re)));
  dudt = (u_s^2*cos(p)*h_r*R - u_b*u_s*h_r*R +F_v)./(u_b*s_b);
  dsdt = (2*u_b*u_s*h_r*R - cos(p)*u_s^2*h_r*R-F_v)./(u_b^2);
  dpdt = (((2*R/(We*sin(p))) -sin(p)*u_s^2*h_r*R)./(u_b^2*s_b))-1;
  drdt = R ./ tan(p);
  du_sdt = ((R-q_j)*u_s*s_1*sin(a)) + ((s_1^3*sin(a).^3)/...
           (9*u_s))*((1/q_j^3)-(1/R^3)) + (R^ + 2*R*sqrt(pi*s_b))...
           -((4*s_1*(R-q_j))/(q_j^2*u_s)) + ...
           (16*sin(a)*s_1^2*(R-q_j).^2/(3*q_j^3*R*Re)) + ...
           ((4*sin(a)^3*s_1^4*(R-q_j)*(R^5-q_j^5))/(15*u_s*q_j^7*R^5*Re))...
           -( (sin(a)*s_1^2*(R-q_j))/u_s);
  dydt = [dudt; dsdt; dpdt; drdt;du_sdt];
end

[Surendra Ratnu]

After doing some simplification i modified the code and i also checked all equations are correct but i could not able sol after 5.881630e-02. but i want to solve 0 to pi (t is in radian). Initial conditions are fixed for my case which is mention in code. Why is variable y(1),y(2) showing incorrect trade with t. Can help to figure out the this problem.

a = deg2rad(45); m = 0.683; We = 277; Re = 75; 
b = (-0.13*m.^3+0.263*m.^2 +0.039*m +0.330)/tan(a);

syms q t
x = b*cos(t)*sin(a)^2; x01 = 4*((1-cos(t)^2*cos(a)^2)); x02 = (b*sin(t)*sin(a))^2;
q_j = -((x -(x01-x02)^0.5)/(x01/4));
d = (b^2*sin(a)^2 + q^2*(1-cos(t)^2*cos(a)^2) + 2*b*q*cos(t)*sin(a)^2)^0.5;
u_j = 2*(m-1)*(4*d^2-1) + m;
S_1 = int(q*u_j, q, 0, q_j);
G = int(q*u_j^3,q,0,q_j);
S_1 = matlabFunction(S_1);
G = matlabFunction(G);

M = eye(5); M(5,5) = 0;
y0=[0.10; 0.65; pi/2; 1.2651; 0.2794]; tspan = [0 pi];
[t,y] = ode15s(@(t,y) ode_vis(t,y,We,Re,b,a,S_1,G), tspan, y0, odeset('Mass',M));

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

for i = 1:numel(t)
  dydt(:,i) = ode_vis(t(i),y(i,:),We,Re,b,a,S_1,G);
end
figure(1)
plot(t,y(:,1),'r',t,y(:,2),'g',t,y(:,3),'k',t,y(:,4),'b',t,y(:,5),'m','LineWidth',1.5);

figure(2)
plot(t,dydt(1,:),'r',t,dydt(2,:),'g',t,dydt(3,:),'k',t,dydt(4,:),'b',t,dydt(5,:),'m','LineWidth',1.5);

function dydt = ode_vis(t,y,We,Re,b,a,S_1,G)
  u_b = y(1); s_b = y(2); p = y(3); R = y(4); u_s = y(5); 
  x = b*cos(t)*sin(a)^2; x01 = 4*((1-cos(t)^2*cos(a)^2)); x02 = (b*sin(t)*sin(a))^2;
  q_j = -((x -(x01-x02)^0.5)/(x01/4));
  s_1 = S_1(t);
  g = G(t);
  h_r = ((s_1^(3/2))/(g^0.5))/R;
  F_v = (((u_b-u_s*cos(p))/sqrt(s_b/pi))*((h_r*R)/(sin(p)*Re)));
  dudt = (u_s^2*cos(p)*h_r*R - u_b*u_s*h_r*R +F_v)/(u_b*s_b);
  dsdt = (2*u_b*u_s*h_r*R - cos(p)*u_s^2*h_r*R-F_v)/(u_b^2);
  dpdt = (((2*R/(We*sin(p))) -sin(p)*u_s^2*h_r*R)/(u_b^2*s_b))-1;
  drdt = R / tan(p);
  du_sdt = (R-q_j)*s_1*sin(a)*u_s^3 +...
          (R.^2+2*R*sqrt(pi*s_b)+(16/(sin(a)*3*q_j^3*R*Re))*...
          ((R-q_j)*s_1*sin(a))^2)*u_s^2 -...
          ((R-q_j)*s_1*sin(a)-((s_1.^3*sin(a)^3)/(9))*((1/q_j^3)-...
          (1/R.^3))+((R-q_j)*s_1*sin(a))*(4/(q_j*sin(a))))*u_s-...
          ((R-q_j)*s_1*sin(a))*((4*sin(a)^2*s_1^3*(R^5-q_j^5))/...
          (15*q_j^7*R^5*Re));
  dydt = [dudt; dsdt; dpdt; drdt;du_sdt];
end

[Torsten]

Why is variable y(1),y(2) showing incorrect trade with t. Can help to figure out the this problem.

I can find MATLAB syntax errors in your code, but how should I be able to help with the model equations ?

[Surendra Ratnu]

I checked my code. It seems correct. Can you just tell me why is it not solve equation upto pi and Is there any other solver in matlab to solve this type of the system of equation ??

[Torsten]

I added the plot of the derivatives that become very large at t = 5.88e-2. There could be a singular point in the solution of your equations (like 1/x at x = 0) (maybe because a denominator becomes 0 in the dydt expressions) or a simple programming error.

MATLAB's ode integrators are made to integrate odes - thus there is no other MATLAB software suited to integrate your equations.

[Surendra Ratnu]

ok @Torsten thanks for your valable suggestion.

[Torsten]

I just noticed that your equation for F_v has sin(p) in the denominator. As you can see from the plot of p = y(3), p reaches pi at t = 5.88e-2. Thus F_v becomes undefined since sin(pi) = 0.

Same for dpdt.

[Surendra Ratnu]

yes @Torsten i checked it. Can you tell me there is any other way to solve these equation (means without using odes and somthing like R-K method and etc)??

[Torsten]

You solved the equations of your model. The solution for your model as it stands explodes at t = 5.88e-2 because you divide by sin(p) which becomes 0. So don't search for new numerical methods because all will give you this solution - search for a modification of your model.

Answer 2

[William Rose]

@Surendra Ratnu,

[edit: added several comments]

Your code includes

y0 = [0.1; 0.3024; 1.5725; 0.10];
t = [0 pi];
[t,y] = ode15s(@ode_vis1,t,y0);

I renamed the function from ode_vis to ode_vis1 so that the function name would not conflict with the main file name.

y0 is 4x1, but ode_vis1() returns dydt which is 6x1. y0 and dydt must be the same size.

You can replace line 1 above as shown:

y0=[0.1; 0.3024; 1.5725; 0.10; 0; 0];

This fixes the "Index exceeds number of array elements..." error.

You must replace lines 2 and 3 above with

tspan = [0 pi];
[t,y] = ode15s(@ode_vis1,tspan,y0);

to avoid a conflict involving t.

Since you pass parameters We and Re to ode_vis1(), the call to ode15s() should look like this:

[t,y] = ode15s(@(t,y) ode_vis1(t,y,We,Re), tspan, y0);

The body of function dydt=ode_vis1() includes variables a and m on the right hand side of some equations. Variables a and m are not defined in the function, and they are not passed to the function, and they are not global. You will need to fix this. You can fix it by passing a and m to ode_vis1(), as shown:

function dydt = ode_vis1(t,y, We,Re,m,a)

Then you need to modify the call to ode15s:

[t,y] = ode15s(@(t,y) ode_vis1(t,y,We,Re,m,a), tspan, y0);

There is still an error because ode_vis1() has the variable q on the right hand side, and q is not provided. You need to provide it.

The last element (6th element) of dydt does not appear to be a derivative. You have

dydt = [dudt; dsdt; dpdt; drdt; dhdt; eq5];

The variable eq5 is defined by a logical equality:

eq5 = f(x) == g(x);

where f(x) and g(x) are functions involving many variables. Therefore dydt(6) will be either 1 or 0 each time ode_vis1() is called.

In function dydt=ode_vis1(), you define

u_s = y(6);

Given that definition, then the 6th element of dydt should be the time derivative of u_6. If eq5 is not the derivative of u_6, then you need to fix something.

If eq5 is your attempt to enforce an algebraic equation, then read carefully the Matlab help for differential algebraic equations. Do the example in the help for ode15s involving the Robertson problem. Be sure you understand what a mass matrix is and how to use it. Do the example involving a mass matrix for a baton thrown into the air to gain familiarity with mass matrices.

Good luck with your work.

[Surendra Ratnu]

u_s is simply algebraic equation. In the first four differential equation, u_s is in equations. u_s is depend on the u_s So 4 differential equation and one algebraic equation want to solve for variable u_b,s_b,p,R,h_r and u_s for the t numerically but i do not know how to deal with that. Here i do not have any equation for h_r. Other than all mention variable, have constant value. So there is any other way to solve these equation numerically ?? if we can defined a difference equation for h_r then solve numerically.

Now i modified the code. Now it did not showing error for q variable. But i do not know how to deal with h_r.

I also tried with DAE solver but wont help for me much.

Please provide me a way to handle this ??

Thank You

[William Rose]

@Torsten has once again gone above and beyond in his effort to assist.