Error solving a neutral DDE

1 view (last 30 days)
Hossein
Hossein on 25 Nov 2021
Commented: Hossein on 27 Nov 2021
Hello,
I am trying to solve a set of neutral DDEs, using the code below
clc;clear;
S=.0725 ; vL =1.0e-3; c = 1481; pvtinf=0 ;
pinf0=1.01e5; k =7/5; rhoL=998;
d=10e-6;
R0=1e-6;
f=1e6;
ps=1e3;
X=[0 d];
tau1=[d/c,d/c,d/c,d/c];
tau2=[d/c,d/c,d/c,d/c];
y0=[R0,0,R0,0];
ff=@(t,x,xdel,xpdel)[x(2);
(((-0.5*x(2)^2)*(3-x(2)/(c))+(1+(1-3*k)*x(2)/c)*((pinf0-pvtinf)/rhoL+2*S/(rhoL*R0))*(R0/x(1))^(3*k)-2*S/(rhoL*x(1))-4*vL*x(2)/(rhoL*x(1))-(1+x(2)/c)*(pinf0-pvtinf+ps*subplus(sin(2*pi*f*t-2*pi*f*X(1)/c))/rhoL)-ps*x(1)*subplus(cos(2*pi*f*t-2*pi*f*X(1)/c)*2*pi*f/(rhoL*c)))-(2*xdel(3)*xdel(4)^2-xdel(3)*xpdel(4)^2)/d)/((1-x(2)/c)*x(1)+4*vL/(rhoL*c));
x(4);
(((-0.5*x(4)^2)*(3-x(4)/(c))+(1+(1-3*k)*x(4)/c)*((pinf0-pvtinf)/rhoL+2*S/(rhoL*R0))*(R0/x(3))^(3*k)-2*S/(rhoL*x(1))-4*vL*x(4)/(rhoL*x(3))-(1+x(4)/c)*(pinf0-pvtinf+ps*sin(2*pi*f*t-2*pi*f*X(2)/c))/rhoL-ps*x(3)*cos(2*pi*f*t-2*pi*f*X(4)/c)*2*pi*f/(rhoL*c))-(2*xdel(1)*xdel(2)^2-xdel(1)*xpdel(2)^2)/d)/((1-x(4)/c)*x(3)+4*vL/(rhoL*c));];
tspan =[0:0.01/f:20/f];
sol = ddensd(ff,tau1,tau2,y0,tspan);
Error using ddensd (line 236)
Index exceeds the number of array elements. Index must not exceed 2.
But it results in the following error which I am not sure where it is referring to. Any help is appreciated.

Accepted Answer

Walter Roberson
Walter Roberson on 25 Nov 2021
Last line of your dde, you have
(((-0.5*x(4)^2)*(3-x(4)/(c))+(1+(1-3*k)*x(4)/c)*((pinf0-pvtinf)/rhoL+2*S/(rhoL*R0))*(R0/x(3))^(3*k)-2*S/(rhoL*x(1))-4*vL*x(4)/(rhoL*x(3))-(1+x(4)/c)*(pinf0-pvtinf+ps*sin(2*pi*f*t-2*pi*f*X(2)/c))/rhoL-ps*x(3)*cos(2*pi*f*t-2*pi*f*X(4)/c)*2*pi*f/(rhoL*c))-(2*xdel(1)*xdel(2)^2-xdel(1)*xpdel(2)^2)/d)/((1-x(4)/c)*x(3)+4*vL/(rhoL*c));];
% ^^^^
At the location I marked, you have X(4) but X is only 1 x 2 .
I would recommend that you rename one of X or x to avoid this confusion.

More Answers (0)

Categories

Find more on Scope Variables and Generate Names in Help Center and File Exchange

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!