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

1 visualización (últimos 30 días)
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!

Respuesta aceptada

Walter Roberson
Walter Roberson el 16 de Ag. de 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 comentarios
Walter Roberson
Walter Roberson el 17 de Ag. de 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)
Yihan Zhang
Yihan Zhang el 17 de Ag. de 2019
Thank you! I see what you mean. The result is wrong indeed. But thank you afterall.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Programming en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by