How to fix the error in integration?

Hi,

While performing an integration using trapz it throws an error "Array indices must be positive integers or logical values."

I'm unable to find the cause for this error.

Your advice will be highly appreciated.

pde2fshear_v4Perturbed_nonlinear()

Warning: The value of local variables may have been changed to match the globals. Future versions of MATLAB will require that you declare a variable to be global before you use that variable.

Warning: The value of local variables may have been changed to match the globals. Future versions of MATLAB will require that you declare a variable to be global before you use that variable.

Array indices must be positive integers or logical values.

Error in solution>pdex2pde (line 463)
I_p=trapz(X(1:n),jb(X(1:n))+jd(X(1:n)));

Error in pdepe (line 242)
[c,f,s] = feval(pde,xi(1),t(1),U,Ux,varargin{:});

Error in solution>pde2fshear_v4Perturbed_nonlinear (line 143)
sol = pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,x,t);

function pde2fshear_v4Perturbed_nonlinear
global x;
global X;
global t;
global chi0;   % declare global variables
global D0;
global chi1;
global D1;
global alpha_chi;
global alpha_D;
global H0;
global S0;
global H;
global S;
global data;
global track;
global xstep;
global tstep;
global count;
global chi_growth;
global lambda_suppress;
global v_e;
global chi_ano;
global D_ano;
global gamma_nonlin;
global alpha_nonlin;
global group_vel;
global drift_vel;
global intensity_diff;
global drift_velFluct;
global D_0;
global safetyFact;
global grad_safetyFact;
global elec_diaVel;
global grad_length;
global mag_shear;
global theta_deg;
global elec_temp;
global T_magFld;
global P_magFld;
global R0;
global I_p;
global jb;
global jd;
global B0;
global xmax;
global xmin;
global STEP_x;




data.constant.critgradpressure =  1.2; 
data.constant.critgraddensity =   1.0;
count = 1;

chi0 = 0.5;
D0 = 0.5;
chi1 = 5.0;
D1 = 5.0;
alpha_chi = 0.1;
alpha_D = 0.1;
H0 = 27;
S0 = 21;
track = 1;
chi_growth=20;
lambda_suppress=0.05;
chi_ano=10;
D_ano=10;




gamma_nonlin = 50;
alpha_nonlin = 1;
D_0 = 10;
STEP_x=1;





xstep = 100;
tstep = 1000; 
xmin = 0;
xmax = 1;
tmin = 0;
tmax = 50;
%tmax = 1;
m = 0; %define type of equation to solve


grad_u1 = zeros(tstep,xstep);
grad_u2 = zeros(tstep,xstep);
grad_u3 = zeros(tstep,xstep);
curve_u1 = zeros(tstep,xstep);
curve_u2 = zeros(tstep,xstep);
curve_u3 = zeros(tstep,xstep);
flowshear_p = zeros(tstep,xstep);
flowshear_n = zeros(tstep,xstep);
Q = zeros(tstep,xstep);
Q0 = zeros(tstep,xstep);
Q1 = zeros(tstep,xstep);
neo_p = zeros(tstep,xstep);
ano_p = zeros(tstep,xstep);
Gam = zeros(tstep,xstep);
Gam0 = zeros(tstep,xstep);
Gam1 = zeros(tstep,xstep);
neo_n = zeros(tstep,xstep);
ano_n = zeros(tstep,xstep);
safetyFact=zeros(tstep,xstep);
elec_diaVel=zeros(tstep,xstep);
grad_length=zeros(tstep,xstep);
mag_shear=zeros(tstep,xstep);
theta_deg=zeros(tstep,xstep);
elec_temp=zeros(tstep,xstep);
T_magFld=zeros(tstep,xstep);
P_magFld=zeros(tstep,xstep);
R0=zeros(tstep,xstep);
I_p=zeros(tstep,xstep);
jb=zeros(tstep,xstep);
B0=zeros(tstep,xstep);






x = linspace(xmin,xmax,xstep);
X = linspace(xmin,xmax,xstep);
t = linspace(tmin,tmax,tstep);



%for i=(xstep/5)+1:xstep
 %   x(i) = x(i-1)+(xmax/2)/(xstep-(xstep/5));
%end

data.variable.x = x;

sol = pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,x,t);
% Extract the first solution component as u1 = pressure
% second solution component as u2 = density
% third solution component as u3 = turbulence intensity



u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3); 


