Can anyone help me fix that error please

function dydt = vdp1(t,y)
% input parameters
Dr = 0.05;
Dc = 0.05;
Dh = 0.05;
Dt = 0.09;
Db = 0.099;
Lr = 0.5;
Lc = 0.5;
Lh = 0.5;
Li = 0.5;
g = 9.81;
rho = 1000;
mu = 1.002e-3;
Lp = 1;
Dp = 0.0414;
H0 = 1.5;
Hs = 0.025;%
C0 = 10;%
Cs = 10;%
Ct = 10;
Pa = 10^5;
D0 = 0.032;
Ds = 0.025;%
Kv0 = 2;%
Kvs = 2;%
Kp = 2992;
mI2 = 0.3;%
Tr = 30;%
Te = 100;
Tc = 20;
%%
Vr = (pi*(Dr/2)^2)*Lr;
Ac = pi*(Dc/2)^2;
Ah = pi*(Dh/2)^2;
Ap = pi*(Dp/2)^2;
A0 = pi*(D0/2)^2;
As = pi*(Ds/2)^2;
At = pi*(Dt/2)^2;
Ab = pi*(Db/2)^2;
%%
P = y(1);
lh = y(2);
dlh = y(3);
lc = y(4);
dlc = y(5);
lt = y(6);
dlt = y(7);





% Pm equation 59
Pm = ((Ab^(2)/(mI2+mf))+(At/(rho*(lt+Ct)))+(Ac/(rho*lc))+(Ah/(rho*lh)))^(-1)...
    *((-Ab/(mI2+mf))*((rho*g*(Li+xp)-Pa)*Ab + Cf*dxp+Hf+Cfd*dxp^(2)-Kp*xp)-...
    (At/(rho*lt))*(-Pa-rho*g*lt-8*mu*pi*(lt-Ct)*dlt/At)-(Ac/(rho*lc))*(-P-rho*g*lc-8*mu*pi*lc*dlc/Ac)...
    -(Ah/(rho*lh))*(P-rho*g*lh-8*pi*mu*lh*dlh/Ah));
%%
dydt = zeros(9,1);

Vc = Ac*(Lc - lc);
Ve = Ah*(Lh - lh);

dydt(1) = (-P/((Vr/Tr)+(Vc/Tc)+(Ve/Te)))*((1/Tc)*(-Ac*dlc)+(1/Te)*(-Ah*dlh));
dydt(2) = y(3);
dydt(3) = (1/(rho*lh))*((Pm - P) - rho*g*lh - 8*pi*mu*(lh/Ah)*dlh);
dydt(4) = y(5);
dydt(5) = (1/(rho*lc))*((Pm-P)-rho*g*lc-8*pi*mu*lc*dlc/Ac);
dydt(6) = y(7);
dydt(7) = (1/(rho*(lt+Ct)))*((Pm-Pa)-rho*g*lt-8*mu*pi*(lt+Ct)*dlt/At);

end


close all;
clear ;
tspan = 0:0.01:20;       
y0 = [1.1025e05 0.5 0 0.4 0 1.5 0 0.04 0]' ;
[t,y] = ode45(@vdp1, tspan, y0);
%%
figure(1)
plot(tspan,y(:,1),'b','LineWidth',2)
title('Variations of working gas pressure')
xlabel('time (Sec)')
ylabel('P (pa)')
%%
figure(2)
plot(tspan,y(:,2),'b','LineWidth',2)
title('Hot liquid piston oscillations')
xlabel('time (Sec)')
ylabel('lh')
%%
figure(3)
plot(tspan,y(:,4),'b','LineWidth',2)
title('Cold liquid piston oscillations')
xlabel('time (Sec)')
ylabel('lc')
%%
figure(4)
plot(tspan,y(:,6),'b','LineWidth',2)
title('Tuning column oscillations')
xlabel('time (Sec)')
ylabel('lc')

Answer 1

[Stephan]

You missed to define several variables, which are used in your equations:

