Integer operands are required for colon operator when used as index.

9 views (last 30 days)
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
Dyuman Joshi on 23 Aug 2022
tspan = [0:0.1:100];
nt=length(tspan)
nt = 1001
%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
ans = 1×10
1.0000 101.1000 201.2000 301.3000 401.4000 501.5000 601.6000 701.7000 801.8000 901.9000
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
ans = 1×11
1 101 201 301 401 501 601 701 801 901 1001
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)
ans = 1×2
1 31

Sign in to comment.

Accepted Answer

Simon Chan
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

More Answers (0)

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by