%grad_u = gradient(u,(x(2)-x(1)));
for j=1:tstep
    for i = 1:xstep/5
        if i == 1
            grad_u1(j,i) = (u1(j,2)-u1(j,1))/(x(2)-x(1));
            grad_u2(j,i) = (u2(j,2)-u2(j,1))/(x(2)-x(1));
            grad_u3(j,i) = (u3(j,2)-u3(j,1))/(x(2)-x(1));
           % grad_u4(j,i) = (u4(j,2)-u4(j,1))/(x(2)-x(1));
                                      
        elseif i == xstep/5
            grad_u1(j,i) = (u1(j,i)-u1(j,i-1))/(x(i)-x(i-1));
            grad_u2(j,i) = (u2(j,i)-u2(j,i-1))/(x(i)-x(i-1));
            grad_u3(j,i) = (u3(j,i)-u3(j,i-1))/(x(i)-x(i-1));
           % grad_u4(j,i) = (u4(j,i)-u4(j,i-1))/(x(i)-x(i-1));                       
        else
            grad_u1(j,i) = (u1(j,i+1)-u1(j,i-1))/(x(i+1)-x(i-1));
            grad_u2(j,i) = (u2(j,i+1)-u2(j,i-1))/(x(i+1)-x(i-1));
            grad_u3(j,i) = (u3(j,i+1)-u3(j,i-1))/(x(i+1)-x(i-1));  
           % grad_u4(j,i) = (u4(j,i+1)-u4(j,i-1))/(x(i+1)-x(i-1));
        end
    end

    for i=(xstep/5)+1:xstep
        if i == xstep/5+1
            grad_u1(j,i) = (u1(j,i+1)-u1(j,i))/(x(i+1)-x(i));
            grad_u2(j,i) = (u2(j,i+1)-u2(j,i))/(x(i+1)-x(i));
            grad_u3(j,i) = (u3(j,i+1)-u3(j,i))/(x(i+1)-x(i));      
           % grad_u4(j,i) = (u4(j,i+1)-u4(j,i))/(x(i+1)-x(i));
        elseif i == xstep
            grad_u1(j,i) = (u1(j,i)-u1(j,i-1))/(x(i)-x(i-1));
            grad_u2(j,i) = (u2(j,i)-u2(j,i-1))/(x(i)-x(i-1));
            grad_u3(j,i) = (u3(j,i)-u3(j,i-1))/(x(i)-x(i-1));   
           % grad_u4(j,i) = (u4(j,i)-u4(j,i-1))/(x(i)-x(i-1));
        else
            grad_u1(j,i) = (u1(j,i+1)-u1(j,i-1))/(x(i+1)-x(i-1));
            grad_u2(j,i) = (u2(j,i+1)-u2(j,i-1))/(x(i+1)-x(i-1));
            grad_u3(j,i) = (u3(j,i+1)-u3(j,i-1))/(x(i+1)-x(i-1));    
           % grad_u4(j,i) = (u4(j,i+1)-u4(j,i-1))/(x(i+1)-x(i-1));
        end
    end