% several used variables are not defined! - i set them all = 1:
mf = 1;
xp = 1;
Cf = 1;
dxp  = 1;
Hf = 1;
Cfd = 1;

The complete and working code is here - just change the values for the variables as needed:

tspan = 0:0.01:20;       
y0 = [1.1025e05 0.5 0 0.4 0 1.5 0 0.04 0]' ;
[t,y] = ode45(@vdp1, tspan, y0);
%%
figure(1)
plot(tspan,y(:,1),'b','LineWidth',2)
title('Variations of working gas pressure')
xlabel('time (Sec)')
ylabel('P (pa)')
%%
figure(2)
plot(tspan,y(:,2),'b','LineWidth',2)
title('Hot liquid piston oscillations')
xlabel('time (Sec)')
ylabel('lh')
%%
figure(3)
plot(tspan,y(:,4),'b','LineWidth',2)
title('Cold liquid piston oscillations')
xlabel('time (Sec)')
ylabel('lc')
%%
figure(4)
plot(tspan,y(:,6),'b','LineWidth',2)
title('Tuning column oscillations')
xlabel('time (Sec)')
ylabel('lc')

function dydt = vdp1(~,y)
% several used variables are not defined! - i set them all = 1:
mf = 1;
xp = 1;
Cf = 1;
dxp  = 1;
Hf = 1;
Cfd = 1;

% input parameters
Dr = 0.05;
Dc = 0.05;
Dh = 0.05;
Dt = 0.09;
Db = 0.099;
Lr = 0.5;
Lc = 0.5;
Lh = 0.5;
Li = 0.5;
g = 9.81;
rho = 1000;
mu = 1.002e-3;
Lp = 1;
Dp = 0.0414;
H0 = 1.5;
Hs = 0.025;%
C0 = 10;%
Cs = 10;%
Ct = 10;
Pa = 10^5;
D0 = 0.032;
Ds = 0.025;%
Kv0 = 2;%
Kvs = 2;%
Kp = 2992;
mI2 = 0.3;%
Tr = 30;%
Te = 100;
Tc = 20;
%%
Vr = (pi*(Dr/2)^2)*Lr;
Ac = pi*(Dc/2)^2;
Ah = pi*(Dh/2)^2;
Ap = pi*(Dp/2)^2;
A0 = pi*(D0/2)^2;
As = pi*(Ds/2)^2;
At = pi*(Dt/2)^2;
Ab = pi*(Db/2)^2;
%%
P = y(1);
lh = y(2);
dlh = y(3);
lc = y(4);
dlc = y(5);
lt = y(6);
dlt = y(7);
% Pm equation 59
Pm = ((Ab^(2)/(mI2+mf))+(At/(rho*(lt+Ct)))+(Ac/(rho*lc))+(Ah/(rho*lh)))^(-1)...
    *((-Ab/(mI2+mf))*((rho*g*(Li+xp)-Pa)*Ab + Cf*dxp+Hf+Cfd*dxp^(2)-Kp*xp)-...
    (At/(rho*lt))*(-Pa-rho*g*lt-8*mu*pi*(lt-Ct)*dlt/At)-(Ac/(rho*lc))*(-P-rho*g*lc-8*mu*pi*lc*dlc/Ac)...
    -(Ah/(rho*lh))*(P-rho*g*lh-8*pi*mu*lh*dlh/Ah));
%%
dydt = zeros(9,1);
Vc = Ac*(Lc - lc);
Ve = Ah*(Lh - lh);
dydt(1) = (-P/((Vr/Tr)+(Vc/Tc)+(Ve/Te)))*((1/Tc)*(-Ac*dlc)+(1/Te)*(-Ah*dlh));
dydt(2) = y(3);
dydt(3) = (1/(rho*lh))*((Pm - P) - rho*g*lh - 8*pi*mu*(lh/Ah)*dlh);
dydt(4) = y(5);
dydt(5) = (1/(rho*lc))*((Pm-P)-rho*g*lc-8*pi*mu*lc*dlc/Ac);
dydt(6) = y(7);
dydt(7) = (1/(rho*(lt+Ct)))*((Pm-Pa)-rho*g*lt-8*mu*pi*(lt+Ct)*dlt/At);
end

