Integer operands are required for colon operator when used as index.
6 views (last 30 days)
Show older comments
University Glasgow
on 23 Aug 2022
Answered: Simon Chan
on 23 Aug 2022
Warning: Integer operands are required for colon operator when used as index.
> In Case1_Ijuptilk_130822 (line 75)
Warning: Integer operands are required for colon operator when used as index.
> In Case1_Ijuptilk_130822 (line 83)
%% initialization
clear
clc
close all
%% Model parameters
global k1 eta1 alpha3 gamma1 d N Phi xi h A B C G
k1 = 6*10^(-12); % elastic constant
eta1 = 0.0240; % viscosity
xi = 0.585; % activity parameter
alpha3 = -0.001104; % viscosity
gamma1 = 0.1093; % viscosity
Theta = 0.0001;
d = 0.0002;
% Plotting r(q)
alpha=1-alpha3^2/gamma1/eta1;
q=0:0.01:(5*pi);
r=q-(1-alpha)*tan(q)+(alpha3*xi*alpha/eta1)./(4*k1*q.^2/d^2-alpha3*xi/eta1).*tan(q);
figure
plot(q,r)
axis([0 5*pi -20 20])
% Simplify model parameters
A = k1/gamma1;
B = alpha3/gamma1;
C = eta1 - alpha3*B;
G = alpha3*A;
N = 30; % N must be even
%% Numerical setup
% step size
h = d/N;
% t span
tspan = [0:0.1:100];
nt=length(tspan);
% range of z
%z=linspace(0,d/2,N+1)
z=linspace(0,d,N+1);
% initial conditions
theta0 = Theta*sin(pi*z/d);
v0 = zeros(1,N+1);
theta0_int=theta0(2:N);
v0_int=v0(2:N);
u0 = [theta0_int'; v0_int'];
% Constant mass matrix M: We use M to represent the left hand side of the system of equations.
M1=eye(N-1,N-1);
M2=zeros(N-1,N-1);
M=[M1 M2;M2 M2];
% Boundary Conditions
Phi = 0;
%% ode solver
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
[t,y] = ode15s(@lcode1,tspan,u0, options);
% Extract the solution for theta and v
theta = [Phi*ones(length(t), 1) y(:,1:N-1) Phi*ones(length(t), 1)];
v = [zeros(length(t), 1) y(:,N:(2*N-2)) zeros(length(t), 1)];
% theta at the middle of the layer (i.e., z=d/2)
theta_middle = theta(:, N/2);
%% Plot the solution.
figure
subplot(2,1,1)
plot(z, theta(1:(nt/10):nt,:))
xlabel('z(\mum)')
ylabel('\theta(z,t)(rad)')
legend('t_0 =0', 't_1=10', 't_2=20', 't_3=30', 't_4=40', 'Location','bestoutside')
%title('Director angle','FontSize',10)
%[xmin(z) xmax(z) ymin(theta) ymax(theta)];
subplot(2,1,2)
plot(z, v(1:(nt/10):nt,:))
xlabel('z(\mum)')
ylabel('v(z,t)(m/s)')
legend('t_0 =0', 't_1=10', 't_2=20', 't_3=30', 't_4=40', 'Location','bestoutside')
%hF=gcf;
%exportgraphics(hF,'thetav007neg.png','Resolution',300)
%title('Flow velocity')
%[xmin(z) xmax(z) ymin(theta) ymax(theta)];
figure
plot(t,theta_middle)
xlabel('t (s)')
ylabel('\theta(d/2,t)(rad)')
%title('Director angle at the middle of the layer', 'FontSize',6)
%hF=gcf;
%exportgraphics(hF,'thetamiddleneg.png','Resolution',300)
%% Functions
% Right hand side of the ODEs: F(U)
function rhsode = lcode1(t,y)
global Phi xi h A B C G N
% initialize theta and v
theta = y(1:(N-1));
v = y(N:(2*N-2));
% theta equations
% theta equations
% for the positon to the right of the left hand boundary z=0
rhsode(1,1) = (A/(h^2))*(theta(2)-2*theta(1)+ Phi) - (B/(2*h))*(v(2)-0);
% for all other internal positions 0<z<d
for i=2:(N-2)
rhsode(i,1) = (A/(h^2))*(theta(i+1)-2*theta(i)+theta(i-1)) - (B/(2*h))*(v(i+1)-v(i-1));
end
% for the positon to the left of the right hand boundary z=d
rhsode(N-1,1) = (A/(h^2))*(Phi-2*theta(N-1)+ theta(N-2)) - (B/(2*h))*(0-v(N-2));
% v equations [REMEMBER theta the RHS index for the v equations is N-1 more than the indices of the theta and v variables in this function]
% for the two positons to the right of the left hand boundary z=0
rhsode(N,1) = (G/(h^3))*(-Phi + 3*theta(1) - 3*theta(2) + theta(3)) + (C/(h^2))*(v(2) -2*v(1)+ 0) + (xi/(2*h))*(theta(2)-Phi);
rhsode(N+1,1) = (G/(2*h^3))*(Phi - 2*theta(1) + 2*theta(3)- theta(4)) + (C/(h^2))*(v(3) -2*v(2) + v(1)) + (xi/(2*h))*(theta(3)-theta(1));
% for all other internal positions 0<z<d
for i=(N+2):(2*N-4)
rhsode(i,1) = (G/(2*h^3))*(theta(i-2-(N-1)) - 2*theta(i-1-(N-1)) + 2*theta(i+1-(N-1))- theta(i+2-(N-1))) +(C/(h^2))*(v(i+1-(N-1)) -2*v(i-(N-1)) + v(i-1-(N-1))) + (xi/(2*h))*(theta(i+1-(N-1))-theta(i-1-(N-1)));
end
% for the two positons to the left of the right hand boundary z=d
rhsode(2*N-3,1) = (G/(2*h^3))*(theta(N-4) - 2*theta(N-3) + 2*theta(N-1) - Phi) +(C/(h^2))*(v(N-1) -2*v(N-2) + v(N-3)) + (xi/(2*h))*(theta(N-1)-theta(N-3));
rhsode(2*N-2,1) = (G/(h^3))*(-theta(N-3) + 3*theta(N-2) - 3*theta(N-1) + Phi) +(C/(h^2))*(0 -2*v(N-1) + v(N-2)) + (xi/(2*h))*(Phi-theta(N-2));
end
1 Comment
Dyuman Joshi
on 23 Aug 2022
tspan = [0:0.1:100];
nt=length(tspan)
%line 75
plot(z, theta(1:(nt/10):nt,:))
%line 83
plot(z, v(1:(nt/10):nt,:))
%What you are trying to do
1:(nt/10):nt
Array indices should be Positive Integers, which as you can see is not the case in here. It is because the value of nt is 1001 and not 1000 as you expected it to be.
%one modification can be
1:(nt-1)/10:nt
Even after this, your plot would give error because there will be a size mismatch
N=30;
d = 0.0002;
z=linspace(0,d,N+1);
size(z)
Accepted Answer
Simon Chan
on 23 Aug 2022
When you plot the solution, you use the following index:
1:(nt/10):nt
However, variable "nt" is equals to 1001 and the result is not an integer when divided by 10 (which is 100.1). So it gives the mentioned warning to you.
Just change the indexing to the following:
1:round(nt/10):nt
0 Comments
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!