end

 for i=1:tstep
     for j=1:xstep
          v_e = -grad_u1(i,j)*grad_u2(i,j)/u2(i,j)^2;             % -g_p*g_n/n^2
          flowshear_p(i,j) = 1+ alpha_chi*v_e^2;
          flowshear_n(i,j) = 1+ alpha_D*v_e^2;
         if abs(grad_u1(i,j)) < abs(data.constant.critgradpressure) %&& abs(grad_u2(i,j)) < abs(data.constant.critgraddensity) 
             Q(i,j) = -grad_u1(i,j)*chi0;
             Q0(i,j) = Q(i,j);
             Q1(i,j) = 0;
             neo_p(i,j) = chi0*(1+grad_u1(i,j))/(1+grad_u1(i,j));
             ano_p(i,j) = 0;
             Gam(i,j) = -grad_u2(i,j)*D0;
             Gam0(i,j) = Gam(i,j);
             Gam1(i,j) = 0;
             neo_n(i,j) = D0*(1+grad_u2(i,j))/(1+grad_u2(i,j));
             ano_n(i,j) = 0;
             intensity_diff(i,j) = D_0*u3(i,j).^alpha_nonlin;
             drift_velFluct(i,j)=D_0*alpha_nonlin*grad_u3(i,j).^(alpha_nonlin-1)*grad_u3(i,j);
             drift_vel(i,j) = group_vel + drift_velFluct(i,j);
            
         elseif abs(grad_u1(i,j)) >= abs(data.constant.critgradpressure) %&& abs(grad_u2(i,j)) < abs(data.constant.critgraddensity)
            Q(i,j) = (chi0*(-grad_u1(i,j)) + chi_ano*(-grad_u1(i,j)+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j))*(-grad_u1(i,j)); 
            Q0(i,j) = -grad_u1(i,j)*chi0;
            Q1(i,j) = (chi_ano*(-grad_u1(i,j)+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j))*(-grad_u1(i,j));
            neo_p(i,j) = chi0*(1+grad_u1(i,j))/(1+grad_u1(i,j));
            ano_p(i,j) = chi_ano*(abs(grad_u1(i,j))+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j)*(-grad_u1(i,j)); 
            Gam(i,j) = -grad_u2(i,j)*D0;
            Gam0(i,j) = Gam(i,j);
            Gam1(i,j) = 0;
            neo_n(i,j) = D0*(1+grad_u2(i,j))/(1+grad_u2(i,j));
            ano_n(i,j) = 0;
            intensity_diff(i,j) = D_0*u3(i,j).^alpha_nonlin;
            drift_velFluct(i,j)=D_0*alpha_nonlin*grad_u3(i,j).^(alpha_nonlin-1)*grad_u3(i,j);
            drift_vel(i,j) = group_vel + drift_velFluct(i,j);
           
       else
             Q(i,j) = (chi0*(-grad_u1(i,j)) + chi_ano*(-grad_u1(i,j)+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j))*(-grad_u1(i,j)); 
             Q0(i,j) = -grad_u1(i,j)*chi0;
             Q1(i,j) = (chi_ano*(-grad_u1(i,j)+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j))*(-grad_u1(i,j));
             neo_p(i,j) = chi0*(1+grad_u1(i,j))/(1+grad_u1(i,j));
             ano_p(i,j) = chi_ano*(abs(grad_u1(i,j))+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j)*(-grad_u1(i,j));
             Gam(i,j) = -grad_u2(i,j)*(D0 + D_ano*(-grad_u2(i,j)+data.constant.critgraddensity)*u3(i,j)/flowshear_n(i,j));
             Gam0(i,j) = -grad_u2(i,j)*D0;
             Gam1(i,j) = -grad_u2(i,j)*D_ano*(-grad_u2(i,j)+data.constant.critgraddensity)*u3(i,j)/flowshear_n(i,j);
             neo_n(i,j) = D0*(1+grad_u2(i,j))/(1+grad_u2(i,j));
             ano_n(i,j) = D_ano*(abs(grad_u2(i,j))+data.constant.critgraddensity)/flowshear_n(i,j);
             intensity_diff(i,j) = D_0*u3(i,j).^alpha_nonlin;
             drift_velFluct(i,j)=D_0*alpha_nonlin*grad_u3(i,j).^(alpha_nonlin-1)*grad_u3(i,j);
             drift_vel(i,j) = group_vel + drift_velFluct(i,j);
            
         end
     end
 end


%curve_u = gradient(grad_u,(x(2)-x(1)));
for j=1:tstep
    for i = 1:xstep/5
        if i == 1
            curve_u1(j,i) = (grad_u1(j,2)-grad_u1(j,1))/(x(2)-x(1));
            curve_u2(j,i) = (grad_u2(j,2)-grad_u2(j,1))/(x(2)-x(1));
            curve_u3(j,i) = (grad_u3(j,2)-grad_u3(j,1))/(x(2)-x(1));
           % curve_u4(j,i) = (grad_u4(j,2)-grad_u4(j,1))/(x(2)-x(1));

        elseif i == xstep/5
            curve_u1(j,i) = (grad_u1(j,i)-grad_u1(j,i-1))/(x(i)-x(i-1));
            curve_u2(j,i) = (grad_u2(j,i)-grad_u2(j,i-1))/(x(i)-x(i-1));
            curve_u3(j,i) = (grad_u3(j,i)-grad_u3(j,i-1))/(x(i)-x(i-1));
           % curve_u4(j,i) = (grad_u4(j,i)-grad_u4(j,i-1))/(x(i)-x(i-1));
                                    
        else
            curve_u1(j,i) = (grad_u1(j,i+1)-grad_u1(j,i-1))/(x(i+1)-x(i-1));
            curve_u2(j,i) = (grad_u2(j,i+1)-grad_u2(j,i-1))/(x(i+1)-x(i-1));
            curve_u3(j,i) = (grad_u3(j,i+1)-grad_u3(j,i-1))/(x(i+1)-x(i-1));   
           % curve_u4(j,i) = (grad_u4(j,i+1)-grad_u4(j,i-1))/(x(i+1)-x(i-1));
        end
    end

    for i=(xstep/5)+1:xstep
        if i == xstep/5+1
            curve_u1(j,i) = (grad_u1(j,i+1)-grad_u1(j,i))/(x(i+1)-x(i));
            curve_u2(j,i) = (grad_u2(j,i+1)-grad_u2(j,i))/(x(i+1)-x(i));
            curve_u3(j,i) = (grad_u3(j,i+1)-grad_u3(j,i))/(x(i+1)-x(i));  
           % curve_u4(j,i) = (grad_u4(j,i+1)-grad_u4(j,i))/(x(i+1)-x(i));
        elseif i == xstep
            curve_u1(j,i) = (grad_u1(j,i)-grad_u1(j,i-1))/(x(i)-x(i-1));
            curve_u2(j,i) = (grad_u2(j,i)-grad_u2(j,i-1))/(x(i)-x(i-1));
            curve_u3(j,i) = (grad_u3(j,i)-grad_u3(j,i-1))/(x(i)-x(i-1));  
           % curve_u4(j,i) = (grad_u4(j,i)-grad_u4(j,i-1))/(x(i)-x(i-1));
        else
            curve_u1(j,i) = (grad_u1(j,i+1)-grad_u1(j,i-1))/(x(i+1)-x(i-1));
            curve_u2(j,i) = (grad_u2(j,i+1)-grad_u2(j,i-1))/(x(i+1)-x(i-1));
            curve_u3(j,i) = (grad_u3(j,i+1)-grad_u3(j,i-1))/(x(i+1)-x(i-1));   
           % curve_u4(j,i) = (grad_u4(j,i+1)-grad_u4(j,i-1))/(x(i+1)-x(i-1));
        end
    end
