What is the meaning of "Index in position 1 exceeds array bounds (must not exceed 1)".

I was trying to solve an integrodifferential equation using MATLAB. I decided to use symbolic to solve ODE, so I wrote the following code.

%% EQUATIONS GOVERNING THE BUBBLE DYNAMICS
% Yihan Zhang
clc;clear;close all

syms t r R(t)

%% parameters
tau_y=8.6; % Pa
gamma=0.07; % Pa m
rho=1000; % kg/m^3
G=81.5; % Pa
P0=1.13e5; % Pa
deltaP=100; % Pa
eta_s=0.001; % Pa s
R0=100e-4; % m（boundary condition) R(t=0)
Rp0=0; % m/s (boundary condition) dR/dt(t=0)
w=((3*P0+4*gamma/R0+4*G)/(rho*R0^2))^0.5; % s^-1

%% expressions of other variables
tau_rr=-4*G.*(R(t).^3-R0^3)./(3.*r.^3);
tau_tt=-tau_rr./2;
Pg=(P0+2*gamma/R0).*(R0./R(t)).^3;
P_inf=P0+deltaP.*sin(w.*t);
P=Pg+tau_rr-eta_s.*4.*R(t)./R(t)-2.*gamma./R(t);

%% differential equation
eqn= [1000*(diff(R(t),t,2)*R(t)+3/2*diff(R(t),t)^2)==P-P_inf-tau_rr+2*int((tau_rr-tau_tt)/r,R(t),9999)]
[eqs,vars]=reduceDifferentialOrder(eqn,R(t))
[M,F]=massMatrixForm(eqs,vars)
M=odeFunction(M,vars)
Fun=odeFunction(F,vars,r)
H=@(t,R) Fun(t,R)
opt=odeset('mass',M,'InitialSlope',Rp0);
ode15s(H,[0,0.0001],R0,opt)

When I ran this code, the error message is followed:

Index in position 1 exceeds array bounds (must not exceed 1).

Error in
symengine>@(t,in2,param2)[sin(t.*1.842156345156404e3).*-1.0e2-in2(2,:).^2.*1.5e3-(7.0./5.0e1)./in2(1,:)+1.0./in2(1,:).^3.*1.13014e-1+1.0./param2.^4.*(in2(1,:).^3.*3.26e2-3.26e-4).*(in2(1,:)-9.999e3)-1.13000004e5;-in2(2,:)]

Error in Bubble_dynamics2>@(t,R)Fun(t,R)

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

Error in ode15s (line 150)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in Bubble_dynamics2 (line 34)
ode15s(H,[0,0.0001],R0,opt)

I have checked the size of H at line 34, but it is 1*1 sym, which does not exceed 1. I am very confused about the error, but I need to find a way to solve it.

Thank you!

[Yihan Zhang]

I amended the code so that it is much easier for you to understand what I am doing.

Answer 1

[Walter Roberson]

Use the flow shown in https://www.mathworks.com/help/symbolic/odefunction.html#bvlo6zg-1

Replace

[M,F]=massMatrixForm(eqs,vars)
M=odeFunction(M,vars)
Fun=odeFunction(F,vars,r)
H=@(t,R) Fun(t,R)

with

[M,F]=massMatrixForm(eqs,vars)
f = M\F;
H = odeFunction(f,vars)

However this will fail because you have the unresolved variable "r" that is distinct from R(t) . Also your equations involve R(t) and DR(t) so you need at least two boundary conditions but your R0 is only a single boundary condition.

[Yihan Zhang]

Thank you for your answer, but I have two boundary conditions in the code, one is for R(t) given by R0, and one is for DR(t) given by Rp0.

r is radial coordinate that I need to use it as a variable and itegrate (tau_rr-tau_tt)/r, to get an expression of R(t).

[Walter Roberson]

ode15s(H,[0,0.0001],R0,opt)

That passes in R0 as the boundary condition.

R0=100e-4; % m（boundary condition) R(t=0)

That is a scalar.

Perhaps you should be passing in [R0, Rp0]

[Walter Roberson]

r is radial coordinate that I need to use it as a variable and itegrate (tau_rr-tau_tt)/r, to get an expression of R(t).

That appears to refer to the sub-expression

int((tau_rr-tau_tt)/r,R(t),9999)

