MATLAB Answers

Plot error: index in position 1 is invalid

12 views (last 30 days)
Hai Nguyen
Hai Nguyen on 29 Mar 2021
Commented: KSSV on 30 Mar 2021
I want to plot Psi(0,t) and save it as a figure (line 285) in the long code below but got error: Index in position 1 is invalid. Array indices must be positive integers or logical values. How to fix this please? The code is also attached.
% colP3D
clear; clf;
% penetration flag 0 = PKN, 1 = P3D, dimensionless: toughness, leak-off;
ipen = 1; K = 0; C = 1;
% # of Colloc pts, Geometric time factor, Max Newton iters, residual Tol
N = 30; rfac = 1.2; Maxit = 30; tol = 1e-9;
% some constants
gamma_0 = 1.0006328466775270;
gamma_mt = 0;
t0 = 1e-6;
tM = 50;
tc = tM;
if C > 0
gamma_mt = 2/pi/C;
Cp = 10/3;
tc = (gamma_mt/gamma_0)^Cp;
t0 = 1e-6/C^Cp;
tM = 1e5/C^Cp;
end
%
lambda_u = (8+sqrt(K^4+64))/K^2; % runaway height growth boundary
% set up mesh and time steps
t = t0;
dt = t0/10;
dx = 1/N;
x = 0:dx:1;
xf = 0:dx/20:1;
xft = 1-dx:dx/50:1;
in = 1:N;
xi = x(in);
xh = 0.5*(x(1:N)+x(2:N+1));
[xph,is] = sort([x xh(1:N-1)]);
% precalculate PHI(OMEGA)
[PHI,LAMBDA,OMEGA] = PHIOMEGA(N,K,lambda_u,ipen);
% initialize the solution to the storage PKN values
[Omega0,P0,lambda0,Psi0,gamma0,gammadot] = init_PKN(x,t0); %#ok<ASGLU> % collocation pts
[Omega0h,P0h,lambda0h,Psi0h,gamma0,gammadot] = init_PKN(xh,t0); % 1/2 pts
% now initialize the collocation soln at channel points
y0 = [Omega0(in);
Psi0(in)];
y = y0;
y0h = [Omega0h(in);
Psi0h(in)];
lambda = ones(N+1,1);
itime = 0;
tauh = 0:t/10:t;
gammah = gamma_0*tauh.^(4/5);
gamma = gamma_0*(t+dt).^(4/5);
xi_l = 0;
fail = false;
while t < tM && lambda(1) < min(5,lambda_u)
t = t + dt;
if t >= tM; break; end
R = [];
tauh = [tauh t]; %#ok<AGROW>
gammah = [gammah gamma]; %#ok<AGROW>
for iter = 1:Maxit
b = getRHS(xi,y,gamma,K,t,tc,y0,y0h,gamma0,dt,OMEGA,PHI,C,tauh,gammah,ipen);
J = getJacN(xi,y,gamma,K,t,tc,y0,y0h,gamma0,dt,b,OMEGA,PHI,C,tauh,gammah,ipen);
dS = J\b;
y = y-reshape(dS(1:2*N),2,N);
gamma = gamma-dS(2*N+1);
gammah(end) = gamma;
Omega = [y(1,:) 0];
lambda = LamEval(OMEGA,LAMBDA,Omega);
if lambda(1) > lambda_u
disp(' lambda_u exceeded ');
break;
end
Psi = [y(2,:) 0];
if iter > 2 && norm(b,2) < tol; break; end
end
if norm(b,2) >= tol
fprintf('No convergence at t = %g out of %g\n', t, tM);
fail = true;
break;
end
[yh,Vol,Leak,Atau,alp] = getyh(xi,y,gamma,K,t,tc,y0,y0h,gamma0,dt,OMEGA,PHI,C,tauh,gammah,ipen);
lambda = LamEval(OMEGA,LAMBDA,Omega);
itime = itime+1;
TS(itime) = t; %#ok<SAGROW>
GAMS(itime) = gamma; %#ok<SAGROW>
ip = find(lambda > 1,1,'last');
if ip > 1
xi_l = min(x(ip-1)+(1-lambda(ip-1))/(lambda(ip)-lambda(ip-1))*(x(ip)-x(ip-1)),1);
end
XIL(itime) = xi_l; %#ok<SAGROW>
plotsol(N,x,Omega,xf,t,K,C,Psi,ipen,xft,Atau,alp,lambda,lambda_u,TS,GAMS,gamma_0,gamma_mt,XIL)
% copy vectors over for the next time step
y0 = y;
y0h = yh;
dgamma = gamma-gamma0;
gamma0 = gamma;
gamma = gamma0+dgamma;
dt = rfac*dt;
end
if fail
disp('Stopped early due to failure');
else
disp('All done!');
end
function [PHI,LAMBDA,OMEGA] = PHIOMEGA(N,K,lambda_u,ipen)
% set up form factor PHI a priori at sample points ready for spline interp
if ipen == 0
OMEGA = (1:N);
LAMBDA = ones(1,N);
PHI = ones(1,N);
else
Omegafun = @(x,k) (k*x.^(3/2)+2*sqrt(x.^2-1))/pi;
if K == 0
lambdaM = 10;
else
lambdaM = lambda_u;
end
Ns = 1500;
ns = 1:Ns;
a = 1.11;
ep = (exp(lambdaM/Ns/a)-1);
LAMBDA = [1:0.001:1.1 a*(1+ep).^ns];
OMEGA = Omegafun(LAMBDA,K);
PHI = Phi(LAMBDA,K);
end
end
function y = Phi(lambda,K)
y = ones(1,length(lambda));
i2 = find(lambda > 1);
lam = lambda(i2);
y(i2) = pi^2*(4*sqrt(lam)-K*sqrt(lam.^2-1)).*IntOmega3(lam,K)./(lam.^2.*(4*sqrt(lam)+3*K*sqrt(lam.^2-1)))./(((K*lam.^(3/2)+2*sqrt(lam.^2-1))/pi).^3)/12;
end
function y = IntOmega3(lambda,K)
y = zeros(1,length(lambda));
for i = 1:length(lambda)
lam = lambda(i);
y(i) = 2*quadgk(@(x)OmegaF3(x,lam,K),0,0.5*lam,'abstol',1e-12,'reltol',1e-8);
end
end
function Omega3 = OmegaF3(x,lambda,K)
Omega3 = ((4/pi)*(((K/pi/sqrt(lambda)-(2/pi)*asin(1/lambda))*sqrt(lambda^2-4*x.^2))+(2/pi)*(sqrt(lambda^2-4*x.^2)*asin(1/lambda)+phi12(x,lambda)))).^3;
end
function y = phi12(x,lambda)
i0 = abs(x) == 0.5;
i1 = find(abs(x) < 0.5);
i2 = find(0.5 < abs(x) & abs(x) <= 0.5*lambda);
y(i0) = log(lambda);
y(i1) = -2*x(i1).*atanh((2*x(i1)*sqrt(lambda^2-1))./sqrt(lambda^2-4*x(i1).^2))+atanh((sqrt(lambda^2-1))./sqrt(lambda^2-4*x(i1).^2));
y(i2) = -2*x(i2).*atanh(sqrt(lambda^2-4*x(i2).^2)./(2*x(i2)*sqrt(lambda^2-1)))+atanh(sqrt(lambda^2-4*x(i2).^2)./sqrt(lambda^2-1));
end
function [Omega0,P0,lambda0,Psi0,gamma0,gammadot] = init_PKN(x,t0)
% initialize the solution at collocation points to the impermeable PKN
gamma_0 = 1.0006328466775270;
Omega_0 = ((12/5)*gamma_0^2)^(1/3);
N = length(x);
gamma0 = gamma_0*t0^(4/5);
gammadot = (4/5)*gamma_0*t0^(-1/5);
Omega0 = Omega_0*t0^(1/5)*(1-x).^(1/3).*(1-(1-x)/96 +23*(1-x).^2/64512);
P0 = Omega0;
lambda0 = ones(N,1);
Psi0 = t0^(1/9)*(1-x).^(1/3);
end
function b = getRHS(x,y,gamma,K,t,tc,y0,y0h,gamma0,dt,OMEGA,PHI,C,tauh,gammah,ipen)
[n,N] = size(y); %#ok<ASGLU>
F = FS(x,y,gamma,K,y0,gamma0,dt,OMEGA,PHI,C,tauh,gammah,ipen);
h = diff(x);
H = spdiags(h',0,N-1,N-1);
ii = 1:N-1;
ip1 = 2:N;
xi = x(ii);
yi = y(:,ii);
Fi = F(:,ii);
xip1 = x(ip1);
yip1 = y(:,ip1);
Fip1 = F(:,ip1);
xih = 0.5*(xi+xip1);
yih = 0.5*(yi+yip1)-(Fip1-Fi)*H/8;
Fih = FS(xih,yih,gamma,K,y0h(:,ii),gamma0,dt,OMEGA,PHI,C,tauh,gammah,ipen);
Phim = yip1-yi-(Fip1+4*Fih+Fi)*H/6;
dgamma = (gamma-gamma0)/dt; %
if t < tc
alp = 1/3;
Atau = (3*dgamma*gamma)^(1/3);
else
alp = 3/8;
Atau = 2*(C/3)^(1/4)*gamma^(3/8)*dgamma^(1/8);
end
Voltip = Atau*h(1)^(1+alp)/(1+alp);
tipVal = Atau*h(1)^alp;
Vol = (yip1(1,:)+4*yih(1,:)+yi(1,:))*h'/6+Voltip;
Leaktip = (2/3)*(gamma*dt/(gamma-gamma0))^(1/2)*h(1)^(3/2);
Leak = (sqrt(t-spline(gammah,tauh,gamma*xip1))+4*sqrt(t-spline(gammah,tauh,gamma*xih))+sqrt(t-spline(gammah,tauh,gamma*xi)))*h'/6+Leaktip;
Phig = gamma*(Vol+2*C*Leak) - t;
Phib = GS(y(:,1), y(:,N), tipVal,t);
b = [Phib(:);
Phim(:);
Phig];
end
function F = FS(x,y,gamma,K,y0,gamma0,dt,OMEGA,PHI,C,tauh,gammah,ipen) %#ok<INUSL>
Ups = Upsilon(OMEGA,PHI,y(1,:),ipen);
tau = tauh(end);
F(1,:) = -gamma*y(2,:)./y(1,:).^3./Ups;
tau0 = spline(gammah,tauh,gamma*x);
F(2,:) = -gamma*(y(1,:)-y0(1,:))/dt-(gamma-gamma0)*gamma*x.*y(2,:)./y(1,:).^3./Ups/dt-C*gamma./sqrt(tau-tau0);
end
function Ups = Upsilon(OMEGA,PHI,Omega,ipen)
if ipen == 0
Ups = ones(1,length(Omega));
else
i1 = find(Omega <= OMEGA(1));
i2 = find(Omega > OMEGA(1));
Ups(i1) = ones(1,length(i1));
Ups(i2) = spline(OMEGA,PHI,Omega(i2));
end %
end
function r = GS(ya, yb, tipVal,t)
alpha=1/9;
r(1) = ya(2)-t.^alpha;
r(2) = yb(1) - tipVal;
end
function J = getJacN(x,y,gamma,K,t,tc,y0,y0h,gamma0,dt,b,OMEGA,PHI,C,tauh,gammah,ipen)
ep = sqrt(eps);
[n,N] = size(y);
I = eye(n*N);
NC = 1:(n*N+1);
J = zeros(n*N+1, n*N);
for j = 1:n*N
yt = reshape(y(:)+ep*I(:,j),n,N);
J(NC,j) = (getRHS(x,yt,gamma,K,t,tc,y0,y0h,gamma0,dt,OMEGA,PHI,C,tauh,gammah,ipen)-b)/ep;
end
J(NC,n*N+1) = (getRHS(x,y,gamma+ep,K,t,tc,y0,y0h,gamma0,dt,OMEGA,PHI,C,tauh,[gammah(1:end-1) gammah(end)+ep],ipen)-b)/ep;
end
function lambda = LamEval(OMEGA,LAMBDA,Omega)
i1 = find(Omega <= OMEGA(1));
i2 = find(Omega > OMEGA(1));
lambda(i1) = ones(1,length(i1));
lambda(i2) = spline(OMEGA,LAMBDA,Omega(i2));
end
function [yih,Vol,Leak,Atau,alp] = getyh(x,y,gamma,K,t,tc,y0,y0h,gamma0,dt,OMEGA,PHI,C,tauh,gammah,ipen) %#ok<INUSL>
[n,N] = size(y); %#ok<ASGLU>
F = FS(x,y,gamma,K,y0,gamma0,dt,OMEGA,PHI,C,tauh,gammah,ipen);
h = diff(x);
H = spdiags(h',0,N-1,N-1);
ii = 1:N-1;
ip1 = 2:N;
xi = x(ii);
yi = y(:,ii);
Fi = F(:,ii);
xip1 = x(ip1);
yip1 = y(:,ip1);
Fip1 = F(:,ip1);
xih = 0.5*(xi+xip1);
yih = 0.5*(yi+yip1)-(Fip1-Fi)*H/8;
dgamma = (gamma-gamma0)/dt;
if t < tc
alp = 1/3;
Atau = (3*dgamma*gamma)^(1/3);
else
alp = 3/8;
Atau = 2*(C/3)^(1/4)*gamma^(3/8)*dgamma^(1/8);
end
Voltip = Atau*h(1)^(1+alp)/(1+alp);
tipVal = Atau*h(1)^alp; %#ok<NASGU>
Vol = (yip1(1,:)+4*yih(1,:)+yi(1,:))*h'/6+Voltip;
Leaktip = (2/3)*(gamma*dt/(gamma-gamma0))^(1/2)*h(1)^(3/2);
Leak = (sqrt(t-spline(gammah,tauh,gamma*xip1))+4*sqrt(t-spline(gammah,tauh,gamma*xih))+sqrt(t-spline(gammah,tauh,gamma*xi)))*h'/6+Leaktip;
end
function plotsol(N,x,Omega,xf,t,K,C,Psi,ipen,xft,Atau,alp,lambda,lambda_u,TS,GAMS,gamma_0,gamma_mt,XIL)
figure(4);
plot(t,Psi(0,t),'bo--','linewidth',2,'markerfacecolor','b')
xlabel('\tau');
ylabel('\Psi(0,\tau) ');
savefig('case2rate');
figure(5);
plot(x,lambda,'ro-',x,ones(1,N+1)*lambda_u,'k','linewidth',2,'markerfacecolor','r','markersize',6)
xlabel('\xi');
ylabel(' \lambda(\xi) ');
lamx =[' \lambda_u = ',num2str(lambda_u)];
legend(' Classical P3D',lamx)
if lambda(1)<lambda_u
ax=axis;ax(3)=1;
ax(4)=1.2*lambda(1);
axis(ax);
end
pause(0.1)
end
function Omega = exact_PKN(x,t)
gamma_0 = 1.0006328466775270;
Omega_0 = ((12/5)*gamma_0^2)^(1/3);
Omega = Omega_0*t^(1/5)*(1-x).^(1/3).*(1-(1-x)/96 +23*(1-x).^2/64512);
end
function Omega = exact_PKNC(x,t,C)
gamma_m0 = 2/pi/C; %#ok<NASGU>
Omega_m0 = 2^(11/8)/pi^(1/2)/(3*C)^(1/4);
Omega = Omega_m0*t^(1/8)*(1-x).^(3/8).*(1+(1-x)/80);
end

Accepted Answer

KSSV
KSSV on 29 Mar 2021
This line:
plot(t,Psi(0,t),'bo--','linewidth',2,'markerfacecolor','b')
t is a scalar value first i.e. t = 1.1000e-06; and Psi is an array of size 1X31....it looks like you are trying to access some value from Psi using: Psi(0,t), this is not correct; index cannot be zero, negative and fraction in MATLAB. You need to rethink on your code. Your code is buggy and got lot ov variables which are not being used.
Learn about debugging this will help you a a lot.
  7 Comments
KSSV
KSSV on 30 Mar 2021
It seems you are plotting point by point, like shown below:
figure
hold on
for i = 1:10
plot(i,rand,'.')
end
To join them, don't plot point by point. Save the points into an array and then plot them at once. Like below:
a = zeros(10,1) ;
for i = 1:10
a(i) = rand ;
end
plot(1:10,a)

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by