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

1 view (last 30 days)
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!

Accepted Answer

Walter Roberson
Walter Roberson on 16 Aug 2019
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.
  8 Comments
Walter Roberson
Walter Roberson on 17 Aug 2019
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)

Sign in to comment.

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!