end

%to save parameters and variables
data.constant.chi0 = chi0;
data.constant.D0 = D0;
data.constant.chi1 = chi1;
data.constant.D1 = D1;
data.constant.alphachi = alpha_chi;
data.constant.alphaD = alpha_D;
data.constant.H0 = H0;
data.constant.S0 = S0;
data.variable.pressure = u1;
data.variable.density = u2;
data.variable.turbulence= u3;
data.variable.intensity_diff=intensity_diff;
data.variable.drift_velFluct=drift_velFluct;
data.variable.drift_vel=drift_vel;
data.variable.gradpressure = grad_u1;
data.variable.graddensity = grad_u2;
data.variable.gradintensity = grad_u3;
data.variable.curvepressure = curve_u1;
data.variable.curvedensity = curve_u2;
data.variable.curveturbulence = curve_u3;
data.variable.x = x;
data.variable.t = t;
data.variable.Q = Q;
data.variable.Gamma = Gam;
data.variable.Q0 = Q0;
data.variable.Gamma0 = Gam0;
data.variable.neo_P = neo_p;
data.variable.neo_n = neo_n;
data.variable.Q1 = Q1;
data.variable.Gam1 = Gam1;
data.variable.ano_p = ano_p;
data.variable.ano_n = ano_n;
data.variable.heatsource = H;
data.variable.particlesource = S;
data.variable.wexb_p = flowshear_p;
data.variable.wexb_n = flowshear_n;
data.control.xgrid = xstep;
data.control.tgrid = tstep;
data.control.xmin = xmin;
data.control.xmax = xmax;
data.control.tmin = tmin;
data.control.tmax = tmax;

figure;
subplot(2,1,1,'FontSize',22)
    plot(x,u1(end,:),'.');
    axis ([0 1 0 20]);
    ylabel('p')
    xlabel('r/a')
    %set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
    head = strcat('Time = ',num2str(t(end)),' s');
    title(head);
    grid on
    
    
    subplot(2,1,2,'FontSize',22)
    plot(x,grad_u1(end,:),'.');
    axis([0 1 -60 0]);
    ylabel('p\prime')
    xlabel('r/a')
    grid on
    
    


    
    

    
  
   

end   
 % --------------------------------------------------------------
function [c,f,s] = pdex2pde(x,t,u,DuDx)
global x;
global t;
global chi0;
global D0;
global chi1;
global D1;
global H0;
global S0;
global data;
global xstep;
global tstep;
global alpha_chi;
global alpha_D;
global chi_growth;         % total growth rate
global length;
global theta_heaviside1;
global lambda_suppress;
global v_e;
global gamma_nonlin;
global alpha_nonlin;
global group_vel;
global drift_vel;
global intensity_diff;
global drift_velFluct;
global D_0;
global safetyFact;
global grad_safetyFact;
global elec_diaVel;
global grad_length;
global mag_shear;
global theta_deg;
global elec_temp;
global T_magFld;
global P_magFld;
global R0;
global I_p;
global jb;
global jd;
global B0;
global X;
global xmax;
global xmin;
global STEP_x;