[Mohamed El kholy]

Thanks Mr. Stephen, I treid now to run it but it still gave me that error.

>> solve

Error using ode45

Too many input arguments.

Error in solve (line 3)

[t,y] = ode45(@vdp1, tspan, y0);

If you run it, may you please send me the graphs you got, it is exptected to give me stable sign waves occialtions.

[Mohamed El kholy]

Also, please check this corrected one after remving the variables that not used.

Your help is greatly apprecaietd regadring that matter !

tspan = 0:0.01:20;       
y0 = [1.1025e05 0.5 0 0.4 0 1.5 0 0.04 0]' ;
[t,y] = ode45(@vdp1, tspan, y0);
%%
figure(1)
plot(tspan,y(:,1),'b','LineWidth',2)
title('Variations of working gas pressure')
xlabel('time (Sec)')
ylabel('P (pa)')
%%
figure(2)
plot(tspan,y(:,2),'b','LineWidth',2)
title('Hot liquid piston oscillations')
xlabel('time (Sec)')
ylabel('lh')
%%
figure(3)
plot(tspan,y(:,4),'b','LineWidth',2)
title('Cold liquid piston oscillations')
xlabel('time (Sec)')
ylabel('lc')
%%
figure(4)
plot(tspan,y(:,6),'b','LineWidth',2)
title('Tuning column oscillations')
xlabel('time (Sec)')
ylabel('lc')
function dydt = vdp1(~,y)

% input parameters
Dr = 0.05;
Dc = 0.05;
Dh = 0.05;
Dt = 0.09;
Lr = 0.5;
Lc = 0.5;
Lh = 0.5;
g = 9.81;
rho = 1000;
mu = 1.002e-3;
Pa = 10^5;
Tr = 30;%
Te = 100;
Tc = 20;
%%
Vr = (pi*(Dr/2)^2)*Lr;
Ac = pi*(Dc/2)^2;
Ah = pi*(Dh/2)^2;

At = pi*(Dt/2)^2;
%%
P = y(1);
lh = y(2);
dlh = y(3);
lc = y(4);
dlc = y(5);
lt = y(6);
dlt = y(7);
% Pm equation 59
Pm = (-(At/(rho*lt))*(-Pa-rho*g*lt-8*mu*pi*(lt-Ct)*dlt/At)-(Ac/(rho*lc))*(-P-rho*g*lc-8*mu*pi*lc*dlc/Ac)...
    -(Ah/(rho*lh))*(P-rho*g*lh-8*pi*mu*lh*dlh/Ah));
%%
dydt = zeros(9,1);
Vc = Ac*(Lc - lc);
Ve = Ah*(Lh - lh);
dydt(1) = (-P/((Vr/Tr)+(Vc/Tc)+(Ve/Te)))*((1/Tc)*(-Ac*dlc)+(1/Te)*(-Ah*dlh));
dydt(2) = y(3);
dydt(3) = (1/(rho*lh))*((Pm - P) - rho*g*lh - 8*pi*mu*(lh/Ah)*dlh);
dydt(4) = y(5);
dydt(5) = (1/(rho*lc))*((Pm-P)-rho*g*lc-8*pi*mu*lc*dlc/Ac);
dydt(6) = y(7);
dydt(7) = (1/(rho*(lt+Ct)))*((Pm-Pa)-rho*g*lt-8*mu*pi*(lt+Ct)*dlt/At);
end

[Stephan]

Your system appears to be highly stiff. I rechanged the formula for Pm to the long version, because i was not able to get a solution with the one you used. Note that this means you need parameter values for the additional parameters. Maybe the problem is more easy to solve when meaningful parameters are used. Also i had to reduce the 'RelTol' Parameter extremly and used a stiff solver. Also note that i have changed tspan to a small value - have a look at what happens when you set the upper limit to 20. Since i have no insight in your problem, i guess there is not much more i can do for you.