However, when you int(expression, a, b) then the integration is with respect to the default variable, chosen by symvar(expression,1)

symvar(s,n) chooses the n symbolic variables in s that are alphabetically closest to x and returns them in alphabetical order.

t is alphabetically closer to x than r is, so t is the default variable, so that becomes

int((tau_rr-tau_tt)/r,t,R(t),9999)

leaving r as an unresolved variable.

[Yihan Zhang]

I changed the

int((tau_rr-tau_tt)/r,R(t),9999)

to

int((tau_rr-tau_tt)/r,r,R(t),9999)

A new error message occurs, which I cannot even understand what is happening.

Error using symengine
Unable to generate code for piecewise for use in anonymous functions.

Error in sym/matlabFunction>mup2mat (line 404)
res = mupadmex('symobj::generateMATLAB',r.s,ano,spa,0);

Error in sym/matlabFunction>mup2matcell (line 374)
    r = mup2mat(c{1},true,sparseMat);

Error in sym/matlabFunction (line 188)
    body = mup2matcell(funs, opts.Sparse);

Error in sym/odeFunction (line 135)
  f = matlabFunction(expr, 'Vars', varsAndParams, matlabFunOptions{:});

Error in Bubble_dynamics2 (line 30)
M=odeFunction(M,vars)

[Walter Roberson]

Use the 'File' option of odeFunction -- code can be generated for piecewise if you do that.

Note that the piecewise that is generated tests whether R(t) <= 0 and says that the result of the integral is infinite in that case, which is a situation that would likely ruin your ode calculations. I would suggest that you look carefully at your equations to figure out whether it makes sense that R(t) could ever be non-positive: if it does then you have a problem with your system of equations. But if R(t) can never be non-positive then assumptions and simplifications can be added that have the effect of getting rid of the piecewise.

[Yihan Zhang]

OK, I used 'File' option of odeFunction as follow

M=odeFunction(M,vars,'File','myfileM','Comments','Version: 1.1')
Fun=odeFunction(F,vars,r,'File','myfileF','Comments','Version: 1.1')

The error message occurs again,

Index in position 1 exceeds array bounds (must not exceed 1).

Error in myfileF (line 9)
DRt = in2(2,:);

Error in Bubble_dynamics2>@(t,R)Fun(t,R)

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

Error in ode15s (line 150)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in Bubble_dynamics2 (line 34)
ode15s(H,[0,0.0001],R0,opt)

[Walter Roberson]

The below does not crash... but the output is not interesting.

%% EQUATIONS GOVERNING THE BUBBLE DYNAMICS
% Yihan Zhang
syms R(t)
syms r real
% parameters
tau_y = 8.6; % Pa
gamma = 0.07; % Pa m
rho = 1000; % kg/m^3
G = 81.5; % Pa
P0 = 1.13e5; % Pa
deltaP = 100; % Pa
eta_s = 0.001; % Pa s
R0 = 100e-4; % m（boundary condition) R(t=0)
Rp0 = 0; % m/s (boundary condition) dR/dt(t=0)
w = ((3*P0+4*gamma/R0+4*G)/(rho*R0^2))^0.5; % s^-1
% expressions of other variables
tau_rr = -4*G.*(R(t).^3-R0^3)./(3.*r.^3);
tau_tt = -tau_rr./2;
Pg = (P0+2*gamma/R0).*(R0./R(t)).^3;
P_inf = P0+deltaP.*sin(w.*t);
P = Pg+tau_rr-eta_s.*4.*R(t)./R(t)-2.*gamma./R(t);
% differential equation
eqn = 1000*(diff(R(t),t,2)*R(t)+3/2*diff(R(t),t)^2)==P-P_inf-tau_rr+2*int((tau_rr-tau_tt)/r,r,R(t),9999);
[eqs,vars] = reduceDifferentialOrder(eqn,R(t));
[M,F] = massMatrixForm(eqs,vars);
f = M\F;
H = odeFunction(f, vars, 'file', 'bubble_ode.m');
Mfun = odeFunction(M, vars, 'file', 'massfun.m');
opt = odeset('mass', Mfun, 'InitialSlope', Rp0);
ode15s(H, [0,0.0001], [R0, Rp0], opt)

[Yihan Zhang]

Thank you! I see what you mean. The result is wrong indeed. But thank you afterall.