%lf FFf F  Fb vfDdength = 0.01 or 1
x_count=1;
length=1;
jb0=1;
jd0=1;
x0=1;
grad_length=DuDx(2);
R0=linspace(1,5,100);
B0=1;


c = [1;1;1];
v_e = -DuDx(1)*DuDx(2)/u(2)^2;             % -g_p*g_n/n^2
flowshear_p = 1+ alpha_chi*v_e^2;
flowshear_n = 1+ alpha_D*v_e^2;
%disp(X)
term1 = abs(DuDx(1))-data.constant.critgradpressure;
drift_velFluct = D_0*DuDx(3)^alpha_nonlin;
intensity_diff = D_0*u(3)^alpha_nonlin;
X = linspace(xmin,xmax,xstep-1);

% Implementing Heaviside function for H-mode in p, n and I equations
       
if term1 > 0
    theta_heaviside1=1;
else
    theta_heaviside1=0;
end

% Determining group velocity of turbulence intensity
if STEP_x > xstep-1
   STEP_x=1;
end

if STEP_x <= xstep-1
    
    X(STEP_x)=x(STEP_x);
   
end
%disp(X(STEP_x))
        num=15;
    if STEP_x==1
        jb=jb0*DuDx(1);
        jd=jd0.*(1-(X(STEP_x)-x0).^2).^num;
        %disp(jb)
        %disp(jd)
        for n=2:5
        I_p=trapz(X(1:n),jb(X(1:n))+jd(X(1:n)));
        end
    end   
        T_magFld=B0./R0;
        P_magFld=B0.*I_p./x;    
        safetyFact=(x./R0).*T_magFld./P_magFld;
        if STEP_x == 1
            grad_safetyFact(STEP_x) = (safetyFact(2)-safetyFact(1))/(x(2)-x(1));
        
        elseif STEP_x == xstep/5
            grad_safetyFact(STEP_x) = (safetyFact(STEP_x)-safetyFact(STEP_x-1))/(x(STEP_x)-x(STEP_x-1));

        else
            grad_safetyFact(STEP_x) = (safetyFact(STEP_x+1)-safetyFact(STEP_x-1))/(x(STEP_x+1)-x(STEP_x-1));
        end
  
            mag_shear=(x/safetyFact)*grad_safetyFact;
           
             
elec_diaVel = -elec_temp/1.6*10^(-19)*T_magFld*grad_length;
group_vel = - (elec_diaVel*(2*grad_length/mag_shear*R0)*sin(theta_deg));
%Turbulence intensity Equation for nonlinear turbulence intensity

s3 = (chi_growth*(term1*theta_heaviside1-lambda_suppress*v_e^2)-gamma_nonlin*u(3)^alpha_nonlin)*u(3)+group_vel*u(3)-(D_0*(1+alpha_nonlin)*alpha_nonlin*u(3)^(alpha_nonlin-1)*u(3)^2); 

 
s = [(H0)*exp(-100*x^2/length)+H0/2; (S0)*exp(-100*(x-0.9)^2/length)+S0/2;s3];


f = [chi0+chi1*u(3)/flowshear_p ; D0+D1*u(3)/flowshear_n ;-(intensity_diff-D_0*(1+alpha_nonlin)*u(3)^alpha_nonlin)].*DuDx;                 % flux term for nonlinear model

 
    
end

% --------------------------------------------------------------
function u0 = pdex2ic(x)
  %u0 = [eps; eps; eps];
 %u0 = [0.01; 0.01;0.1*exp(-100*(x-1)^2)];
  u0 = [0.1*(1-x^2); 0.1*(1-x^2); 0.5*exp(-100*(x-1)^2)];
end
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t)

pl = [0; 0; 0];
ql = [1; 1; 1];


pr = [ur(1); ur(2);0];
qr = [0.01; 0.1; 1];


%---------------------------------------------------------------
end

Answer 1

[Walter Roberson]

xmin = 0;
xmax = 1;

jb=zeros(tstep,xstep);

X = linspace(xmin,xmax,xstep-1);

The X values are going to be between 0 and 1.

 I_p=trapz(X(1:n),jb(X(1:n))+jd(X(1:n)));

You are indexing the 2D array jb at a vector, which is not necessarily impossible but should be a sign of caution.

You are indexing it at X values. X values are between 0 and 1, and so are not integers.

[Rahul]

Hi Walter,

I cannot increase xmax more than 1 as this is my max system length.

Can you advice me how this integration logic can be corrected? I'm not having any clue.

Rgds,

rc

[Walter Roberson]

 I_p = trapz(X(1:n),jb(1:n)+jd(1:n));