tspan = [0 0.5];       
y0 = [1.1025e05, 0.5, 0, 0.4, 0, 1.5, 0, 0.04, 0]';
opts = odeset('RelTol',1e-2);
[t,y] = ode15s(@vdp1, tspan, y0, opts);
%%
figure(1)
plot(t,y(:,1),'b','LineWidth',2)
title('Variations of working gas pressure')
xlabel('time (Sec)')
ylabel('P (pa)')
%%
figure(2)
plot(t,y(:,2),'b','LineWidth',2)
title('Hot liquid piston oscillations')
xlabel('time (Sec)')
ylabel('lh')
%%
figure(3)
plot(t,y(:,4),'b','LineWidth',2)
title('Cold liquid piston oscillations')
xlabel('time (Sec)')
ylabel('lc')
%%
figure(4)
plot(t,y(:,6),'b','LineWidth',2)
title('Tuning column oscillations')
xlabel('time (Sec)')
ylabel('lc')

function dydt = vdp1(~,y)
% Missing parameters:
mf = 1;

% input parameters
Dr = 0.05;
Dc = 0.05;
Dh = 0.05;
Dt = 0.09;
Lr = 0.5;
Lc = 0.5;
Lh = 0.5;
g = 9.81;
rho = 1000;
mu = 1.002e-3;
Pa = 10^5;
Tr = 30;%
Te = 100;
Tc = 20;
%%
% Dp = 0.0414;
% D0 = 0.032;
% Ds = 0.025;%
Ct = 10;
Db = 0.099;
mI2 = 0.3;%
Li = 0.5;
xp = 1;
Cf = 1;
dxp  = 1;
Hf = 1;
Cfd = 1;
Kp = 2992;

Vr = (pi*(Dr/2)^2)*Lr;
Ac = pi*(Dc/2)^2;
Ah = pi*(Dh/2)^2;
% Ap = pi*(Dp/2)^2;
% A0 = pi*(D0/2)^2;
% As = pi*(Ds/2)^2;
% At = pi*(Dt/2)^2;
Ab = pi*(Db/2)^2;
At = pi*(Dt/2)^2;


%%
P = y(1);
lh = y(2);
dlh = y(3);
lc = y(4);
dlc = y(5);
lt = y(6);
dlt = y(7);
% Pm equation 59
Pm = ((Ab^(2)/(mI2+mf))+(At/(rho*(lt+Ct)))+(Ac/(rho*lc))+(Ah/(rho*lh)))^(-1)...
    *((-Ab/(mI2+mf))*((rho*g*(Li+xp)-Pa)*Ab + Cf*dxp+Hf+Cfd*dxp^(2)-Kp*xp)-...
    (At/(rho*lt))*(-Pa-rho*g*lt-8*mu*pi*(lt-Ct)*dlt/At)-(Ac/(rho*lc))*(-P-rho*g*lc-8*mu*pi*lc*dlc/Ac)...
    -(Ah/(rho*lh))*(P-rho*g*lh-8*pi*mu*lh*dlh/Ah));
%%
dydt = zeros(9,1);
Vc = Ac*(Lc - lc);
Ve = Ah*(Lh - lh);
dydt(1) = (-P/((Vr/Tr)+(Vc/Tc)+(Ve/Te)))*((1/Tc)*(-Ac*dlc)+(1/Te)*(-Ah*dlh));
dydt(2) = y(3);
dydt(3) = (1/(rho*lh))*((Pm - P) - rho*g*lh - 8*pi*mu*(lh/Ah)*dlh);
dydt(4) = y(5);
dydt(5) = (1/(rho*lc))*((Pm-P)-rho*g*lc-8*pi*mu*lc*dlc/Ac);
dydt(6) = y(7);
dydt(7) = (1/(rho*(lt+Ct)))*((Pm-Pa)-rho*g*lt-8*mu*pi*(lt+Ct)*dlt/At);
end

Answer 2

[Mohamed El kholy]

Thanks so much for your help. Highly appreciated!

I will run it and see.